重温数学基础——矩阵求逆

目标：给定任意一个矩阵A，求它的逆矩阵。

方案一：直接求逆

【例1】：简单的二维矩阵

设$A=\left[\begin{array}{cc}a& b\\ c& d\end{array}\right]$，则$A^{-1}=\frac{1}{det(A)}\left[\begin{array}{cc}d& -b\\ -c& a\end{array}\right]$。
【编程实现】：
注意1.det(A)=ad-bc的值等于0时，矩阵A没有逆；
注意2.det(A)=ad-bc的值非常接近于0时，求得的逆$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=\left[\begin{array}{cccc}{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{array}\right]$，同时设它的逆矩阵$A^{-1}=\left[\begin{array}{cccc}{a}_{00}^{\prime }& a_{01}^{\prime }& a_{02}^{\prime }& a_{03}^{\prime }\\ a_{10}^{\prime }& a_{11}^{\prime }& a_{12}^{\prime }& a_{13}^{\prime }\\ a_{20}^{\prime }& a_{21}^{\prime }& a_{22}^{\prime }& a_{23}^{\prime }\\ a_{30}^{\prime }& a_{31}^{\prime }& a_{32}^{\prime }& a_{33}^{\prime }\end{array}\right]$，根据$A{A}^{-1}=I$，将$I$第一列对应的四个方程取出有：
$$\begin{array}{rl}{a}_{00}\ast a_{00}^{\prime }+a_{01}\ast a_{10}^{\prime }+a_{02}\ast a_{20}^{\prime }+a_{03}\ast a_{30}^{\prime }& =1\\ a_{11}\ast a_{10}^{\prime }+a_{12}\ast a_{20}^{\prime }+a_{13}\ast a_{30}^{\prime }& =0\\ a_{22}\ast a_{20}^{\prime }+a_{23}\ast a_{30}^{\prime }& =0\\ a_{33}\ast a_{30}^{\prime }& =0\end{array}$$
则，（四个方程四个未知量）我们可以求出$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分解