线性方程组的敏度分析与消去法的舍入误差分析

向量范数和矩阵范数

向量范数

定义 2.1.1 非负函数 \(\|\cdot\|:\mathbb{R}^n\to\mathbb{R}\) 称为 \(\mathbb{R}^n\) 上的向量范数，如果有

• 正定性： \(\forall x\in\mathbb{R}^n,\ \|x\|\ge 0\) ，等号成立当且仅当 \(x = 0\) .
• 齐次性： \(\forall x\in\mathbb{R}^n,\ \forall \alpha\in\mathbb{R},\ \|\alpha x\| = |\alpha|\|x\|\) .
• 三角不等式： \(\forall x,y\in\mathbb{R}^n,\ \|x+y\|\le \|x\|+\|y\|\) .

容易验证向量范数是连续的.

最常用的向量范数是 \(p\) 范数

\[\|x\|_p = \left(|x_1|^p+\cdots+|x_n|^p\right)^{\frac{1}{p}},\quad p\ge 1 \]

其中最重要的是 \(p=1,2,\infty\) ，即

\[\begin{aligned} &\|x\|_1 = |x_1|+\cdots+|x_n|\\ &\|x\|_2 = \left(|x_1|^2+\cdots+|x_n|^2\right)^{\frac{1}{2}} = \sqrt{x^Tx}\\ &\|x\|_{\infty} = \max\{|x_i|:i=1,\cdots,n\} \end{aligned} \]

分别称为 \(1\) 范数， \(2\) 范数， \(\infty\) 范数，这三个范数的正定性和齐次性是显然的，容易证明 \(\|\cdot\|_1\) 和 \(\|\cdot\|_{\infty}\) 满足三角不等式，而对于 \(\|\cdot\|_2\) 需要用到 Cauchy-Schwartz 不等式

\[|x^Ty| \le \|x\|_2\|y\|_2,\quad x,y\in\mathbb{R}^n \]

这是 Holder 不等式

\[|x^Ty| \le \|x\|_p\|y\|_p,\quad \dfrac{1}{p} + \dfrac{1}{q} = 1 \]

的特殊形式.

定理 2.1.1 设 \(\|\cdot\|_{\alpha},\|\cdot\|_{\beta}\) 是 \(\mathbb{R}^n\) 上任意两个范数，则存在 \(c_1,c_2\in\mathbb{R}\) ，使得

\[\forall x\in\mathbb{R}^n,\quad c_1\|x\|_{\alpha} \le \|x\|_{\beta} \le c_2\|x\|_{\alpha} \]

这一定理的证明较为复杂，它说明所有范数在此意义下都是等价的.

定理 2.1.2 设 \(x_k\in\mathbb{R}^n\) ，则 \(\lim_{k\to\infty}\|x_k-x\| = 0\) 的充要条件为

\[\lim_{k\to\infty}\left|x_i^{(k)}-x_i\right| = 0,\quad i=1,\cdots,n \]

即向量序列的范数收敛等价于其分量收敛.

矩阵范数

我们可以将全体 \(n\times n\) 实矩阵构成的空间看做 \(n^2\) 维空间，然后直接推广向量范数到矩阵。但这样并未考虑到矩阵乘法，因此

定义 2.1.2 非负函数 \(\|\cdot\|:\mathbb{R}^{n\times n}\to\mathbb{R}\) 称为 \(\mathbb{R}^{n\times n}\) 上的矩阵范数，如果有

• 正定性： \(\forall A\in\mathbb{R}^{n\times n},\ \|A\|\ge 0\) ，等号成立当且仅当 \(A = 0\) .
• 齐次性： \(\forall A\in\mathbb{R}^{n\times n},\ \forall \alpha\in\mathbb{R},\ \|\alpha A\| = |\alpha|\|A\|\) .
• 三角不等式： \(\forall A,B\in\mathbb{R}^{n\times n},\ \|A+B\|\le \|A\|+\|B\|\) .
• 相容性： \(\forall A,B\in\mathbb{R}^{n\times n},\ \|AB\|\le \|A\|\|B\|\) .

类似的，有

• \(\mathbb{R}^{n\times n}\) 上的矩阵范数等价.
• 矩阵序列的范数收敛等价于其元素收敛.

定义 2.1.3 若矩阵范数 \(\|\cdot\|_M\) 和向量范数 \(\|\cdot\|_v\) 满足

\[\|Ax\|_v \le \|A\|_M\|x\|_v,\quad A\in\mathbb{R}^{n\times n},\ x\in\mathbb{R}^n \]

则称矩阵范数 \(\|\cdot\|_M\) 与向量范数 \(\|\cdot\|_v\) 相容.

定理 2.1.3 设 \(\|\cdot\|\) 是 \(\mathbb{R}^n\) 上的向量范数，若定义

\[\|A\|_M = \max_{\|x\|=1}\|Ax\|,\quad A\in\mathbb{R}^{n\times n} \]

则 \(\|\cdot\|_M\) 是 \(\mathbb{R}^{n\times n}\) 上的矩阵范数.

证明 我们只证明两点：首先由于 \(\mathcal{D}=\{x\in\mathbb{R}^n:\|x\|=1\}\) 为有界闭集，而 \(\|\cdot\|\) 是连续函数，因此 \(\|Ax\|\) 最大值存在，这样定义的矩阵范数是有意义的；其次，对于任意 \(x\in\mathbb{R}^n,\ x\neq 0\) ，则有

\[\dfrac{\|Ax\|}{\|x\|} = \left\|A\dfrac{x}{\|x\|}\right\| \le \|A\|_M \]

从而有

\[\|Ax\| \le \|A\|_M\|x\|,\quad x\in\mathbb{R}^n \]

利用此性质可以证明 \(\|\cdot\|_M\) 满足矩阵范数的性质.

定理 2.1.4 设 \(A=[a_{ij}]\in\mathbb{R}^{n\times n}\) ，则有

\[\begin{aligned} &\|A\|_1 = \max_{1\le j\le n}\sum_{i=1}^n|a_{ij}|\\ &\|A\|_{\infty} = \max_{1\le i\le n}\sum_{j=1}^n|a_{ij}|\\ &\|A\|_2 = \sqrt{\lambda_{\max}(A^TA)} \end{aligned} \]

其中 \(\lambda_{\max}(A^TA)\) 表示 \(A^TA\) 的最大特征值.

证明 \(A=0\) 时显然成立；对于 \(1\) 范数，我们从结论出发，发现 \(\|A\|_1\) 是它所有列向量的 \(1\) 范数的最大值，因此我们取

\[A= [a_1,\cdots,a_n],\quad \delta = \max_{1\le j\le n}\|a_j\|_1 \]

则任意 \(x\in\mathbb{R}^n,\ \|x\|_1=1\) 有

\[\|Ax\|_1 = \left\|\sum_{j=1}^nx_ja_j\right\|_1 \le \sum_{j=1}^n|x_j|\|a_j\|_1 \le \left(\sum_{j=1}^n|x_j|\right)\max_{1\le j\le n}\|a_j\|_1 = \delta \]

设 \(\delta\) 是第 \(j_0\) 列的范数，可以取 \(e_{j_0}\) ，这样

\[\|e_{j_0}\| = 1,\quad \|Ae_{j_0}\| = \delta \]

因此可以取得最大值 \(\delta\) .

对于 \(\infty\) 范数，我们从结论出发，发现 \(\|A\|_{\infty}\) 是它的每行元素绝对值和的最大值，可以这样考虑：

\[\|Ax\|_{\infty} = \left\|\left(\sum_{j=1}^na_{1j}x_j,\cdots,\sum_{j=1}^na_{nj}x_j\right)^T\right\|_{\infty} \]

因为 \(x_j\) 的符号可以根据 \(a_{ij}\) 调整，使得 \(a_{ij}x_j\ge 0\) ，这样我们只需要取 \(\sum_{j=1}^na_{ij}x_j\) 中的最大值，记为 \(\eta\) ， 这意味着

\[\|Ax\|_{\infty} \le \max_{1\le i\le n}\sum_{j=1}^na_{ij}x_j \le \max_{1\le i\le n}\sum_{j=1}^n|a_{ij}| = \eta \]

设最大值在 \(i=i_0\) 处取得，则我们取 \(x = (\mathrm{sgn}(a_{i_0,1}),\cdots,\mathrm{sgn}(a_{i_0,n}))^T\) ，则有

\[\|x\|_{\infty} = 1,\quad \|Ax\|_{\infty} = \eta \]

因此可以取得最大值 \(\eta\) .

对于 \(2\) 范数，应用向量范数有

\[\|A\|_2 = \max_{\|x\|_2=1}\|Ax\|_2 = \max_{\|x\|_2=1}\left(x^TA^TAx\right)^{\frac{1}{2}} \]

注意到 \(A^TA\) 对称半正定，因而有特征值

\[\lambda_1 \ge \lambda_2\ge \cdots \ge \lambda_n \ge 0 \]

和对应的正交标准特征向量 \(v_1,v_2,\cdots,v_n\in\mathbb{R}^n\) ，则任意 \(\|x\|_2=1,\ x\in\mathbb{R}^n\) 有

\[x = \sum_{i=1}^n\alpha_iv_i,\quad \sum_{i=1}^n\alpha_i^2 = 1 \]

从而有

\[x^TA^TAx = \sum_{i=1}^n\lambda_i\alpha_i^2 \le \lambda_1 \]

另一方面，取 \(x=v_1\) 有

\[x^TA^TAx = v_1^TA^TAv_1 = v_1^T\lambda_1v_1 = \lambda_1 \]

因此可以取得最大值 \(\lambda_1\) .

定理 2.1.5 设 \(A\in\mathbb{R}^{n\times n}\) ，则

• \(\|A\|_2 = \max\left\{\left|y^TAx\right|:x,y\in\mathbb{C}^n,\ \|x\|_2=\|y\|_2=1\right\}\)

• \(\|A^T\|_2 = \|A\|_2 = \sqrt{\|A^TA\|_2}\)

• 对任意 \(n\) 阶正交阵 \(U,V\) ，有 \(\|UA\|_2=\|AV\|_2 = \|A\|_2\)

证明 我们先证明小于关系，再证明可以取得最值

\[\left|y^TAx\right| \le \left|y^T\right|\|Ax\|_2 = \|Ax\|_2 \]

设 \(x_0\in\mathbb{R}^n,\ \|x_0\|_2 = 1,\ \|A\|_2 = \|Ax_0\|_2\) ，取 \(y=Ax_0/\|A\|_2,\ \|y\|_2=1\) ，则有

\[\left|y^TAx_0\right| = \left|(Ax_0)^TAx_0\right|/\|A\|_2 = \|A\|_2 \]

从而可以取得最大值 \(\|A\|_2\) .

在 \(\mathbb{R}^{n\times n}\) 上另一个常用的矩阵范数为

\[\|A\|_F = \left(\sum_{i,j=1}^n|a_{ij}|^2\right)^{\frac{1}{2}} \]

称为 Frobenius 范数，它是向量 \(2\) 范数的自然推广.

定义 2.1.4 设 \(A\in\mathbb{R}^{n\times n}\) ，则称

\[\rho(A) = \max\{|\lambda|:\lambda\in\lambda(A)\} \]

为 \(A\) 的谱半径，这里 \(\lambda(A)\) 表示 \(A\) 的特征值全体.

定理 2.1.6 设 \(A\in\mathbb{R}^{n\times n}\) ，则有

• 对 \(\mathbb{C}^{n\times n}\) 上的任意矩阵范数 \(\|\cdot\|\) 有

\[\rho(A) \le \|A\| \]

• 对任意 \(\epsilon>0\) ，存在 \(\mathbb{C}^{n\times n}\) 上的算子范数 \(\|\cdot\|\) 使得

\[\|A\|\le \rho(A) + \epsilon \]

这说明谱半径是所有矩阵范数的下确界.

证明 注意到 \(\rho(A)\) 是特征值 \(\lambda\) 的绝对值，我们取特征向量 \(x\in\mathbb{C}^{n\times n}\) ，则有

\[\rho(A)\left\|xe_1^T\right\| = \left\|\lambda xe_1^T\right\| = \left\|Axe_1^T\right\| \le \|A\|\left\|xe_1^T\right\| \]

从而有

\[\rho(A) \le \|A\| \]

第二部分应用若当标准型定义满足条件的算子范数，事实上有非奇异矩阵 \(P\in\mathbb{C}^{n\times n}\) 使得

\[P^{-1}AP = \left( \begin{matrix} \lambda_1 & \delta_2\\ & \lambda_2 & \delta_3\\ & & \ddots & \ddots\\ & & & \lambda_{n-1} & \delta_{n}\\ & & & & \lambda_n \end{matrix} \right) \]

其中 \(\delta_i = 0,1\) ，即为若当标准型。对任意 \(\epsilon>0\) ，我们作微小扰动 \(D_{\epsilon} = \mathrm{diag}(1,\epsilon,\cdots,\epsilon^{n-1})\) ，则有

\[D_{\epsilon}^{-1}P^{-1}APD_{\epsilon} = \left( \begin{matrix} \lambda_1 & \epsilon\delta_2\\ & \lambda_2 & \epsilon\delta_3\\ & & \ddots & \ddots\\ & & & \lambda_{n-1} & \epsilon\delta_{n}\\ & & & & \lambda_n \end{matrix} \right) \]

定义范数

\[\|G\|_{\epsilon} = \left\|D_{\epsilon}^{-1}P^{-1}GPD_{\epsilon}\right\|_{\infty},\quad G\in\mathbb{C}^{n\times n} \]

则容易验证这是矩阵范数，并且

\[\|A\|_{\epsilon} = \left\|D_{\epsilon}^{-1}P^{-1}APD_{\epsilon}\right\|_{\infty} = \max_{1\le i\le n}(|\lambda_i|+|\epsilon\delta_i|) \le \rho(A) + \epsilon \]

其中假定 \(\delta_0 = 0\) .

定理 2.1.7 设 \(A\in\mathbb{C}^{n\times n}\) ，则

\[\lim_{k\to\infty} A^k = 0\ \Leftrightarrow\ \rho(A) < 1 \]

证明 首先，如果 \(A^k\to 0\) ，考虑到 \(\lambda\in\lambda(A)\) 有 \(\lambda^k\in\lambda\left(A^k\right)\) ，则

\[\rho(A)^k \le \rho\left(A^k\right) \le \left\|A^k\right\|_2 \to 0 \]

从而只能有 \(\rho(A)<1\) ；反之，取 \(\epsilon = 1-\rho(A)>0\) ，则有算子范数 \(\|\cdot\|\) 使得

\[\|A\| \le \rho(A) + \dfrac{\epsilon}{2} = 1 \]

从而有

\[\left\|A^k\right\| \le \|A\|^k \to 0 \]

即证.

定理 2.1.8 设 \(A\in\mathbb{C}^{n\times n}\) ，则

• \(\sum_{k=0}^{\infty}A^k\) 收敛当且仅当 \(\rho(A)<1\) .
• 若 \(\sum_{k=0}^{\infty}A^k\) 收敛，则

\[\sum_{k=0}^{\infty}A^k = (I-A)^{-1} \]

且存在 \(\mathbb{C}^{n\times n}\) 上的算子范数 \(\|\cdot\|\) 使得

\[\left\|(I-A)^{-1}-\sum_{k=0}^mA^k\right\| \le\dfrac{\|A\|^{m+1}}{1-\|A\|} \]

对任意 \(m\in\mathbb{N}\) 成立.

证明 先计算矩阵和

\[S_n = \sum_{k=0}^nA^k\ \Rightarrow\ (I-A)S_n = I - A^{n+1}\ \Rightarrow\ S_n = (I-A)^{-1}\left(I - A^{n+1}\right) \]

则第一部分显然；当 \(S_n\) 收敛，必有 \(S = (I-A)^{-1}\) 存在，并且

\[\left\|(I-A)^{-1}-\sum_{k=0}^mA^k\right\| = \|S-S_m\| = \left\|(I-A)^{-1}A^{m+1}\right\| \le \dfrac{\|A\|^{m+1}}{1-\|A\|} \]

即证.

推论 2.1.1 设 \(\mathbb{C}^{n\times n}\) 上的算子范数 \(\|\cdot\|\) 满足 \(\|I\|=1\) ，且 \(A\in\mathbb{C}^{n\times n}\) 有 \(\|A\|<1\) ，则 \(I-A\) 可逆，且有

\[\left\|(I-A)^{-1}\right\| \le\dfrac{1}{1-\|A\|} \]

证明 由 \(\rho(A)\le\|A\|<1\) ，则 \((I-A)^{-1}\) 存在，满足上面定理条件，令 \(m=0\) 则有

\[\left\|(I-A)^{-1}\right\|-\|I\| \le \left\|(I-A)^{-1}-I\right\| \le \dfrac{\|A\|}{1-\|A\|} \]

从而有

\[\left\|(I-A)^{-1}\right\| \le 1 + \dfrac{\|A\|}{1-\|A\|} = \dfrac{1}{1-\|A\|} \]

即证.

线性方程组的敏度分析

考虑非奇异线性方程组 \(Ax=b\) 对系数扰动的敏感性，设有

\[(A+\delta A)(x+\delta x) = b + \delta b \]

则代入 \(b = Ax\) 得

\[(A+\delta A)\delta x = \delta b - \delta Ax \]

由于 \(A\) 非奇异，因此 \(A+\delta A\) 非奇异，这是由 \(\det\) 的连续性得到的；事实上，只要 $ \left|A^{-1}\right||\delta A|<1 $ ，就有 \(A+\delta A\) 可逆，因为

\[A+\delta A = A\left(I+A^{-1}\delta A\right) \]

即有 \(A+\delta A\) 可逆等价于 \(I+A^{-1}\delta A\) 可逆；而且

\[\left\|\left(I+A^{-1}\delta A\right)^{-1}\right\| \le \dfrac{1}{1-\left\|A^{-1}\right\|\|\delta A\|} \]

因此此时 \(A+\delta A\) 可逆；我们考虑

\[\delta x = \left(A+\delta A\right)^{-1}(\delta b-\delta Ax) = \left(I+A^{-1}\delta A\right)^{-1}A^{-1}(\delta b-\delta Ax) \]

取范数就有

\[\|\delta x\| \le \left\|\left(I+A^{-1}\delta A\right)^{-1}\right\|\left\|A^{-1}\right\|(\|\delta b\|+\|\delta A\|\|x\|) \le \dfrac{\left\|A^{-1}\right\|}{1-\left\|A^{-1}\right\|\|\delta A\|}(\|\delta b\|+\|\delta A\|\|x\|) \]

注意到 \(\|b\|\le \|A\|\|x\|\) ，则有

\[\dfrac{\|\delta x\|}{\|x\|} \le \dfrac{\left\|A^{-1}\right\|\|A\|}{1-\left\|A^{-1}\right\|\|\delta A\|}\left(\dfrac{\|\delta A\|}{\|A\|}+\dfrac{\|\delta b\|}{\|b\|}\right) \]

从而我们有如下定理：

定理 2.2.1 设 \(\|\cdot\|\) 是 \(\mathbb{R}^{n\times n}\) 上的一个满足条件 \(\|I\|=1\) 的矩阵范数，且 \(A\in\mathbb{R}^{n\times n}\) 可逆， \(b\in\mathbb{R}^n\) 非零；再假定 \(\delta A\in\mathbb{R}^{n\times n}\) 满足 $ \left|A^{-1}\right||\delta A|<1 $ ，若 \(x\) 和 \(x+\delta x\) 分别是线性方程组

\[Ax = b,\quad (A+\delta A)(x+\delta x) = b + \delta b \]

的解，则

\[\dfrac{\|\delta x\|}{\|x\|} \le \dfrac{\kappa(A)}{1-\kappa(A)\dfrac{\|\delta A\|}{\|A\|}}\left(\dfrac{\|\delta A\|}{\|A\|}+\dfrac{\|\delta b\|}{\|b\|}\right) \]

其中 \(\kappa(A)=\left\|A^{-1}\right\|\|A\|\) .

? 当 \(\|\delta A\|/\|A\|\) 较小时，有

\[\dfrac{\kappa(A)}{1-\kappa(A)\dfrac{\|\delta A\|}{\|A\|}} \approx \kappa(A) \]

从而有

\[\dfrac{\|\delta x\|}{\|x\|} \lessapprox \kappa(A)\left(\dfrac{\|\delta A\|}{\|A\|}+\dfrac{\|\delta b\|}{\|b\|}\right) \]

因此相对误差由 \(\kappa(A)\) 和初始条件的相对扰动限制.

定义 2.2.1 称数 \(\kappa(A)=\|A\|\left\|A^{-1}\right\|\) 为线性方程组 \(Ax = b\) 的条件数.

由于范数的等价性，因此不同范数下的条件数等价.

推论 2.2.1 设 \(\|\cdot\|\) 是 \(\mathbb{R}^{n\times n}\) 上的一个满足条件 \(\|I\|=1\) 的矩阵范数，且 \(A\in\mathbb{R}^{n\times n}\) 可逆， \(b\in\mathbb{R}^n\) 非零；再假定 \(\delta A\in\mathbb{R}^{n\times n}\) 满足 $ \left|A^{-1}\right||\delta A|<1 $ ，则有

\[\dfrac{\left\|\left(A+\delta A\right)^{-1}-A^{-1}\right\|}{\left\|A^{-1}\right\|} \le \dfrac{\kappa(A)}{1-\kappa(A)\dfrac{\|\delta A\|}{\|A\|}}\dfrac{\|\delta A\|}{\|A\|} \]

证明 利用之前得到的

\[\left\|(I-A)^{-1}-I\right\| \le \dfrac{\|A\|}{1-\|A\|} \]

我们有

\[\left\|(A+\delta A)^{-1}-A^{-1}\right\| \le \left\|\left(I+A^{-1}\delta A\right)^{-1}-I\right\|\left\|A^{-1}\right\| \le \dfrac{\left\|A^{-1}\delta A\right\|}{1-\left\|A^{-1}\delta A\right\|}\left\|A^{-1}\right\| \le \dfrac{\left\|A^{-1}\right\|\|\delta A\|}{1-\left\|A^{-1}\right\|\|\delta A\|}\left\|A^{-1}\right\| \]

两边同时除以 \(\left\|A^{-1}\right\|\) 得到

\[\dfrac{\left\|\left(A+\delta A\right)^{-1}-A^{-1}\right\|}{\left\|A^{-1}\right\|} \le \dfrac{\kappa(A)}{1-\kappa(A)\dfrac{\|\delta A\|}{\|A\|}}\dfrac{\|\delta A\|}{\|A\|} \]

即证.

定理 2.2.2 设 \(A\in\mathbb{R}^{n\times n}\) 可逆，则

\[\min\left\{\dfrac{\|\delta A\|_2}{\|A\|_2}:\exist(A+\delta A)^{-1}\right\} = \dfrac{1}{\|A\|_2\left\|A^{-1}\right\|_2} = \dfrac{1}{\kappa_2(A)} \]

这说明当 \(A\) 十分病态时，它与一个奇异矩阵十分靠近.

证明 由前容易发现当 $ \left|A^{-1}\right||\delta A|<1 $ 时， \(A+\delta A\) 可逆，因而 \(1/\kappa_2(A)\) 是一个下界，我们说明可以取到这个下界：取 \(x\in\mathbb{R}^n\) ，使得 \(\left\|A^{-1}x\right\|_2 = \left\|A^{-1}\right\|_2,\ \|x\|_2=1\) ，令

\[y = \dfrac{A^{-1}x}{\left\|A^{-1}x\right\|_2},\quad \delta A = -\dfrac{xy^T}{\left\|A^{-1}\right\|_2} \]

则有 \(\|y\|_2=1\) ，且

\[(A+\delta A)y = \dfrac{x}{\left\|A^{-1}x\right\|_2} - \dfrac{x}{\left\|A^{-1}\right\|_2} = 0\\ \|\delta A\|_2 = \max_{\|z\|_2=1}\left\|\dfrac{xy^T}{\left\|A^{-1}\right\|_2}z\right\|_2 = \dfrac{\|x\|_2}{\left\|A^{-1}\right\|_2}\max_{\|z\|_2=1}\left|y^Tz\right| = \dfrac{1}{\left\|A^{-1}\right\|_2} \]

因而我们取到了这个下界，即证.

基本运算的舍入误差分析

关于浮点数系统的部分我们在数值分析中进行讨论，在这里只列出相关定理.

定理 2.3.1 设 \(m\le \|x\|\le M\) ，其中 \(m,M\) 是有效数字，则

\[\mathrm{fl}(x) = x(1+\delta),\quad |\delta| \le \epsilon_u \]

其中 \(\epsilon_u\) 为舍入单位，

\[\epsilon_u = \left\{ \begin{aligned} &\dfrac{1}{2}\beta^{1-t},\quad &舍入法\\ &\beta^{1-t},\quad &截断法 \end{aligned} \right. \]

这里利用相对误差分析即可.

定理 2.3.2 设 \(a,b\in\mathcal{F}\) ，则

\[\mathrm{fl}(a\odot b) = (a\odot b)(1+\delta),\quad |\delta|\le \epsilon_u \]

四则运算的误差界.

定理 2.3.3 若 \(|\delta_i|\le \epsilon_u,\ n\epsilon_u\le 0.01\) ，则

\[1-n\epsilon_u \le \prod_{i=1}^n(1+\delta_i) \le 1+1.01n\epsilon_u \]

此定理在数值分析中有详细证明.

列主元 Gauss 消去法的舍入误差分析

这里进行向后误差分析：首先考虑三角分解的舍入误差

引理 2.4.1 设 \(n\times n\) 浮点数矩阵 \(A=[a_{ij}]\) 有三角分解，且 \(1.01n\epsilon_u\le 0.01\) ，则三角分解满足

\[\widetilde{L}\widetilde{U} = A + E,\quad |E| \le 2.05n\epsilon_u|\widetilde{L}||\widetilde{U}| \]

分析较为繁琐，我们省略这部分.

注意到交换矩阵的行或列不会引入误差，因此有

推论 2.4.1 设 \(n\times n\) 浮点数矩阵 \(A\) 非奇异，且 \(1.01n\epsilon_u\le 0.01\) ，则列主元 Gauss 消去法计算得到的分解有

\[\widetilde{L}\widetilde{U} = \widetilde{P}A + E,\quad |E| \le 2.05n\epsilon_u|\widetilde{L}||\widetilde{U}| \]

这是直接推论，接着就是求解三角形方程组.

引理 2.4.2 设 \(n\times n\) 浮点数三角矩阵 \(S\) 非奇异，且 \(1.01n\epsilon_u\le 0.01\) ，则三角形方程组 \(Sx=b\) 的计算解 \(\widetilde{x}\) 满足

\[(S+H)\widetilde{x} = b,\quad |H| \le 1.01n\epsilon_u|S| \]

这种证明看一看就好.

应用上述引理到三角形方程组

\[\widetilde{L}y = \widetilde{P}b,\quad \widetilde{U}x = y \]

得到计算解满足

\[(\widetilde{L}+F)(\widetilde{U}+G)\widetilde{x} = \widetilde{P}b,\quad |F|\le 1.01n\epsilon_u|\widetilde{L}|,\ |G|\le 1.01n\epsilon_u|\widetilde{U}| \]

通过代入估计得到

\[|\delta A|\le 4.09n\epsilon_u\widetilde{P}^T|\widetilde{L}||\widetilde{U}| \]

注意到 \(\widetilde{L}\) 元素绝对值不大于 \(1\) ，则 \(\|\widetilde{L}\|_{\infty}\le n\) ；为了给出 \(\|\widetilde{U}\|_{\infty}\) 的估计，我们定义

\[\rho = \max_{i,j}|\widetilde{u}_{ij}|/\max_{i,j}|a_{ij}| \]

称为列主元 Gauss 消去法的增长因子.

定理 2.4.1 设 \(n\times n\) 浮点数矩阵 \(A\) 非奇异，且 \(1.01n\epsilon_u\le 0.01\) ，则列主元 Gauss 消去法的计算解满足

\[(A+\delta A)\widetilde{x} = b,\quad \|\delta A\|_{\infty}/\|A\|_{\infty} \le 4.09n^3\rho\epsilon_u \]

这是上面过程的直接推论，从而我们可以由 \(A\) 的敏度分析来得到解的敏度.

计算解的精度估计和迭代改进

假设一种计算线性方程组 \(Ax=b\) 的算法得到的计算解为 \(\hat{x}\) ，令

\[r = b - A\hat{x}\ \Rightarrow\ r = Ax - A\hat{x} = A(x-\hat{x}) \]

从而有

\[\|x-\hat{x}\| = \left\|A^{-1}r\right\| \le \left\|A^{-1}\right\|\|r\| \]

又有

\[\|b\| \le \|A\|\|x\|\ \Rightarrow\ \dfrac{1}{\|x\|} \le \dfrac{\|A\|}{\|b\|} \]

我们得到

\[\dfrac{\|x-\hat{x}\|}{\|x\|} \le \dfrac{\left\|A^{-1}\right\|\|r\|}{\|x\|} \le \left\|A^{-1}\right\|\|A\|\dfrac{\|r\|}{\|b\|} \]

在上式中取 \(\infty\) 范数有

\[\dfrac{\|x-\hat{x}\|_{\infty}}{\|x\|_{\infty}} \le \kappa_{\infty}(A)\dfrac{\|r\|_{\infty}}{\|b\|_{\infty}} \]

这样我们可以通过计算右式来进行精度估计.

注意到 \(\left\|A^{-1}\right\|_{\infty}\) 是较难计算的，因此我们给出 LAPACK 采用的一种优化方法：

盲人爬山法

? 设 \(B\in\mathbb{R}^{n\times n}\) ，我们估计 \(\|B\|_1\) ，定义

\[f(x) = \|Bx\|_1,\quad \mathcal{D} = \{x\in\mathbb{R^{n}}:\|x\|_1\le 1\} \]

则可以证明 \(f\) 是凸函数， \(\mathcal{D}\) 是凸集，并且 \(\|B\|_1\) 就是 \(f\) 在 \(\mathcal{D}\) 上的最大值.

设 \(f\) 在 \(x\) 处的梯度 \(\nabla f(x)\) 存在，则由凸函数的性质有

\[f(y) \ge f(x) + \nabla f(x)(y-x),\quad y\in\mathbb{R}^n \]

即凸函数总是在它的切线上方；若有 \(x_0\in\mathbb{R}^n,\ \|x_0\|_1=1\) ，使得 \(Bx_0\) 的任意分量非零，令

\[\xi_i = \mathrm{sgn}\left(\sum_{j=1}^nb_{ij}x_j\right) \]

则在 \(x_0\) 附近将 \(f(x)\) 表示为

\[f(x) = \sum_{i=1}^n\sum_{j=1}^n\xi_ib_{ij}x_j \]

因此有

\[\nabla f(x_0) = \left(\dfrac{\partial f(x_0)}{\partial x_1},\cdots,\dfrac{\partial f(x_0)}{\partial x_n}\right) = v^TB = \left(B^Tv\right)^T,\quad v = (\xi_1,\cdots,\xi_n)^T \]

这里 \(v\) 是 \(\omega\) 的符号向量；再令

\[\omega = Bx_0,\quad z = B^Tv \]

分别代表 \(f(x_0)\) 和 \(\nabla f(x_0)\) 的转置.

定理 2.5.1 根据上面得到的关系，我们有

• 若 \(\|z\|_{\infty}\le z^Tx_0\) ，则 \(\|\omega\|_1=\|Bx_0\|_1\) 是局部极大值
• 若 \(\|z\|_{\infty} > z^Tx_0\) ，则 \(\|Be_j\|_1 > \|Bx_0\|_1\) ，其中 \(j\) 满足

\[|z_j| = \|z\|_{\infty} \]

证明 在 \(x_0\) 附近可将 \(f\) 看作线性函数

\[f(x) = \|\omega\|_1 + z^T(x-x_0) \]

若 \(\|z\|_{\infty}\le z^Tx_0\) ， 对于\(\|x\|_1\le 1\) ，有

\[z^T(x-x_0) = z^Tx-z^Tx_0 \le \|z\|_{\infty}\|x\|_1 - z^Tx_0 \le \|z\|_{\infty} - z^Tx_0 \le 0 \]

因而在 \(x_0\) 附近 \(f\) 都不大于 \(\|\omega\|_1\) ；反之，对于满足的 \(j\) 有

\[\begin{aligned} f(e_j) &\ge \|\omega\|_1 + z^T(e_j-x_0) = \|\omega\|_1 + z^Te_j-z^Tx_0 \\ &= \|\omega\|_1 + |z_j| - z^Tx_0\\ &= \|\omega\|_1 + \|z\|_{\infty} - z^Tx_0 > \|\omega\|_1 \end{aligned} \]

即证.

我们可以应用该定理设计算法求解 \(\left\|A^{-1}\right\|_{\infty}\) 的估计值：设有分解 \(PA = LU\) ，只需计算 \(\|B\|_{1},\ B = A^{-T}\) 即可，此时

\[\omega = Bx,\quad z = B^Tv\ \Leftrightarrow\ A^T\omega = x,\quad Az = v \]

因此，给定初始的 \(x_0\) ，可以解出 \(\omega\) ，然后计算 \(\omega\) 的符号向量 \(v\) ，解出梯度转置 \(z\) ，根据上述定理判断是否达到最大值，若未达到最大值就调整 \(x\) ，循环计算直到达到最大值.

注意到一共只有 \(n\) 种 \(e_j\) ，因此循环必然在有限步后终止.

迭代改进

如果计算解 \(\hat{x}\) 精度太低，同时 \(A\) 不是非常病态，则可应用 Newton 迭代法改进精度

• 计算 \(r = b-A\hat{x}\)
• 求解 \(Az = r\)
• 计算 \(x = \hat{x} + z\)

重复迭代直到相对误差足够小.