对称矩阵求逆用CBLAS / LAPACKç [英] Symmetric Matrix Inversion in C using CBLAS/LAPACK
问题描述
我C语言编写的算法，需要矩阵与向量乘法。我有一个矩阵 问：的（W x宽），这是由转乘以一个向量的创建 Ĵ的（ 1 x宽），与自身和增加单位矩阵 I 的，使用标量的缩放的在的。
Q = [（J ^ T）* J + AI]
我接下来要乘Q 的逆是矢量G 以获得向量的 M
M =（Q ^（ - 1））。* G
我使用的 cblas 的和的 CLAPACK 的发展我的算法。当矩阵问：使用随机数（float类型）填充，并使用程序倒置的 sgetrf _ 的和的 sgetri _ 的，计算的逆为正确
但是，当矩阵Q是对称，这是当你乘（J ^ T）×j中的情况下，计算出的逆是错的！
我知道阵列的行主（C语言）和列为主（在FORTRAN）格式，同时调用的 LAPACK 的程序从C，但对于一个对称矩阵这不应该是一个问题为A ^ T = A。
我重视我的C函数code下面矩阵求逆。
我肯定有一个更好的办法来解决这个问题。谁能帮我这个？
使用cblas一个解决方案将是巨大的......
感谢。
无效InverseMatrix_R（浮动*矩阵，INT W）
{
诠释LDA = W;
INT IPIV [W]
INT ERR_INFO;
INT LWORK = W *宽;
浮动工作区[LWORK] // - 计算N矩阵A中的M的LU分解
sgetrf _（安培; W，放大器; W，矩阵，与放大器; LDA，IPIV，＆安培; ERR_INFO）; // - 生成矩阵的逆给予其LU decompsotion
sgetri _（安培; W，矩阵，与放大器; LDA，IPIV，工作区，和放大器; LWORK，＆安培; ERR_INFO）; // - 显示逆矩阵
PrintMatrix（矩阵，W，W）;}无效PrintMatrix（浮动*矩阵，诠释行，诠释科尔姆）
{
INT I，K; 对于（i = 0; I＆LT;排;我++）
{
对于（K = 0; K＆LT;科尔姆; k ++）
{
的printf（％G，矩阵[我*科尔姆+ K]）;
} 的printf（\\ n）;
}
}
我不知道BLAS或LAPACK，所以我不知道什么可能导致此行为。
但是，对于给定的形式的矩阵，计算逆是很容易的。造成这种情况的重要事实是
（J ^ T * j）^ 2 =（J ^ T * j）*（J ^ T * j）= J（1）T *（J * J（1）T）* J =＆LT;Ĵ| J＆GT; *（J ^ T * j）
其中，＆LT; U | V＆GT; 表示内积（如果组件是真正的 - 规范的双线性的形式复杂的组件，但那么你可能考虑不转，但共轭转，你会回来的内积）。
要概括，
（J ^ T * j）^ N =（小于Ĵ| J＆GT;）^（N-1）*（J ^ T * J）对于n＆GT; = 1。
让我们表示对称方阵（J ^ T * j）按取值和标量＆LT;Ĵ| J＆GT; 按①。然后，对于一般 A = 0 的足够大的绝对值（ | A |＆GT;问）：
（A * I + S）^（ - 1）= 1 / A *（I + A ^（ - 1）* S）^（ - 1）
= 1 / A *（I +Σ（-1）^ K *一个^（ - k）的* S ^ k）的
K＆大于0
= 1 / A *（I +（Σ（-1）^ K *一个^（ - 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 /一个* S - 1 /（A + q）* S - 1 /（A *（A + q））* S ^ 2
= I + 1 /一个* S - 1 /（A + Q）* S - Q /（A *（A + Q））* S
= I +（（A + Q） - 一 - Q）/（A *（A + Q））* S
= I
使用秒2 = Q * S 。
即计算也比先查找LU分解更简单和更有效
I am writing an algorithm in C that requires Matrix and Vector multiplications. I have a matrix Q (W x W) which is created by multiplying the transpose of a vector J(1 x W) with itself and adding Identity matrix I, scaled using scalar a.
Q = [(J^T) * J + aI].
I then have to multiply the inverse of Q with vector G to get vector M.
M = (Q^(-1)) * G.
I am using cblas and clapack to develop my algorithm. When matrix Q is populated using random numbers (type float) and inverted using the routines sgetrf_ and sgetri_ , the calculated inverse is correct.
But when matrix Q is symmetrical, which is the case when you multiply (J^T) x J, the calculated inverse is wrong!!.
I am aware of the row-major (in C) and column-major (in FORTRAN) format of arrays while calling lapack routines from C, but for a symmetrical matrix this should not be a problem as A^T = A.
I have attached my C function code for matrix inversion below.
I am sure there is a better way to solve this. Can anyone help me with this?
A solution using cblas would be great...
Thanks.
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.
这篇关于对称矩阵求逆用CBLAS / LAPACKç的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
