Boost库，如何从lu_factorize()中获取行列式？

Question

我正在尝试使用boost c++库计算行列式。我找到了函数InvertMatrix()的代码，如下所示。每次计算逆矩阵时，我都想得到行列式。我知道如何计算，通过从LU分解的U矩阵对角线向下相乘。但是有一个问题，我能够正确地计算行列式，除了符号。根据枢轴选取的不同，有一半的时间我会得到错误的符号。有人有建议如何每次都得到正确的符号吗？提前感谢。

template<class T>
bool InvertMatrix(const ublas::matrix<T>& input, ublas::matrix<T>& inverse)
{
 using namespace boost::numeric::ublas;
 typedef permutation_matrix<std::size_t> pmatrix;
 // create a working copy of the input
 matrix<T> A(input);
 // create a permutation matrix for the LU-factorization
 pmatrix pm(A.size1());

 // perform LU-factorization
 int res = lu_factorize(A,pm);
 if( res != 0 ) return false;

这里是我尝试计算行列式的最佳方法。

 T determinant = 1;

 for(int i = 0; i < A.size1(); i++)
 {
  determinant *= A(i,i);
 }

结束我的代码部分。

 // create identity matrix of "inverse"
 inverse.assign(ublas::identity_matrix<T>(A.size1()));

 // backsubstitute to get the inverse
 lu_substitute(A, pm, inverse);

 return true;
}

- phoganuci

Answer 1

排列矩阵 pm 包含了您需要确定符号变化的信息：您需要将行列式乘以排列矩阵的行列式。

查看源文件lu.hpp，我们可以找到一个名为 swap_rows 的函数，它告诉我们如何将排列矩阵应用于矩阵。很容易修改该函数以产生排列矩阵的行列式（排列的符号），假设每个实际交换都会贡献一个因子-1。

template <typename size_type, typename A>
int determinant(const permutation_matrix<size_type,A>& pm)
{
    int pm_sign=1;
    size_type size=pm.size();
    for (size_type i = 0; i < size; ++i)
        if (i != pm(i))
            pm_sign* = -1; // swap_rows would swap a pair of rows here, so we change sign
    return pm_sign;
}

另一种选择是使用lu_factorize和lu_substitute方法，它们不执行任何枢轴操作（请参阅源代码，但基本上在调用lu_factorize和lu_substitute时删除pm）。 这个变化将使您的行列式计算按原样工作。 但是要小心：去除枢轴操作会使算法的数值稳定性降低。

- Managu

Answer 2

根据http://qiangsong.wordpress.com/2011/07/16/lu-factorisation-in-ublas/：

只需将determinant *= A(i,i)更改为determinant *= (pm(i) == i ? 1 : -1) * A(i,i)。 我尝试了这种方法，它有效。

我知道，实际上它与Managu的答案非常相似，思路也相同，但我相信它更简单（如果在InvertMatrix函数中使用，则是“2合1”）。

- jacek.ciach