为什么要少用逆矩阵
我经常会在线性代数教材以及论坛讨论中看到不建议使用逆矩阵$\mathbf{A}^{-1}$来求解线性方程$\mathbf{A}\mathbf{x}=\mathbf{b}$,尽管我一直遵循这样的原则(实践中逆矩阵确实不够稳定),但仍然不明白不使用逆矩阵的理由。本文总结了我在网上看到的一些关于逆矩阵的讨论,希望能解释为什么要少用逆矩阵来求解线性方程。
数学原理
解线性方程组(逆矩阵)
考虑求解线性方程组中的$\mathbf{x} \in \mathbb{R}^n$:
\[\begin{equation} \mathbf{A} \mathbf{x} = \mathbf{b} \end{equation}\]其中$\mathbf{A} \in \mathbb{R}^{n \times n},\mathbf{b} \in \mathbb{R}^n$。
显然一种简单直观的求解方法是计算$\mathbf{A}$的逆矩阵$\mathbf{A}^{-1}$:
\[\begin{equation} \mathbf{x} = \mathbf{A}^{-1} \mathbf{b} \end{equation}\]求解逆矩阵通常需要$2n^3$次浮点数操作(floating-point operations,flops),此外计算$\mathbf{A}^{-1}\mathbf{b}$需要$2n^2$ flops,总共$2n^3+2n^2$ flops。
不过实践中不推荐逆矩阵求解,因为可能在准确性上产生问题。求解线性方程更常用的是分解方法,比如接下来的LU分解。
解线性方程组(LU分解)
LU分解(lower–upper decomposition)是一种常用的矩阵分解方法,基本思路是利用高斯消元法将$\mathbf{A}$转化为上三角矩阵$\mathbf{U}$,主要过程如下所示:
\[\begin{equation}\nonumber \overbrace{\begin{bmatrix} \times & \times & \times\\ \times & \times & \times\\ \times & \times & \times\\ \end{bmatrix}}^{\mathbf{A}} \rightarrow \overbrace{\begin{bmatrix} \times & \times & \times\\ 0 & \times & \times\\ 0 & \times & \times\\ \end{bmatrix}}^{\mathbf{L}_1\mathbf{A}} \rightarrow \overbrace{\begin{bmatrix} \times & \times & \times\\ 0 & \times & \times\\ 0 & 0 & \times\\ \end{bmatrix}}^{\mathbf{L}_2\mathbf{L}_1\mathbf{A}} \end{equation}\]其中$\mathbf{L}_1,\mathbf{L}_2$是一系列下三角矩阵,用来表示第n行加上第n-1行乘以某个数的消元操作,所以矩阵$\mathbf{A}$可以消元成为一个上三角矩阵$\mathbf{U}$:
\[\begin{equation} \mathbf{L}_{n-1}\mathbf{L}_{n-2} \cdots \mathbf{L}_{1}\mathbf{A} = \mathbf{U} \end{equation}\]矩阵$\mathbf{A}$的LU分解为:
\[\begin{equation} \begin{split} \mathbf{A} &= \left(\mathbf{L}_{1}^{-1} \cdots \mathbf{L}_{n-2}^{-1} \mathbf{L}_{n-1}^{-1}\right) \mathbf{U}\\ &= \mathbf{L} \mathbf{U} \end{split} \end{equation}\]下三角矩阵的逆依然是下三角矩阵,并且一系列下三角矩阵的乘积也是一个三角矩阵,因此$\mathbf{L}$是下三角矩阵。此外,虽然涉及到求$\mathbf{L}_i$的逆,但是这一步骤实施起来是非常简单的,$\mathbf{L}_i$主对角线元素不变、下三角元素全部乘以-1即可得到$\mathbf{L}_i^{-1}$。因为$\mathbf{L}_i$每次都只对第$i$列上的主元做消元操作,所以下三角部分只会在第$i$列存在不为0的元素。以$3 \times 3$矩阵为例,$\mathbf{L}_1$可以写成如下的形式:
\[\begin{equation}\nonumber \begin{bmatrix} 1 & 0 & 0\\ a & 1 & 0\\ b & 0 & 1\\ \end{bmatrix} \end{equation}\]很自然的可以证明$\mathbf{L}_1^{-1}$为:
\[\begin{equation}\nonumber \begin{bmatrix} 1 & 0 & 0\\ -a & 1 & 0\\ -b & 0 & 1\\ \end{bmatrix} \end{equation}\]上述步骤的LU分解需要的浮点数操作近似为$\frac{2}{3}n^3$ flops。
接下来我们看看如何利用LU分解解线性方程组。约定$\mathbf{A} \backslash \mathbf{b}$代表$\mathbf{A}\mathbf{x}=\mathbf{b}$的解$\mathbf{x}$(这也是MATLAB求解线性方程组的代码),按照以下步骤求解出$\mathbf{x}$:
\[\begin{equation} \begin{split} \mathbf{L} \mathbf{U} &= \mathbf{A}\\ \mathbf{L} \backslash \mathbf{b} &= \mathbf{y}\\ \mathbf{U} \backslash \mathbf{y} &= \mathbf{x}\\ \end{split} \end{equation}\]所以先求解$\mathbf{L}\mathbf{y}=\mathbf{b}$得到中间变量$\mathbf{y}$,再求解$\mathbf{U}\mathbf{x}=\mathbf{y}$得到$\mathbf{x}$,看上去比直接求解$\mathbf{A}\mathbf{x} = \mathbf{b}$复杂了不少,为什么要这么做呢?
原因在于$\mathbf{L}$和$\mathbf{U}$作为三角矩阵能显著简化线性方程组求解,比如$\mathbf{L}\mathbf{y}=\mathbf{b}$,我们有:
\[\begin{equation} \begin{bmatrix} 1 & & & & \\ \vdots&\ddots& & &\\ l_{i,1}&\cdots&1& &\\ \vdots&\vdots&\vdots&\ddots&\\ l_{n,1}&\cdots&\cdots&\cdots&1\\ \end{bmatrix} \begin{bmatrix} y_1\\ \vdots\\ y_i\\ \vdots\\ y_n\\ \end{bmatrix} = \begin{bmatrix} b_1\\ \vdots\\ b_i\\ \vdots\\ b_n\\ \end{bmatrix} \end{equation}\]可以用递归方法求解$\mathbf{y}$:
\[\begin{equation} \begin{split} y_1 &= b_1\\ y_i &= b_i - l_{i,1}y_1 - \cdots - l_{i, i-1}y_{i-1} \end{split} \end{equation}\]这种思路叫forward substitution,即从上往下求解$\mathbf{y}$。之后可以用相似的方式从下往上求解$\mathbf{x}$,叫backward substitution。forward substitution的成本为(每一个元素各需要$i-1$次乘法和减法)$n^2-n$ flops,同理可得backward substitution的计算成本为$n^2+n$ flops,所以用LU分解求解线性方程组的总计算成本为$\frac{2}{3}n^3+2n^2$ flops。
解线性方程组(QR分解)
QR分解(QR decomposition)也常用于求解线性方程组。QR分解将$\mathbf{A} \in \mathbb{R}^{n \times n}$分解为正交矩阵$\mathbf{Q} \in \mathbb{R}^{n \times n}$和上三角矩阵$\mathbf{R} \in \mathbb{R}^{n \times n}$:
\[\begin{equation} \mathbf{A} = \mathbf{Q} \mathbf{R} \end{equation}\]对于正交矩阵$\mathbf{Q}$,它的逆矩阵$\mathbf{Q}^{-1}$实际就是$\mathbf{Q}^T$本身,所以对于线性方程组可以简化为:
\[\begin{equation} \begin{split} \mathbf{A} \mathbf{x} &= \mathbf{b}\\ \mathbf{Q}\mathbf{R}\mathbf{x} &= \mathbf{b}\\ \mathbf{R}\mathbf{x} &= \mathbf{Q}^T\mathbf{b}\\ \end{split} \end{equation}\]因为$\mathbf{R}$是上三角矩阵,对于$\mathbf{x}$的求解可以采用LU分解中介绍的backward substitution方式求解,从而避免了求逆的操作。
QR分解可以采用Gram-Schmidt正交化过程计算得到,这一过程与正交的几何表示密切相关,因此易于理解,不过这种正交化过程容易数值不稳定,因此很少直接在实践中使用。实践中更常用的是使用Givens rotations方法来实现QR分解。
Givens rotations的步骤与LU分解相似,对于矩阵$\mathbf{A}$,采用一系列旋转矩阵$\mathbf{G}_i^j$将$\mathbf{A}$逐步转化为上三角矩阵:
\[\begin{equation}\nonumber \overbrace{\begin{bmatrix} \times & \times & \times\\ \times & \times & \times\\ \times & \times & \times\\ \end{bmatrix}}^{\mathbf{A}} \rightarrow \overbrace{\begin{bmatrix} \times & \times & \times\\ \times & \times & \times\\ 0 & \times & \times\\ \end{bmatrix}}^{\mathbf{G}_3^1\mathbf{A}} \rightarrow \overbrace{\begin{bmatrix} \times & \times & \times\\ 0 & \times & \times\\ 0 & \times & \times\\ \end{bmatrix}}^{\mathbf{G}_2^1\mathbf{G}_3^1\mathbf{A}} \rightarrow \overbrace{\begin{bmatrix} \times & \times & \times\\ 0 & \times & \times\\ 0 & 0 & \times\\ \end{bmatrix}}^{\mathbf{G}_3^2\mathbf{G}_2^1\mathbf{G}_3^1\mathbf{A}} \end{equation}\]对于需要置0的第$a_{ij}$个元素$(i \gt j)$,$\mathbf{G}_i^j$构造为:
\[\begin{equation} \mathbf{G}_i^j = \begin{bmatrix} 1 & \cdots & 0 & \cdots & 0 & \cdots & 0\\ \vdots & \ddots & \vdots & & \vdots & & \vdots\\ 0 & \cdots & c_{jj} & \cdots & -s_{ji} & \cdots & 0\\ \vdots & & \vdots & \ddots & \vdots & & \vdots\\ 0 & \cdots & s_{ij} & \cdots & c_{ii} & \cdots & 0\\ \vdots & & \vdots & & \vdots & \ddots & \vdots\\ 0 & \cdots & 0 & \cdots & 0 & \cdots & 1\\ \end{bmatrix} \end{equation}\]其中$c_{jj} = c_{ii} = \cos(\theta)$,$s_{ij} = s_{ji} = \sin(\theta)$。显而易见这是一个对$ij$平面上的向量进行旋转的矩阵,对于任意向量,若想将其在第$i$个轴上的分量置0,则需要将其按$\theta$角旋转至对应的坐标轴上,旋转变化不会改变除第$i$、$j$行外的其他元素,因此通过反复运用旋转变换,可以将$\mathbf{A}$的下三角区域置0,有:
\[\begin{equation} \mathbf{G}_n^{n-1} \cdots \mathbf{G}_{n-1}^1 \mathbf{G}_n^1 \mathbf{A} = \mathbf{R} \end{equation}\]又因为旋转矩阵本身就是正交矩阵,所以其逆就是矩阵转置本身,所以有:
\[\begin{equation} \begin{split} \mathbf{A} &= \left(\mathbf{G}_n^{n-1} \cdots \mathbf{G}_{n-1}^1 \mathbf{G}_n^1\right)^T \mathbf{R}\\ &= \mathbf{Q}\mathbf{R} \end{split} \end{equation}\]Givens rotations实现的QR分解的复杂度大约是$\frac{7}{3}n^3$ flops(我也不确定),优势在于可以较为容易的实现并行化操作。
复杂度分析
从上面的分析中可以看出,LU分解大概需要$\frac{2}{3}n^3$ flops,逆矩阵大概需要$2n^3$ flops,是LU分解的3倍,因此线性方程组求解采用逆矩阵方法会更慢一些。尽管逆矩阵求解存在更快的算法,但是通常来说LU分解还是更有效率,所以一般而言还是应当直接求解线性方程组而不是计算逆矩阵。
此外,逆矩阵计算往往需要更大的内存。如果$\mathbf{A}$是稀疏矩阵(大多数元素是0),则可以使用更加有效率的稀疏结构存储。但是$\mathbf{A}^{-1}$通常是稠密的,因此存储逆矩阵需要更多的内存。
准确性分析
逆矩阵的另一个缺点在于求解的误差更大,特别在$\mathbf{A}$是病态矩阵的情况下尤为明显。
数值分析领域常用前向误差(forward error)和反向误差(backward error)来分析误差。比如对于函数$y=f(x)$,用$\hat{y}$来近似输出$y$,那么前向误差被定义为$|\hat{y}-y|$(这似乎就是我们常说的误差),而反向误差则衡量问题的敏感性,即为了产生近似解,输入数据同真实数据的偏移情况。$\hat{y}$的反向误差是满足$\hat{y}=f(x+\Delta x)$的最小$\Delta x$。
前向误差和反向误差联合起来可以分析一个问题是良态(well-conditioned)还是病态(ill-conditioned)的。对于良态问题,输入较小的改变会导致输出产生较小的改变;而病态问题,输入较小的改变会导致输出较大的改变。条件数(condition number)就是用来判断良态和病态的工具,条件数越大,输入很小的差异就会产生较大的输出误差,更倾向于病态问题。矩阵的条件数$\kappa(\mathbf{A})$定义为:
\[\begin{equation} \kappa(\mathbf{A}) = \|\mathbf{A}\| \|\mathbf{A}^{-1}\| \end{equation}\]如果采用矩阵的2范数,上式可以进一步简化为最大特征值$\lambda_1$与最小特征值$\lambda_n$之比:
\[\begin{equation} \kappa(\mathbf{A}) = \frac{\lambda_1}{\lambda_n} \end{equation}\]stackexchange和这篇blog详细分析了前向误差$|\mathbf{x}-\hat{\mathbf{x}}|$和反向误差$|\mathbf{b}-\mathbf{A}\hat{\mathbf{x}}|$(注意这里$\mathbf{x}$是我们求解问题的输出),有如下结论:
- 对于良态矩阵$\mathbf{A}$,LU分解的前向误差和逆矩阵的前向误差是很接近的
- 对于病态矩阵$\mathbf{A}$,逆矩阵的反向误差会远大于LU分解的反向误差
- 即使对于良态矩阵$\mathbf{A}$,逆矩阵的反向误差也比LU分解的反向误差更大
此外,前向误差和反向误差有如下不等式关系:
\[\begin{equation} \|\mathbf{x}-\hat{\mathbf{x}}\| \le \kappa(\mathbf{A}) \|\mathbf{b}-\mathbf{A}\hat{\mathbf{x}}\| \end{equation}\]综合来看的话,利用LU分解、QR分解等方法直接求解线性方程组比逆矩阵求解要更加准确。