我使用C编写的算法，需要矩阵与向量乘法。 我有一个矩阵Q（W×宽），这是由矢量Ĵ（1×W）与自身的转置相乘并添加单位矩阵I创建的，缩放的用标量一个 。
Q = [（j ^ T）* J + Al]。
然后，我要乘Q的倒数与矢量G来获得向量M。
M =（Q ^（ - 1））* G.
我使用cblas和CLAPACK发展我的算法。 当矩阵Q被使用的随机数（float类型）填充和倒置使用例程sgetrf_和sgetri_，所计算出的逆是正确的 。
但是，当矩阵Q是对称的 ，这是当你乘法（j ^ T）×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");
}
}
I don't know BLAS or LAPACK, so I have no idea what may cause this behaviour.
But, for matrices of the given form, calculating the inverse is quite easy. The important fact for this is
(J^T*J)^2 = (J^T*J)*(J^T*J) = J^T*(J*J^T)*J = <J|J> * (J^T*J)
where <u|v> denotes the inner product (if the components are real - the canonical bilinear form for complex components, but then you'd probably consider not the transpose but the conjugate transpose, and you'd be back at the inner product).
Generalising,
(J^T*J)^n = (<J|J>)^(n-1) * (J^T*J), for n >= 1.
Let us denote the symmetric square matrix (J^T*J) by S and the scalar <J|J> by q. Then, for general a != 0 of sufficiently large absolute value (|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.
你可能想尝试犰狳 ，这是一个易于使用的C ++包装的LAPACK。 它提供了一些逆相关的功能：
- INV（） ，一般逆，具有对称正定矩阵的可选加速
- PINV（） ，伪逆
- 解决（） ，求解线性方程系统（可以是过表达或低测定），而不进行实际的逆
例如，对于3×3的矩阵求逆，访问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);
}
