病态矩阵与条件数

您所在的位置：网站首页 › 矩阵的条件数cond1怎么计算 › 病态矩阵与条件数

病态矩阵与条件数

2024-05-08 12:30| 来源: 网络整理| 查看: 265

本文的阅读等级：高级

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

$\Vert A+B\Vert\le\Vert A\Vert+\Vert B\Vert$

$\Vert AB\Vert\le\Vert A\Vert\cdot\Vert B\Vert$

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

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

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

与原方程式 $A\mathbf{x}=\mathbf{b}$ 相减可得 $A(\delta\mathbf{x})=\delta\mathbf{b}$ ，方程解的误差为 $\delta\mathbf{x}=A^{-1}(\delta\mathbf{b})$ 。利用不等式 $\Vert\delta\mathbf{x}\Vert=\Vert A^{-1}(\delta\mathbf{b})\Vert\le\Vert A^{-1}\Vert\cdot\Vert\delta\mathbf{b}\Vert$ 和 $\Vert\mathbf{b}\Vert=\Vert A\mathbf{x}\Vert\le\Vert A\Vert\cdot\Vert\mathbf{x}\Vert$ ，可得

$\displaystyle\frac{\Vert\delta\mathbf{x}\Vert}{\Vert\mathbf{x}\Vert}=\frac{\Vert A^{-1}(\delta\mathbf{b})\Vert}{\Vert\mathbf{x}\Vert}\le\frac{\Vert A^{-1}\Vert\cdot\Vert\delta\mathbf{b}\Vert}{\Vert\mathbf{x}\Vert}\le\left(\Vert A\Vert\cdot\Vert A^{-1}\Vert\right)\frac{ \Vert\delta\mathbf{b}\Vert}{\Vert\mathbf{b}\Vert}$

另一方面，我们考虑不等式 $\Vert\delta\mathbf{b}\Vert=\Vert A(\delta\mathbf{x})\Vert\le\Vert A\Vert\cdot\Vert\delta\mathbf{x}\Vert$ 和 $\Vert\mathbf{x}\Vert=\Vert A^{-1}\mathbf{b}\Vert\le\Vert A^{-1}\Vert\cdot\Vert\mathbf{b}\Vert$ ，就有

$\displaystyle \frac{\Vert\delta\mathbf{x}\Vert}{\Vert\mathbf{x}\Vert}\ge\frac{\Vert\delta\mathbf{b}\Vert}{\Vert A\Vert\cdot\Vert\mathbf{x}\Vert}\ge\left(\Vert A\Vert\cdot\Vert A^{-1}\Vert\right)^{-1}\frac{\Vert\delta\mathbf{b}\Vert}{\Vert\mathbf{b}\Vert}$

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

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

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

$\displaystyle \kappa(A)^{-1}\frac{\Vert\delta\mathbf{b}\Vert}{\Vert\mathbf{b}\Vert}\le\frac{\Vert\delta\mathbf{x}\Vert}{\Vert\mathbf{x}\Vert}\le\kappa(A)\frac{\Vert\delta\mathbf{b}\Vert}{\Vert\mathbf{b}\Vert}$

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

$\displaystyle \sigma_1=\max_{\Vert\mathbf{x}\Vert=1}\Vert\Sigma\mathbf{x}\Vert=\max_{\Vert\mathbf{x}\Vert=1}\Vert A\mathbf{x}\Vert=\Vert A\Vert_2$

$\displaystyle \sigma_n=\min_{\Vert\mathbf{x}\Vert=1}\Vert\Sigma\mathbf{x}\Vert=\min_{\Vert\mathbf{x}\Vert=1}\Vert A\mathbf{x}\Vert=\frac{1}{\Vert A^{-1}\Vert_2}$

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

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

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

$A^2=AA^T=(U\Sigma V^T)(V\Sigma U^T)=U\Sigma^2U^T$

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

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

又若 $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$ 称为右奇异向量，左右奇异向量的关系如下：

$A\mathbf{v}_j=\sigma_j\mathbf{u}_j,~~~j=1,\ldots,n$

考虑 $\mathbf{b}=\alpha\mathbf{u}_1$ ， $\delta\mathbf{b}=\beta\mathbf{u}_n$ ，方程解和误差可分别表示为

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

$\displaystyle \delta\mathbf{x}=A^{-1}(\delta\mathbf{b})=A^{-1}(\beta\mathbf{u}_n)=\frac{\beta}{\sigma_n}\mathbf{v}_n$

计算向量长度并相除，

$\displaystyle\frac{\Vert\delta\mathbf{x}\Vert}{\Vert\mathbf{x}\Vert}=\left(\frac{\sigma_1}{\sigma_n}\right)\frac{\vert\beta\vert}{\vert\alpha\vert}=\kappa(A)\frac{\Vert\delta\mathbf{b}\Vert}{\Vert\mathbf{b}\Vert}$

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

$\displaystyle\frac{\Vert\delta\mathbf{x}\Vert}{\Vert\mathbf{x}\Vert}=\left(\frac{\sigma_n}{\sigma_1}\right)\frac{\vert\beta\vert}{\vert\alpha\vert}=\kappa(A)^{-1}\frac{\Vert\delta\mathbf{b}\Vert}{\Vert\mathbf{b}\Vert}$

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

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

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

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

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

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

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

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

为简化问题，以下假设 $\frac{\Vert\delta A\Vert}{\Vert A\Vert}\le\epsilon$ ， $\frac{\Vert\delta\mathbf{b}\Vert}{\Vert\mathbf{b}\Vert}\le\epsilon$ 。令 $k=\epsilon\kappa(A)$ ，当 $\epsilon\frac{1}{\kappa(A)}$ 时，下面三个性质成立。

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

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

$\Vert A^{-1}(\delta A)\Vert\le\Vert A^{-1}\Vert\cdot\Vert\delta A\Vert\le\epsilon\Vert A^{-1}\Vert\cdot\Vert A\Vert=k1$

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

性质二： $\displaystyle\frac{\Vert\mathbf{x}+\delta\mathbf{x}\Vert}{\Vert\mathbf{x}\Vert}\le\frac{1+k}{1-k}$

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

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

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

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

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

$\displaystyle \Vert\mathbf{x}+\delta\mathbf{x}\Vert\le\Vert(I+A^{-1}(\delta A))^{-1}\Vert\cdot(\Vert\mathbf{x}\Vert+\epsilon\Vert A^{-1}\Vert\cdot\Vert\mathbf{b}\Vert)\le\frac{1}{1-k}\left(\Vert\mathbf{x}\Vert+k\frac{\Vert\mathbf{b}\Vert}{\Vert A\Vert}\right)$

但 $\Vert\mathbf{b}\Vert=\Vert A\mathbf{x}\Vert\le\Vert A\Vert\cdot\Vert\mathbf{x}\Vert$ ，所以

$\displaystyle \Vert\mathbf{x}+\delta\mathbf{x}\Vert\le\frac{1}{1-k}(\Vert\mathbf{x}\Vert+k\Vert\mathbf{x}\Vert)=\frac{1+k}{1-k}\Vert\mathbf{x}\Vert$

性质三： $\displaystyle\frac{\Vert\delta\mathbf{x}\Vert}{\Vert\mathbf{x}\Vert}\le\frac{2\epsilon}{1-k}\kappa(A)$

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

$\Vert\delta\mathbf{x}\Vert\le\Vert A^{-1}\Vert\cdot\Vert\delta\mathbf{b}-(\delta A)(\mathbf{x}+\delta\mathbf{x})\Vert$

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

$\Vert\delta\mathbf{x}\Vert\le\epsilon\Vert A^{-1}\Vert\cdot\Vert\mathbf{b}\Vert+\epsilon\Vert A^{-1}\Vert\cdot\Vert A\Vert\cdot\Vert\mathbf{x}+\delta\mathbf{x}\Vert$

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

$\displaystyle\frac{\Vert\delta\mathbf{x}\Vert}{\Vert\mathbf{x}\Vert}\le\epsilon\kappa(A)\frac{\Vert\mathbf{b}\Vert}{\Vert A\Vert\cdot\Vert\mathbf{x}\Vert}+\epsilon\kappa(A)\frac{\Vert\mathbf{x}+\delta\mathbf{x}\Vert}{\Vert\mathbf{x}\Vert}\le\epsilon\kappa(A)\left(1+\frac{1+k}{1-k}\right)=\frac{2\epsilon}{1-k}\kappa(A)$

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