2

行列とベクトルの乗算を必要とするアルゴリズムを C で作成しています。ベクトルJ (1 x W)の転置をそれ自体で乗算し、スカラーaを使用してスケーリングされた単位行列Iを追加することによって作成される行列Q (W x W) があります。

Q = [(J^T) * J + aI]。

次に、 Q の逆数にベクトルGを掛けて、ベクトルMを取得する必要があります。

M = (Q^(-1)) * G.

私はcblasとclapackを使用してアルゴリズムを開発しています。行列Qが乱数 (タイプ float) を使用して設定され、ルーチンsgetrf_およびsgetri_を使用して反転される場合、計算された逆行列は正しいです。

しかし、行列 Q が対称の場合、つまり (J^T) x J を掛けた場合、計算された逆数は間違っています!! .

C からlapackルーチンを呼び出す際に、配列の行優先 (C の場合) および列優先 (FORTRAN の場合) の形式を認識していますが、対称行列の場合、これは A^T = A であるため問題になりません。

行列反転の C 関数コードを以下に添付しました。

これを解決するためのより良い方法があると確信しています。誰でもこれで私を助けることができますか？

cblas を使用したソリューションは素晴らしいでしょう...

ありがとう。

void InverseMatrix_R(float *Matrix, int W)
{   
    int     LDA = W;
    int     IPIV[W];
    int     ERR_INFO;
    int     LWORK = W * W;
    float   Workspace[LWORK];

    // - Compute the LU factorization of a M by N matrix A
    sgetrf_(&W, &W, Matrix, &LDA, IPIV, &ERR_INFO);

    // - Generate inverse of the matrix given its LU decompsotion
    sgetri_(&W, Matrix, &LDA, IPIV, Workspace, &LWORK, &ERR_INFO);

    // - Display the Inverted matrix
    PrintMatrix(Matrix, W, W);

}

void PrintMatrix(float* Matrix, int row, int colm)
{
    int i,k;

    for (i =0; i < row; i++) 
    {
        for (k = 0; k < colm; k++) 
        {
            printf("%g, ",Matrix[i*colm + k]);
        }

        printf("\n");
    }
}
4

3 に答える 3

4

私は BLAS や LAPACK を知らないので、何がこの動作を引き起こすのかわかりません。

しかし、与えられた形式の行列の場合、逆行列を計算するのは非常に簡単です。そのために重要な事実は

(J^T*J)^2 = (J^T*J)*(J^T*J) = J^T*(J*J^T)*J = <J|J> * (J^T*J)

where<u|v>は内積を示します (コンポーネントが実数の場合 - 複雑なコンポーネントの正準双一次形式ですが、転置ではなく共役転置を考慮すると、内積に戻ることになります)。

一般化して、

(J^T*J)^n = (<J|J>)^(n-1) * (J^T*J), for n >= 1.

対称正方行列(J^T*J)をSで表し、スカラー<J|J>をで表しqます。次に、a != 0絶対値が十分に大きい一般 ( |a| > q) の場合:

(a*I + S)^(-1) = 1/a * (I + a^(-1)*S)^(-1)
               = 1/a * (I + ∑ (-1)^k * a^(-k) * S^k)
                           k>0
               = 1/a * (I + (∑ (-1)^k * a^(-k) * q^(k-1)) * S)
                            k>0
               = 1/a * (I - 1/(a+q)*S)
               = 1/a*I - 1/(a*(a+q))*S

その式は (分析により) と を除くすべてにa適用され、計算によって検証できます。a = 0a = -q

(a*I + S) * (1/a*I - 1/(a*(a+q))*S) = I + 1/a*S - 1/(a+q)*S - 1/(a*(a+q))*S^2
                                    = I + 1/a*S - 1/(a+q)*S - q/(a*(a+q))*S
                                    = I + ((a+q) - a - q)/(a*(a+q))*S
                                    = I

を使用してS^2 = q*S。

この計算は、最初に LU 分解を見つけるよりもはるかに単純で効率的です。

于 2012-05-23T14:27:34.750 に答える
0

3x3行列反転の例、詳細についてはsgetri.fにアクセスしてください

//__CLPK_integer is typedef of int
//__CLPK_real is typedef of float


__CLPK_integer ipiv[3];
{
    //Compute LU lower upper factorization of matrix
    __CLPK_integer m=3;
    __CLPK_integer n=3;
    __CLPK_real *a=(float *)this->m1;
    __CLPK_integer lda=3;
    __CLPK_integer info;
    sgetrf_(&m, &n, a, &lda, ipiv, &info);
}

{
    //compute inverse of a matrix
    __CLPK_integer n=3;
    __CLPK_real *a=(float *)this->m1;
    __CLPK_integer lda=3;
    __CLPK_real work[3];
    __CLPK_integer lwork=3;
    __CLPK_integer info;
    sgetri_(&n, a, &lda, ipiv, work, &lwork, &info);
}
于 2012-07-24T09:11:24.170 に答える
0

LAPACK の使いやすい C++ ラッパーであるArmadilloを試してみてください。いくつかの逆関連関数を提供します。

  • inv()、一般的な逆行列、対称正定行列のオプションの高速化
  • pinv()、疑似逆
  • solve()、実際の逆行列を実行せずに、一次方程式系 (過剰決定または過少決定の可能性がある) を解きます。
于 2012-08-01T08:25:13.510 に答える