病态矩阵与条件数

本文的阅读等级：高级

当一线性系统受到极微小的扰动即可引发方程解剧烈变化时，我们将无从信任计算结果，便称它是病态系统(见“ 病态系统 ”)。条件数(condition number) 是矩阵运算误差分析的基本工具，它可以度量矩阵对于数值计算的敏感性和稳定性，也可以用来检定病态系统。本文通过一个简单的线性方程扰动问题介绍条件数的推导过程，基本推演工具是矩阵范数 $Vert AVert$ 的定义所含的两个不等式(见“ 矩阵范数 ”)：

$Vert A+BVertleVert AVert+Vert BVert$

$Vert ABVertleVert AVertcdotVert BVert$

问题一 ：考虑线性方程 $Amathbf{x}=mathbf{b}$ ，其中 $A$ 为 $n imes n$ 阶可逆矩阵。假设 $mathbf{b}$ 受到微小扰动成为 $mathbf{b}+deltamathbf{b}$ ，方程解随之由 $mathbf{x}$ 变为 $mathbf{x}+deltamathbf{x}$ ，试计算相对误差 $frac{Vertdeltamathbf{x}Vert}{Vertmathbf{x}Vert}$ 的可能范围。

考虑包含扰动及误差的线性方程：

$A(mathbf{x}+deltamathbf{x})=mathbf{b}+deltamathbf{b}$

与原方程式 $Amathbf{x}=mathbf{b}$ 相减可得 $A(deltamathbf{x})=deltamathbf{b}$ ，方程解的误差为 $deltamathbf{x}=A^{-1}(deltamathbf{b})$ 。利用不等式 $Vertdeltamathbf{x}Vert=Vert A^{-1}(deltamathbf{b})VertleVert A^{-1}VertcdotVertdeltamathbf{b}Vert$ 和 $Vertmathbf{b}Vert=Vert Amathbf{x}VertleVert AVertcdotVertmathbf{x}Vert$ ，可得

$displaystylefrac{Vertdeltamathbf{x}Vert}{Vertmathbf{x}Vert}=frac{Vert A^{-1}(deltamathbf{b})Vert}{Vertmathbf{x}Vert}lefrac{Vert A^{-1}VertcdotVertdeltamathbf{b}Vert}{Vertmathbf{x}Vert}leleft(Vert AVertcdotVert A^{-1}Vert ight)frac{ Vertdeltamathbf{b}Vert}{Vertmathbf{b}Vert}$

另一方面，我们考虑不等式 $Vertdeltamathbf{b}Vert=Vert A(deltamathbf{x})VertleVert AVertcdotVertdeltamathbf{x}Vert$ 和 $Vertmathbf{x}Vert=Vert A^{-1}mathbf{b}VertleVert A^{-1}VertcdotVertmathbf{b}Vert$ ，就有

$displaystyle frac{Vertdeltamathbf{x}Vert}{Vertmathbf{x}Vert}gefrac{Vertdeltamathbf{b}Vert}{Vert AVertcdotVertmathbf{x}Vert}geleft(Vert AVertcdotVert A^{-1}Vert ight)^{-1}frac{Vertdeltamathbf{b}Vert}{Vertmathbf{b}Vert}$

上面二式指出相对误差 $frac{Vertdeltamathbf{x}Vert}{Vertmathbf{x}Vert}$ 的上下界由矩阵范数乘积 $Vert AVert cdotVert A^{-1}Vert$ 决定，根据这个结果我们定义方阵 $A$ 的条件数为

$kappa(A)=Vert AVertcdotVert A^{-1}Vert$

也就是说， $deltamathbf{b}$ 引起的相对误差大小由内部因素 $kappa(A)$ 和外部因素 $frac{Vertdeltamathbf{b}Vert}{Vertmathbf{b}Vert}$ 共同决定：

$displaystyle kappa(A)^{-1}frac{Vertdeltamathbf{b}Vert}{Vertmathbf{b}Vert}lefrac{Vertdeltamathbf{x}Vert}{Vertmathbf{x}Vert}lekappa(A)frac{Vertdeltamathbf{b}Vert}{Vertmathbf{b}Vert}$

接着讨论条件数的计算方法。设 $n imes n$ 阶矩阵 $A$ 的奇异值分解为 $A=USigma V^T$ ，其中 $Sigma=mathrm{diag}(sigma_1,ldots,sigma_n)$ ， $sigma_1gecdotsgesigma_nge 0$ ， $U$ 和 $V$ 为正交矩阵(见“ 奇异值分解(SVD) ”)。若 $A$ 为可逆矩阵， $mathrm{rank}A=n$ ，所以 $sigma_i>0$ ， $i=1,ldots,n$ 。因为正交矩阵 $U$ 和 $V$ 不改变向量长度，最大奇异值和最小奇异值分别满足

$displaystyle sigma_1=max_{Vertmathbf{x}Vert=1}VertSigmamathbf{x}Vert=max_{Vertmathbf{x}Vert=1}Vert Amathbf{x}Vert=Vert AVert_2$

$displaystyle sigma_n=min_{Vertmathbf{x}Vert=1}VertSigmamathbf{x}Vert=min_{Vertmathbf{x}Vert=1}Vert Amathbf{x}Vert=frac{1}{Vert A^{-1}Vert_2}$

以2-范数(2-norm) 计算的条件数也就为

$displaystyle kappa(A)=Vert AVert_2cdotVert A^{-1}Vert_2=frac{sigma_1}{sigma_n}ge 1$

若 $A$ 为对称可逆矩阵，计算

$A^2=AA^T=(USigma V^T)(VSigma U^T)=USigma^2U^T$

因此 $sigma_i=vertlambda_ivert$ ， $i=1,ldots,n$ ， $lambda_i$ 是 $A$ 的特征值，条件数可表示为

$displaystyle kappa(A)=left|frac{lambda_1}{lambda_n} ight|$

又若 $A$ 为正交矩阵， $A^T=A^{-1}$ ，则 $sigma_1=cdots=sigma_n=1$ ，故 $kappa(A)=1$ ，这说明正交矩阵具有最佳的数值稳定性。

下面我们进一步讨论奇异值分解与条件数上下界的关系。奇异值分解的 $U$ 矩阵其行向量 $mathbf{u}_1,ldots,mathbf{u}_n$ 称为左奇异向量， $V$ 矩阵的行向量 $mathbf{v}_1,ldots,mathbf{v}_n$ 称为右奇异向量，左右奇异向量的关系如下：

$Amathbf{v}_j=sigma_jmathbf{u}_j,~~~j=1,ldots,n$

考虑 $mathbf{b}=alphamathbf{u}_1$ ， $deltamathbf{b}=etamathbf{u}_n$ ，方程解和误差可分别表示为

$displaystyle mathbf{x}=A^{-1}mathbf{b}=A^{-1}(alphamathbf{u}_1)=frac{alpha}{sigma_1}mathbf{v}_1$

$displaystyle deltamathbf{x}=A^{-1}(deltamathbf{b})=A^{-1}(etamathbf{u}_n)=frac{eta}{sigma_n}mathbf{v}_n$

计算向量长度并相除，

$displaystylefrac{Vertdeltamathbf{x}Vert}{Vertmathbf{x}Vert}=left(frac{sigma_1}{sigma_n} ight)frac{vertetavert}{vertalphavert}=kappa(A)frac{Vertdeltamathbf{b}Vert}{Vertmathbf{b}Vert}$

得知 $mathbf{b}=alphamathbf{u}_1$ ， $deltamathbf{b}=etamathbf{u}_n$ 满足相对误差的上界。同样方式也可以证明 $mathbf{b}=alphamathbf{u}_n$ ， $deltamathbf{b}=etamathbf{u}_1$ 满足相对误差的下界，亦即

$displaystylefrac{Vertdeltamathbf{x}Vert}{Vertmathbf{x}Vert}=left(frac{sigma_n}{sigma_1} ight)frac{vertetavert}{vertalphavert}=kappa(A)^{-1}frac{Vertdeltamathbf{b}Vert}{Vertmathbf{b}Vert}$

考虑下面两个取自“ 病态系统 ”一文的例子：

$A=egin{bmatrix} .540&.387\ .647&.323 end{bmatrix},~B=egin{bmatrix} .540&.323\ .647&.387 end{bmatrix}$

计算奇异值后代入条件数公式：

$egin{aligned} kappa(A)&=.978920/.077605=12.614\ kappa(B)&=.981991/.000001=981,991approx 10^{6}end{aligned}$

此例中， $kappa(A)$ 小(指 $10^{gamma}$ 的幂 $gamma$ )， $A$ 是良置的，小扰动 $deltamathbf{b}$ 不至于对方程解造成大误差。但 $kappa(B)$ 很大， $B$ 是病态的，小扰动 $deltamathbf{b}$ 可能引起大 $deltamathbf{x}$ ，但大扰动 $deltamathbf{b}$ 也可能产生微小的 $deltamathbf{x}$ 。由于难以掌握 $deltamathbf{b}$ 的变化方向，在面对病态系统时我们总是要预防最坏的情形发生。

问题二 ：如果线性方程 $Amathbf{x}=mathbf{b}$ 的系数矩阵 $A$ 和常数向量 $mathbf{b}$ 都受到杂音干扰，这对方程解 $mathbf{x}$ 的误差有何影响？

考虑 $A$ 和 $mathbf{b}$ 分别受到 $delta A$ 和 $deltamathbf{b}$ 干扰，则 $deltamathbf{x}$ 满足

$(A+delta A)(mathbf{x}+deltamathbf{x})=mathbf{b}+deltamathbf{b}$

为简化问题，以下假设 $frac{Vertdelta AVert}{Vert AVert}leepsilon$ ， $frac{Vertdeltamathbf{b}Vert}{Vertmathbf{b}Vert}leepsilon$ 。令 $k=epsilonkappa(A)$ ，当 $epsilon<frac{1}{kappa(A)}$ 时，下面三个性质成立。

性质一 ： $A+delta A$ 为可逆矩阵。

已知 $A$ 是可逆的，由假设条件可得

$Vert A^{-1}(delta A)VertleVert A^{-1}VertcdotVertdelta AVertleepsilonVert A^{-1}VertcdotVert AVert=k<1$

根据Neumann 无穷级数性质可知 $I+A^{-1}(delta A)$ 是可逆的，且 $Vert(I+A^{-1}(delta A))^{-1}Vertlefrac{1}{1-k}$ (见“ Neumann无穷级数 ”)，故 $A(I+A^{-1}(delta A))= A+delta A$ 亦为可逆矩阵。

性质二 ： $displaystylefrac{Vertmathbf{x}+deltamathbf{x}Vert}{Vertmathbf{x}Vert}lefrac{1+k}{1-k}$

将 $(A+delta A)(mathbf{x}+deltamathbf{x})=mathbf{b}+deltamathbf{b}$ 等号两边左乘 $A^{-1}$ ，

$(I+A^{-1}(delta A))(mathbf{x}+deltamathbf{x})=mathbf{x}+A^{-1}(deltamathbf{b})$

性质一指出 $I+A^{-1}(delta A)$ 是可逆的，就有

$mathbf{x}+deltamathbf{x}=( I+A^{-1}(delta A))^{-1}(mathbf{x}+A^{-1}(deltamathbf{b}))$

利用性质一的不等式，可得

$displaystyle Vertmathbf{x}+deltamathbf{x}VertleVert(I+A^{-1}(delta A))^{-1}Vertcdot(Vertmathbf{x}Vert+epsilonVert A^{-1}VertcdotVertmathbf{b}Vert)lefrac{1}{1-k}left(Vertmathbf{x}Vert+kfrac{Vertmathbf{b}Vert}{Vert AVert} ight)$

但 $Vertmathbf{b}Vert=Vert Amathbf{x}VertleVert AVertcdotVertmathbf{x}Vert$ ，所以

$displaystyle Vertmathbf{x}+deltamathbf{x}Vertlefrac{1}{1-k}(Vertmathbf{x}Vert+kVertmathbf{x}Vert)=frac{1+k}{1-k}Vertmathbf{x}Vert$

性质三 ： $displaystylefrac{Vertdeltamathbf{x}Vert}{Vertmathbf{x}Vert}lefrac{2epsilon}{1-k}kappa(A)$

性质二的关系式可写为 $deltamathbf{x}=A^{-1}(deltamathbf{b})-A^{-1}(delta A)(mathbf{x}+deltamathbf{x})$ ，计算 $deltamathbf{x}$ 长度：

$Vertdeltamathbf{x}VertleVert A^{-1}VertcdotVertdeltamathbf{b}-(delta A)(mathbf{x}+deltamathbf{x})Vert$

利用已知条件与三角不等式，可得

$Vertdeltamathbf{x}VertleepsilonVert A^{-1}VertcdotVertmathbf{b}Vert+epsilonVert A^{-1}VertcdotVert AVertcdotVertmathbf{x}+deltamathbf{x}Vert$

将上式除以 $Vertmathbf{x}Vert$ 并利用性质二，

$displaystylefrac{Vertdeltamathbf{x}Vert}{Vertmathbf{x}Vert}leepsilonkappa(A)frac{Vertmathbf{b}Vert}{Vert AVertcdotVertmathbf{x}Vert}+epsilonkappa(A)frac{Vertmathbf{x}+deltamathbf{x}Vert}{Vertmathbf{x}Vert}leepsilonkappa(A)left(1+frac{1+k}{1-k} ight)=frac{2epsilon}{1-k}kappa(A)$

方程解的相对误差来自 $A$ 和 $mathbf{b}$ 的相对扰动量和 $2epsilon$ ，条件数 $kappa(A)$ 再将此误差放大。换句话说，执行数值计算时，因舍入误差而造成的方程解误差有两个来源：一是线性系统自身的敏感性，以 $kappa(A)$ 表示，另一则是真实误差 $delta A$ 和 $deltamathbf{b}$ 。当我们用LU 分解计算时，受限于电脑的精准度，实际得到的是近似矩阵 $ilde{L}$ 和 $ilde{U}$ ，因此无可避免地使用了错误矩阵 $A+delta A= ilde{L} ilde{U}$ 。英国数学家威尔金森(James H. Wilkinson)证明部分轴元法(partial pivoting，见“ PA=LU分解 ”)确实能够有效控制 $delta A$ ，故舍入误差的主要决定因素为系数矩阵的条件数。感谢威尔金森提供的定心丸，除非碰上病态系统，否则我们大可放心地使用高斯消去法。

本文参考：
Gene H. Golub, Charles F. Van Loan，Matrix Computations，2 ^nd ed.，1989。

____________________________________________________________

1. 病态系统

现在有线性系统： Ax = b，解方程

很容易得到解为： x1 = -100, x2 = -200. 如果在样本采集时存在一个微小的误差，比如，将 A 矩阵的系数 400 改变成 401：

则得到一个截然不同的解： x1 = 40000, x2 = 79800.

当解集 x 对 A 和 b 的系数高度敏感，那么这样的方程组就是病态的 (ill-conditioned).

2. 条件数

那么，如何评价一个方程组是病态还是非病态的呢？在此之前，需要了解矩阵和向量的 norm，这里具体是计算很简单的 infinity norm，即找行绝对值之和最大，举个例子：

infinity norm 具有三角性质：||x+y|| <= ||x|| + ||y||. 理解了这些概念，下面讨论一下衡量方程组病态程度的条件数，首先假设向量 b 受到扰动，导致解集 x 产生偏差：

即有：

同时，由于

综合上面两个不等式：

即得到最终的关系：

如果是矩阵 A 产生误差，同样可以得到：

其中，条件数定义为：

一般来说，方程组解集的精度大概是个十进制的位的误差。比如，IEEE 标准表示的双精度浮点数的有效位是 16 位，如果条件数是 1e+10, 那么得到的结果中只有 6 位是精确的。所以，只有当方程组是良态时，残差 R = Ax - b 才能准确指示解的精度。

3. 病态的由来

自己的看法：

线性系统 Ax = b 为什么会病态？归根到底是由于 A 矩阵列向量线性相关性过大，表示的特征太过于相似以至于容易混淆所产生的。举个例子, 现有一个两个十分相似的列向量组成的矩阵 A：

在二维空间上，这两个列向量夹角非常小。假设第一次检测得到数据 b = [1000, 0]^T, 这个点正好在第一个列向量所在的直线上，解集是 [1, 0]^T。现在再次检测，由于有轻微的误差，得到的检测数据是 b = [1000, 0.001]，这个点正好在第二个列向量所在的直线上，解集是 [0, 1]^T。两次求得到了差别迥异的的解集。

4. 与特征值和 SVD 的关系

特征值

假设 A 的两个单位特征向量是 x1, x2, 根据特征向量的性质：

$Ax_1 = \lambda _1x_1, Ax_2 = \lambda _2x_2$

上述矩阵 A 的特征值和特征向量分别为：

对于平面上的某一个向量 b，可以分解为两个特征向量的线性组合：

把上式带入， $b = mx_1 + nx_2 = \frac{m}{\lambda _1}\lambda _1x_1 + \frac{n}{\lambda _2}\lambda _2x_2 = A(\frac{m}{\lambda _1}x_1 + \frac{n}{\lambda _2}x_2) = Ax$

如果 $\lambda _1$ 远远大于 $\lambda _2$ ，当 b 点在 x1 方向发生移动, m 值改变，解集 x 变化不明显，反之，如果在 x2 方向移动， n 值改变，解集 x 变化非常大！可以看到，特征值对解集起到了一个 scaling 的作用。反过来说，如果一个特征值比其它特征值在数量级上小很多，x在对应特征向量 (x2) 方向上很大的移动才能产生b微小的变化.

2. SVD

SVD 分解：

联系上次学到的 SVD 知识，将 A 分解成三个矩阵的乘积，中间的对角线矩阵也起到了 scaling 的作用。我们按照正向思维来考虑这个问题，现在来了一个解集 x 向量，左乘 A 矩阵等价与左乘 USV^T, x 向量正好等于 V^T 最后一行向量，经过 S 矩阵的 scaling 缩小之后对 b 的影响非常小。也就是说，解集 x 在 V^T 最后一行的行向量方向自由度最大！自由度越大，越不稳定，极端情况是该方向奇异值为 0, 解集可以在该方向取任意值，这也正好对应了矩阵 A 有零特征值， Ax 在对应特征向量的方向上移动不改变 Ax 的值。

在不同的 norm 下，条件数又可以由最大奇异值与最小奇异值之间的比值，或者最大特征值和最小特征值之间比值的绝对值来表示，详情请参考维基百科

最后， A 的条件数究竟等于多少呢？ cond(A) = 2e+06

5. 病态矩阵处理方法

真正的自由是建立在规范的基础上的。病态矩阵解集的不稳定性是由于解集空间包含了自由度过大的方向，解决这个问题的关键就是将这些方向去掉，而保留 scaling 较大的方向，从而把解集局限在一个较小的区域内。在上面的讨论中， A 矩阵的特征向量不一定正交，不适合做新基， SVD 分解正好分解出了正交基，可以选前 k 个 v^T 向量作为正交基。

比如，现在只选取前一个 (0.707, 0.707) 方向作为基，解集局限咋 y = x 这条直线上。直观的解释就是， A 矩阵的两个列向量过于类似，我们就可以将它们等同看待，第一次 b = (1000, 0), 解集是(0.5, 0.5), 第二次 b = (1000, 0.001), 解集还是 (0.5, 0.5).

总结起来，解决 A 病态就是将解集限定在一组正交基空间内，即对于坐标 y，选择 k 个正交基 Zk，解决问题：

这个就是 reduce-rank model. 具体方法有 truncated SVD 和 Krylov subspace method。

参考资料：

(1) ILL-Conditioned Systems

作者：daniel-D 出处：http://www.cnblogs.com/daniel-D/ 欢迎转载或分享，但请务必声明文章出处。

分类: Mathematics

标签: matrix