5

我正在尝试使用 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;
}
4

2 回答 2

3

置换矩阵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_factorizeand方法（咨询源代码，但基本上放弃对and的调用）。该更改将使您的行列式计算按原样工作。但是要小心：移除旋转会使算法在数值上不太稳定。lu_substitutepmlu_factorizelu_substitute

于 2009-09-14T06:15:38.940 回答
1

根据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”）。

于 2013-04-06T19:32:42.017 回答