boost::ublas如何获取int矩阵的行列式？

Question

我找到了计算boost::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;
}

此函数工作正常，但仅适用于非整数类型（它适用于 float 和 double）。当我尝试传递 same 矩阵但类型为 int 时，我收到编译错误：

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

是boost的bug还是我的功能有问题？我可以做些什么来避免这个错误？

Answer 1

当然这不是错误，因为像这样的检查遍及 uBLAS。我想这是因为它的大部分操作对非浮点类型没有意义。

您可以使用

禁用类型检查

#define BOOST_UBLAS_TYPE_CHECK 0

在包括之前。但请三思！看例子

#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;
}

很明显，矩阵 m 是奇异的，如果将类型从 int 更改为 double 函数 returns 为零。但是 int 它 returns -48.

编辑#1

此外，您可以毫无错误地创建 ublas::matrix<int> 并将其分配给 ublas::matrix<float>。因此，正确计算行列式的一种方法是重载 det_fast 并删除 define 语句。

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

编辑#2

看看 table，其中 average time 是 5 个程序运行的完整行列式计算时间（对于 int 具有复制操作）。

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

您可能认为 int 它更快，但事实并非如此。平均而言，对于更多的运行，我相信你会得到 float 的轻微加速（我猜大约 0-15 毫秒，这是时间测量误差）。此外，如果我们测量复制时间，我们将看到，对于小于 3000 x 3000 的矩阵大小，它接近于零（它也大约为 0-15 毫秒，这是时间测量误差）。对于较大的矩阵，复制时间增加（例如，对于 5000 x 5000 矩阵，它是 125 毫秒）。但是有一个重要的注意事项！ int 类型矩阵的所有行列式计算的复制时间小于 1.0%，并且随着大小的增加大大减少！

所有这些措施都是针对使用 Visual Studio 2013 在发布配置中使用 boost 版本 1.58 编译的程序进行的。时间是用 clock 函数测量的。 CPU 是 Intel Xeon E5645 2.4GHz。

此外，性能主要取决于您的方法将如何使用。如果您要为小矩阵多次调用此函数，那么复制时间可能会变得很重要。

boost::ublas如何获取int矩阵的行列式？

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

c++

math

boost-ublas