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

Posted

技术标签:

【中文标题】使用 CBLAS/LAPACK 在 C 中进行对称矩阵求逆【英文标题】:Symmetric Matrix Inversion in C using CBLAS/LAPACK 【发布时间】:2012-05-23 11:58:27 【问题描述】:

我正在用 C 语言编写一个需要矩阵和向量乘法的算法。我有一个矩阵 Q (W x W),它是通过将向量 J( 1 x W) 与自身并添加单位矩阵 I,使用标量 a 进行缩放。

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

然后我必须将 Q 的倒数向量 G 相乘以得到向量 M

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

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

【问题讨论】:

【参考方案1】:

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

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

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

其中&lt;u|v&gt; 表示内积(如果组件是实数 - 复杂组件的规范 双线性 形式,但是您可能不会考虑转置而是共轭转置,并且您'会回到内部产品)。

概括,

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

让我们用S 表示对称方阵(J^T*J),用q 表示标量&lt;J|J&gt;。那么,对于绝对值足够大的一般a != 0|a| &gt; 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 = 0a = -q 外,该公式对所有a 均成立(通过分析),可以通过计算验证

(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 分解更简单、更有效。

【讨论】:

非常感谢,您的回答非常有帮助。你们以前的 cmets,让我思考。当我第一次发布这个问题时,我只是在测试我的 C 函数。在测试阶段,我使用 (J^T) * J 创建了矩阵 QQ 的这个值被传递给我的反函数,正如您在早期的 cmets,它不能倒置。我应该完全计算 Q = (J^T) * J + aI 并将其传递给我的函数。我刚刚在matlab中尝试了这个理论,结果和你建议的一样!!! 我认为这应该可以解决我的代码问题。再次感谢您的帮助.. :)【参考方案2】:

您可能想尝试Armadillo,这是一个易于使用的 LAPACK 的 C++ 包装器。它提供了几个逆相关函数:

inv(),一般逆,具有对称正定矩阵的可选加速 pinv(),伪逆 solve(),求解一个线性方程组(可以超定或欠定),而不进行实际的逆运算

【讨论】:

【参考方案3】:

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

【讨论】:

以上是关于使用 CBLAS/LAPACK 在 C 中进行对称矩阵求逆的主要内容,如果未能解决你的问题,请参考以下文章

区块链科普篇——非对称加密

使用OpenSSL进行RSA非对称加密(C++版本)

使用OpenSSL进行RSA非对称加密(C++版本)

使用OpenSSL进行RSA非对称加密(C++版本)

使用OpenSSL进行RSA非对称加密(C++版本)

c ++ - 使用存储在堆上的数组填充对称矩阵