boost::ublas how to get determinant of int matrix?

Question

I found the function that calculates the determinant of boost::ublas matrix:

template<typename ValType>
ValType det_fast(const ublas::matrix<ValType>& matrix)
{
    // create a working copy of the input 
    ublas::matrix<ValType> mLu(matrix);
    ublas::permutation_matrix<std::size_t> pivots(matrix.size1());

    auto isSingular = ublas::lu_factorize(mLu, pivots);
    if (isSingular)
        return static_cast<ValType>(0);

    ValType det = static_cast<ValType>(1);
    for (std::size_t i = 0; i < pivots.size(); ++i) 
    {
        if (pivots(i) != i)
            det *= static_cast<ValType>(-1);

        det *= mLu(i, i);
    }

    return det;
}

This function works fine, but only with non-integer types (it works fine with float and double). When I try to pass the same matrix but with type of int, I received compilation error:

Check failed in file c:\boost\boost_1_58_0\boost\numeric\ublas\lu.hpp at line 167: singular != 0 || detail::expression_type_check (prod (triangular_adaptor (m), triangular_adaptor (m)), cm) unknown location(0): fatal error in "BaseTest": std::logic_error: internal logic

It is the bug of boost or my function is wrong? What change I could do for avoid this error?

Nikolay K · Accepted Answer · 2015-04-28T15:49:43.933

Of course this is not a bug, since checks like this one is all over uBLAS. I guess this is because most of its operations would make no sense for non-float types.

You can disable type checks using

#define BOOST_UBLAS_TYPE_CHECK 0

before includes. But think twice! Take a look at example

#include <iostream>
#define BOOST_UBLAS_TYPE_CHECK 0
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/lu.hpp>
namespace ublas = boost::numeric::ublas;

template<typename ValType>
ValType det_fast(const ublas::matrix<ValType>& matrix)
{
    // create a working copy of the input 
    ublas::matrix<ValType> mLu(matrix);
    ublas::permutation_matrix<std::size_t> pivots(matrix.size1());

    auto isSingular = ublas::lu_factorize(mLu, pivots);
    if (isSingular)
        return static_cast<ValType>(0);

    ValType det = static_cast<ValType>(1);
    for (std::size_t i = 0; i < pivots.size(); ++i) 
    {
        if (pivots(i) != i)
            det *= static_cast<ValType>(-1);

        det *= mLu(i, i);
    }

    return det;
}

int main()
{
    ublas::matrix<int> m (3, 3);
    for (unsigned i = 0; i < m.size1 (); ++ i)
        for (unsigned j = 0; j < m.size2 (); ++ j)
            m (i, j) = 3 * i + j;
    std::cout << det_fast(m) << std::endl;
}

It is obvious that matrix m is singular, and if you change type from int to double functions returns zero. But with int it returns -48.

EDIT #1

Also, you can create ublas::matrix<int> without any errors and assign it to ublas::matrix<float>. So one way to correctly calculate determinant is to overload det_fast and remove define statement.

int det_fast(const ublas::matrix<int>& matrix)
{
    return (int)det_fast(ublas::matrix<float>(matrix));
}

EDIT #2

Take a look at the table, where average time is time of full determinant calculation (for int with copy operation) of 5 program runs.

size |  average time, ms
------------------------
     |    int      float
------------------------
100  |    9.0        9.0
200  |   46.6       46.8
300  |  156.4      159.4
400  |  337.4      331.4
500  |  590.8      609.2
600  |  928.2     1009.4
700  | 1493.8     1525.2
800  | 2162.8     2231.0
900  | 3184.2     3025.2

You might think that for int it is faster, but it's not true. In average for more runs I'm sure that you will get slight acceleration with float (I guess about 0-15 ms, which is time measure error). Also if we measure copy time we will see, that for matrix size less than 3000 x 3000 it is near zero (it is also about 0-15 ms, which is time measure error). For larger matrices copy time increases (for example for 5000 x 5000 matrix it is 125 ms). But there is one important note! Copy time is less than 1.0% of all determinant calculation for int type matrix and it is greatly decreasing with size grow!

All of this measures were made for a program compiled with Visual Studio 2013 in release configuration with boost version 1.58. Time were measured with clock function. CPU is Intel Xeon E5645 2.4GHz.

Also perfomance mostly depends on how your methods will be used. If you're going to call this function many times for small matrices than the copy time may become significant.

But there is a bug! Hmmm, so boost ublas didn`t aplicable for integer types? And there are no any way to fix it? It`i really sad, so I must implement own data types and operations for integer matrices.. — AeroSun, Apr 27 '15 at 11:52
Nikolay, at your example the matrix is really singular, but why it fails, for example on this matrix: {{4,5,6}{3,9,1}{7,8,2}} ? This matrix is non-singular but it fails with the same error, that "singular !=0". — AeroSun, Apr 27 '15 at 11:55
@AeroSun It fails not because it is singular, but because of type check. Look at [live](http://rextester.com/IHPBTB3552) version. Determinant is calculated correctly and it's ``-189`` for ``float`` and ``double``. For ``int`` the value is obviously wrong - ``-378``. If you remove ``define`` statement you will recive error since uBLAS is for floating types only. Also LU decomposition makes no sense for integral types, so it's not a bug. — Nikolay K, Apr 27 '15 at 12:12
Ok, this works (its too simple:)). But except that in this case created copy of int matrix....nothing free in this world...thx!) — AeroSun, Apr 27 '15 at 12:52
I made sequence of performance measurement with different matrix dimensions (form 10x10 to 1000x1000) and always you det_fast function for int matrix works faster then det_fast directly for float matrix. I didn`t get why....could you explain? — AeroSun, Apr 28 '15 at 13:20

boost::ublas how to get determinant of int matrix?

1 Answers1