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

此函数工作正常，但仅适用于非整数类型（它适用于浮点和双精度）。当我尝试传递相同的矩阵但类型为 int 时，我收到编译错误：

在第 167 行的文件 c:\boost\boost_1_58_0\boost\numeric\ublas\lu.hpp 中检查失败：singular != 0 || detail::expression_type_check (prod (triangular_adaptor (m), triangular_adaptor (m)), cm) 未知位置 (0)：“BaseTest”中的致命错误：std::logic_error：内部逻辑

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

4

1 回答 1

2

当然，这不是错误，因为像这样的检查遍布 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 次程序运行average time的完整行列式计算（用于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 毫秒）。但是有一个重要的注意事项！复制时间不到类型矩阵所有行列式计算的1.0%int ，并且随着大小的增长而大大减少！

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

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

于 2015-04-26T16:35:49.443 回答