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

Question

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

4

我找到了计算 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类型下正常工作）。但是，当我尝试传递相同的矩阵但类型为int时，我收到了编译错误：

“检查失败，位于c:\boost\boost_1_58_0\boost\numeric\ublas\lu.hpp的167行： singular != 0 || detail::expression_type_check (prod (triangular_adaptor (m), triangular_adaptor (m)), cm) 未知位置(0)："BaseTest"中的致命错误：std::logic_error: internal logic”

这是boost的bug还是我的函数有误？有什么改动可以避免这个错误？

- AeroSun

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Nikolay K · Accepted Answer

当然，这不是一个错误，因为在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类型，函数返回零。但是对于int类型，它返回-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

请看表格，其中平均时间是对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。

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