5

BoostC++ライブラリを使用して行列式を計算しようとしています。以下にコピーした関数InvertMatrix（）のコードを見つけました。この逆数を計算するたびに、行列式も必要です。LU分解からU行列の対角線を掛けて計算する方法がわかります。1つの問題があります。符号を除いて、行列式を正しく計算できます。ピボットに応じて、半分の時間で間違ったサインが表示されます。誰かが毎回正しいサインを取得する方法についての提案がありますか？前もって感謝します。

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

もう1つの方法は、ピボットを実行しないlu_factorizeandメソッドを使用することです（ソースを参照しますが、基本的にとへの呼び出しをドロップします）。この変更により、行列式の計算がそのまま機能するようになります。ただし、注意してください。ピボットを削除すると、アルゴリズムの数値的安定性が低下します。lu_substitutepmlu_factorizelu_substitute

于 2009-09-14T06:15:38.940 に答える