2

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

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

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

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

我正在使用cblas和clapack来开发我的算法。当矩阵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");
    }
}
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

That formula holds (by analyticity) for all a except a = 0 and a = -q, as can be verified by calculating

(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

using S^2 = q*S.

That calculation is also much simpler and more efficient than first finding the LU decomposition.

于 2012-05-23T14:27:34.750 回答
0

Example for 3x3 matrix inversion, visit sgetri.f for more

//__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

You may want to try Armadillo, which is an easy to use C++ wrapper for LAPACK. It provides several inverse related functions:

  • inv(), general inverse, with an optional speedup for symmetric positive definite matrices
  • pinv(), pseudo-inverse
  • solve(), solve a system of linear equations (that can be over- or under-determined), without doing the actual inverse
于 2012-08-01T08:25:13.510 回答