目标:给定任意一个矩阵A,求它的逆矩阵。
方案一:直接求逆
【例1】:简单的二维矩阵
设
A
=
[
a
b
c
d
]
A=\begin{bmatrix} a & b \\ c & d\end{bmatrix}
A=[acbd],则
A
−
1
=
1
d
e
t
(
A
)
[
d
−
b
−
c
a
]
A^{-1}=\frac{1}{det(A)} \begin{bmatrix} d & -b \\ -c & a\end{bmatrix}
A−1=det(A)1[d−c−ba]。
【编程实现】:
注意1.det(A)=ad-bc
的值等于0时,矩阵A没有逆;
注意2.det(A)=ad-bc
的值非常接近于0时,求得的逆
A
−
1
A^{-1}
A−1可能是不准确的(由于一些数值舍入误差?)。
int QRInvM2(double A[2][2]) {
double detA = A[0][0] * A[1][1] - A[0][1] * A[1][0];
if ((detA < 1.0e-20) && (detA > -1.0e-20))
return(0);
double t = 1/detA ;
A[0][0] = A[0][0] * t;
A[1][1] = A[1][1] * t;
A[0][1] = -A[0][1] * t;
A[1][0] = -A[1][0] * t;
t = A[0][0];
A[0][0] = A[1][1];
A[1][1] = t;
printf("%lf\n", detA);
printf("%lf, %lf\n", A[0][0], A[0][1]);
printf("%lf, %lf\n", A[1][0], A[1][1]);
return(1);
}
【例2】:特殊的上三角矩阵
设矩阵
A
=
[
a
00
a
01
a
02
a
03
0
a
11
a
12
a
13
0
0
a
22
a
23
0
0
0
a
33
]
A=\begin{bmatrix} a_{00} & a_{01} & a_{02} & a_{03} \\ 0 & a_{11} & a_{12} & a_{13} \\ 0 & 0 & a_{22} & a_{23} \\ 0 & 0 & 0 & a_{33} \\ \end{bmatrix}
A=⎣⎢⎢⎡a00000a01a1100a02a12a220a03a13a23a33⎦⎥⎥⎤,同时设它的逆矩阵
A
−
1
=
[
a
00
′
a
01
′
a
02
′
a
03
′
a
10
′
a
11
′
a
12
′
a
13
′
a
20
′
a
21
′
a
22
′
a
23
′
a
30
′
a
31
′
a
32
′
a
33
′
]
A^{-1}=\begin{bmatrix} a'_{00} & a'_{01} & a'_{02} & a'_{03} \\ a'_{10} & a'_{11} & a'_{12} & a'_{13} \\ a'_{20} & a'_{21} & a'_{22} & a'_{23} \\ a'_{30} & a'_{31} & a'_{32} & a'_{33} \\ \end{bmatrix}
A−1=⎣⎢⎢⎡a00′a10′a20′a30′a01′a11′a21′a31′a02′a12′a22′a32′a03′a13′a23′a33′⎦⎥⎥⎤,根据
A
A
−
1
=
I
AA^{-1}=I
AA−1=I,将
I
I
I第一列对应的四个方程取出有:
a
00
∗
a
00
′
+
a
01
∗
a
10
′
+
a
02
∗
a
20
′
+
a
03
∗
a
30
′
=
1
a
11
∗
a
10
′
+
a
12
∗
a
20
′
+
a
13
∗
a
30
′
=
0
a
22
∗
a
20
′
+
a
23
∗
a
30
′
=
0
a
33
∗
a
30
′
=
0
\begin{aligned} a_{00}*a'_{00}+a_{01}*a'_{10}+a_{02}*a'_{20}+a_{03}*a'_{30} & =1 \\ a_{11}*a'_{10}+a_{12}*a'_{20}+a_{13}*a'_{30} &=0 \\ a_{22}*a'_{20}+a_{23}*a'_{30} &=0 \\ a_{33}*a'_{30} &=0 \end{aligned}
a00∗a00′+a01∗a10′+a02∗a20′+a03∗a30′a11∗a10′+a12∗a20′+a13∗a30′a22∗a20′+a23∗a30′a33∗a30′=1=0=0=0
则,(四个方程四个未知量)我们可以求出
A
−
1
A^{-1}
A−1的第一列。后面三列求取的思路与此类似。
【编程实现】:
int UInvM4(double U[4][4]) {
double detU = U[0][0] * U[1][1] * U[2][2] * U[3][3];
if ((detU < 1.0e-20) && (detU > -1.0e-20))
return(0);
double invU[4][4];
for (int i = 0; i < 4; i++) {
for (int j = 0; j < 4; j++) {
invU[i][j] = 0.0;
}
}
double b[4];
for (int k = 0; k < 4; k++) {
for (int i = 0; i < 4; i++) {
b[i] = 0.0;
}
b[k] = 1.0;
invU[3][k] = b[3] / U[3][3];
invU[2][k] = (b[2] - U[2][3] * invU[3][k]) / U[2][2];
invU[1][k] = (b[1] - U[1][3] * invU[3][k] - U[1][2] * invU[2][k]) / U[1][1];
invU[0][k] = (b[0] - U[0][3] * invU[3][k] - U[0][2] * invU[2][k] - U[0][1] * invU[1][k]) / U[0][0];
}
return(1);
}
方案二:LU分解