c - 使用 CBLAS/LAPACK 在 C 中进行对称矩阵求逆

Question

我正在用 C 语言编写一个需要矩阵和 vector 乘法的算法。我有一个矩阵 问 (W x W) 通过将 vector 的转置相乘 创建J (1 x W) 与自身并添加单位矩阵 我 , 使用标量 缩放一个 .

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

然后我必须乘以 Q的倒数与 vector G 获取 vector 手机 .

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

我正在使用 cblas 和 clapack 来开发我的算法。当矩阵问 使用随机数(浮点类型)填充并使用例程 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");
    }
}

Answer 1

我不知道 BLAS 或 LAPACK，所以我不知道是什么导致了这种行为。

但是，对于给定形式的矩阵，计算逆很容易。重要的事实是

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

在哪里 <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 = 0和 a = -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 分解更简单、更有效。