高斯消去法的误差分析

两步法获得 $Ax=b$ 解的误差界：

1) 向后误差分析

求解线性代数方程组

$Ax = b$

设计一个算法 $alg(A, b)$, 求出一个数值解

$$\hat x = alg(A, b)$$

数值解 $\hat x$，相当于求解下面扰动后的方程：

$(A + \delta A)\hat x = b + \delta b$

其中 $\delta A$ 和 $\delta b$ 称为算法的 向后误差。

而 $\delta x = \hat{x} - x$ 的界是我们希望界定的 计算误差。

2）利用扰动分析理论分析, 如前面得到的 ( 2.2 ) 和 ( 2.5 )。

$$\frac{ \| \delta x \|}{ \| \hat x \|} \leq \| A^{-1} \|\cdot \| A \|\cdot \left(\frac{ \| \delta A \|}{ \| A \|} + \frac{ \| \delta b \| }{ \| A \| \| \hat x \| } \right)\quad ( 2.3 ).$$

• (2.3) 算出的误差界往往太大。

$$\| \delta x \| = \| A^{-1}r \|\leq \| A^{-1} \|\cdot \| A\hat{x} - b \|.\quad ( 2.5 )$$

• ( 2.5 ) 更易计算, 在实践中更常使用。

本节的两个目标：

• 展示如何实施高斯消去法，使得向后误差 $\delta A$ 和 $\delta b$ 尽量小：

$$\begin{align} \frac{ \| \delta A \| }{ \| A \| } = O(\varepsilon), \quad \frac{ \| \delta b \| }{ \| b \| } = O(\varepsilon). \end{align}$$

• 获得同时易于计算(cheapness)和紧致(tightness)的误差界。
  • 紧致就是说算出的误差界尽量接近真实的误差。
  • 易于计算和紧致的误差界可能不存在。

为什么要选主元？

给定一个矩阵, 这里用 4 位的十进制的浮点数表示，

$$A = \begin{bmatrix} 0.0001 & 1 \\ 1 & 1 \\ \end{bmatrix}$$

易知 $\| A \|_\infty \| A^{-1} \|_\infty \approx 4$。

$$L = \begin{bmatrix} 1 & 0 \\ fl(\frac{1}{10^{-4}}) & 1 \end{bmatrix}$$ 其中 $fl(\frac{1}{10^{-4}}) = 10^4$

$$U = \begin{bmatrix} 10^{-4} & 1 \\ & fl(1 - 10^4*1) \end{bmatrix}$$ 其中 $fl(1 - 10^4*1) = -10^4$ ( 大数吃小数 )

$$LU = \begin{bmatrix} 1 & 0 \\ 10^4 & 1 \\ \end{bmatrix} \begin{bmatrix} 10^{-4} & 1 \\ & -10^4 \end{bmatrix} = \begin{bmatrix} 10^{-4} & 1 \\ 1 & 0 \end{bmatrix} \not=\begin{bmatrix} 0.0001 & 1 \\ 1 & 1 \\ \end{bmatrix}$$

• $a_{22}$ 换成其它小的数, $LU$ 分解的结果是一样的， 因为 $fl(a_{22} - 10^4*1) = -10^4$。
• 这种现象称为 数值不稳定 ( numerical instability )。
• 取 $b = [1, 2]^T$， 用 $LU$ 分解求解 $Ax=b$。
• $\| A-LU \|$ 太大。
• $\kappa(A) \approx 4$， 而 $\kappa(LU) \approx 10^8$
• 列主元可以解决数值不稳定的问题。

$$L = \begin{bmatrix} 1 & 0 \\ fl( .0001 / 1) & 1 \end{bmatrix} \begin{bmatrix} 1 & 0 \\ 0.0001 & 1 \end{bmatrix}$$

$$U = \begin{bmatrix} 1 & 1 \\ 0 & fl(1 - .0001\cdot 1) \end{bmatrix} \begin{bmatrix} 1 & 1 \\ 0 & 1 \end{bmatrix} .$$

注意 如果 $LU$ 分解过程中出现很大的量，则 $A$ 中的元素减去这些量时， 就会出现大数吃小数的现象。

正式的高斯消去法误差分析

因为浮点数运算中舍入误差的存在， 用 $LU$ 分解算法来计算 $Ax=b$ 会带来向 后误差 $\delta A$ ，这个向后误差有两部分组成

1. $LU$ 分解算法带来的误差， 即 $A = LU + E$。
2. 求解 $LU x = b$ 带来的向后误差， 即 $( L+\delta L ) ( U + \delta U ) \hat{x} = b$。

本小节就来分析这两种误差的界。

首先，分析 $A$ 的元素 $a_{jk}$ 与 $L\cdot U$ 的元素 $$\sum_{i=1}^{k} l_{ji}u_{ik}$$ 的差别。

首先，给定一个一般的 $LU$ 分解 $$\begin{align} &\begin{bmatrix} a_{11} & a_{12} & a_{13} & \cdots & a_{1n} \\ a_{21} & a_{22} & a_{23} & \cdots & a_{2n} \\ a_{31} & a_{32} & a_{33} & \cdots & a_{3n} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ a_{n1} & a_{n 2} & a_{n 3} & \cdots & a_{nn} \end{bmatrix}\\ =& \begin{bmatrix} 1 & & & & \\ l_{21} & 1 & & & \\ l_{31} & l_{32} & 1 & & \\ \vdots & \vdots & \vdots & \ddots & \\ l_{n 1} & l_{n 2 } & l_{n 3} & \cdots & 1 \end{bmatrix} \begin{bmatrix} u_{11} & u_{12} & u_{13} & \cdots & u_{1n} \\ & u_{22} & u_{23} & \cdots & u_{2n} \\ & & u_{33} & \cdots & u_{3n} \\ & & & \ddots & \vdots \\ & & & & u_{nn} \end{bmatrix} \end{align}$$

上述三个矩阵的元素有如下的关系成立：

$$u_{j k}=a_{j k}-\sum_{i=1}^{j-1} l_{j i} u_{i k} \text{ with } j < k\quad ( 1 ).$$

$$l_{j k}=\frac{a_{j k}-\sum_{i=1}^{k-1} l_{j i} u_{i k}}{u_{k k}} \text{ with } j \geq k \quad ( 2 )$$

两个向量做内积运算，浮点运算结果满足如下的关系式 ( Question 1.10 )：

$$\mathrm{f}\left(\sum_{i=1}^{d} x_{i} y_{i}\right)=\sum_{i=1}^{d} x_{i} y_{i}\left(1+\delta_{i}\right) \quad \text { with }\left|\delta_{i}\right| \leq d \varepsilon$$

应用到关系式 ( 1 )， 可得：

$$u_{j k}=\left(a_{j k}-\sum_{i=1}^{j-1} l_{j i} u_{i k}\left(1+\delta_{i}\right)\right)\left(1+\delta^{\prime}\right)$$

其中 $|\delta_i| \leq ( j-1 )\varepsilon$, $|\delta'| \leq \varepsilon$ 。

由上式解出 $a_{jk}$ $$\begin{aligned} a_{j k}=& \frac{1}{1+\delta^{\prime}} u_{j k} \cdot l_{j j}+ \sum_{i=1}^{j-1} l_{j i} u_{i k}\left(1+\delta_{i}\right) \quad \text { since } l_{j j}=1 \\ =& \sum_{i=1}^{j} l_{j i} u_{i k}+\sum_{i=1}^{j} l_{j i} u_{i k} \delta_{i} \\ & \text { with }\left|\delta_{i}\right| \leq(j-1) \varepsilon \text { and } 1+\delta_{j} \equiv \frac{1}{1+\delta^{\prime}} \\ & \equiv \sum_{i=1}^{j} l_{j i} u_{i k}+E_{j k} \end{aligned}$$

从而可以界定 $E_{jk}$ $$\left|E_{j k}\right|=\left|\sum_{i=1}^{j} l_{j i} \cdot u_{i k} \cdot \delta_{i}\right| \leq \sum_{i=1}^{j}\left|l_{j i}\right| \cdot\left|u_{i k}\right| \cdot n \varepsilon=n \varepsilon(|L| \cdot|U|)_{j k}$$

对 $l_{jk}$ 做类似分析， 可得 $$l_{j k}=\left(1+\delta^{\prime \prime}\right)\left(\frac{\left(1+\delta^{\prime}\right)\left(a_{j k}-\sum_{i=1}^{k-1} l_{j i} u_{i k}\left(1+\delta_{i}\right)\right)}{u_{k k}}\right)$$ 其中 $|\delta_i|\leq ( k-1 )\varepsilon$, $|\delta''| \leq \varepsilon$。 同样解出 $a_{jk}$, $$\begin{aligned} a_{j k} &=\frac{1}{\left(1+\delta^{\prime}\right)\left(1+\delta^{\prime \prime}\right)} u_{k k} l_{j k}+\sum_{i=1}^{k-1} l_{j i} u_{i k}\left(1+\delta_{i}\right) \\ &=\sum_{i=1}^{k} l_{j i} u_{i k}+\sum_{i=1}^{k} l_{j i} u_{i k} \delta_{i} & \text { with } 1+\delta_{k} \equiv \frac{1}{\left(1+\delta^{\prime}\right)\left(1+\delta^{\prime \prime}\right)} \\ & \equiv \sum_{i=1}^{k} l_{j i} u_{i k}+E_{j k} \end{aligned}$$ 其中 $|\delta_i| \leq n\varepsilon$， 同样可得：

$$E_{jk} \leq n\varepsilon(|L|\cdot|U|)_{jk}$$

经过上述分析，可以得到如下结论：

$$\begin{align} A =& LU + E \text{ with } |E| \leq n\varepsilon |L|\cdot|U|,\\ \| E \|\leq& n\varepsilon\| |L| \|\cdot\| |U| \|,\\ \| E \|\leq& n\varepsilon\| L \|\cdot\| U \| \text{ for F, inf and one norm, but not two norm}\\ \end{align}$$

再考虑求解 $LUx=b$ 产生的向后误差， 可得(Question 1.11)

$$\begin{aligned} b &=(L+\delta L) \hat{y} \\ &=(L+\delta L)(U+\delta U) \hat{x} \\ &=(L U+L \delta U+\delta L U+\delta L \delta U) \hat{x} \\ &=(A-E+L \delta U+\delta L U+\delta L \delta U) \hat{x} \\ & \equiv(A+\delta A) \hat{x}, \quad \text { where } \quad \delta A=-E+L \delta U+\delta L U+\delta L \delta U \end{aligned}$$

最终可得 $\delta A$ 误差界： $$\begin{align} |\delta A|=&|-E+L \delta U+\delta L U+\delta L \delta U| \\ \leq&|E|+|L \delta U|+|\delta L U|+|\delta L \delta U| \\ \leq&|E|+|L| \cdot|\delta U|+|\delta L| \cdot|U|+|\delta L| \cdot|\delta U| \\ \leq& n \varepsilon|L| \cdot|U|+n \varepsilon|L| \cdot|U|+n \varepsilon|L| \cdot|U|+n^{2} \varepsilon^{2}|L| \cdot|U| \\ \approx& 3 n \varepsilon|L| \cdot|U| \end{align}$$

进而可得

$$\| \delta A \| \leq 3n\varepsilon \| L \|\cdot \| U \|\text{ for F, inf and one norm, but not two norm}$$

如果要高斯消去法向后稳定，需要求 $$\| L \|\cdot \| U \| = O(\varepsilon)\| A \|$$ 进而有 $$\frac{ \| \delta A \| }{ \| A \| } \leq O(\varepsilon).$$

实践经验表明

• $l_{ij} \leq 1$
• 主元增长因子 $g_{PP} = \frac{ \| U \|_{max} }{ \| A \|_{max} }$
• $g_{PP}$ 要比较小或者增长缓慢才能使得高斯消去法是向后稳定的。

估计条件数

• 估计条件数的关键是估计 $\| A^{-1} \|$。
• 直接计算 $A^{-1}$ 的计算复杂度为 $2n^3 > \frac{2}{3} n^3$。
• 避免直接计算 $A^{-1}$ 且代价更低的算法称为条件数估计子， 要具备如下性质
  1. 已经计算出 $L$ 和 $U$ 时，它的计算代价为$O(n^2)$。
  2. 估计子的大小为$C\| A \|$， 则其中 $C$ 要小于 10, 相当于之差了一个 十进制位。

$$\begin{aligned} \| B \|_1 =& \max_{j} \sum_{i=1}^{n} |b_{ij}| \\ \| B \|_1 =& \max_{x\not=0} \frac{ \| Bx \|_1 }{ \| x \|_1 } \end{aligned}$$

• 两个等式代表两种等价的计算方式。
• 但计算效率不等效。

$$f(x) = \max_{\| x \|_1\leq 1}\| Bx \|_1$$

• $\| x \|_1 \leq 1$ 是一个凸集。
• $f(x)$ 是一个凸函数。

上述求极大值的问题可以用梯度上升法解决, 这个方法的关键是计算 $\nabla f$ 。

因为有绝对值的存在，在零点处函数不可导，所以可假设 $$f(x) = \sum_i|\sum_j b_{ij}x_j| \text{ with all } \sum_{j} b_{ij}x_j \not=0.$$

注意，这个假设是几乎总是成立的。

引入 $$\zeta_i = \text{sign}(\sum_j b_{ij}x_j) = \pm 1$$

可得

$$f(x) = \sum_i\sum_j \zeta_i b_{ij}x_j.$$

从而可得 $$\frac{ \partial f }{ \partial x_k } = \sum_i \zeta_i b_{ik}$$

$$\nabla f = \zeta^T B = (B^T\zeta)^T$$

算法 2.5 Hager's condition estimator returns a lower bound $\| w \|_1$ on $\| B \|_1$: $$\begin{align} &\text{choose any } x \text{ such that } \| x \|_1 = 1\\ &\text{repeat }\\ &\quad\text{if }\|z\|_\infty \leq z^Tx \text{ then}\\ &\quad\quad\text{return }\|w\|_1\\ &\quad\text{else}\\ &\quad\quad x = e_j \text{ where } |z_j| = \|z\|_\infty\\ &\quad text{endif}\\ &\text{end repeat}\\ \end{align}$$

• 初值取 $x_i=\frac{1}{n}$
• $z^T=\nabla f$

定理 2.6

1. 如果返回 $\|w\|_1$，则$\|w\|_1 = \|Bx\|_1$ 是 $\|Bx\|_1$ 的局部最大值 。
2. 否则 $\|Be_j\|_1 (\text{at end of loop}) > \|Bx\|(\text{at start})$

证明：

1. 目标是证明 $x$ 是局部极大值点。易知在 $x$ 附近， $f(x)$ 是线性函数，则 有 $$f(y) = f(x) + \nabla f(x)\cdot(y-x) = f(x)+z^T(y-x)$$ 其中 $z^T = \nabla f(x)$。取 $\|y\|_1=1$, $$\begin{aligned} z^T(y-x) =& z^Ty - z^Tx = \sum_i z_i\cdot y_i - z^Tx \leq \sum_i |z_i|\cdot |y_i| - z^Tx \\ \leq & \|z\|_\infty\cdot \|y\|_1 - z^Tx = \|z\|_\infty - z^Tx\leq 0. \end{aligned}$$

2. 如果 $\|z\|_\infty > z^Tx$。设 $|z_j| = \|z\|_\infty$，取 $\tilde x = e_j\cdot\text{sign}(z_j)$， 则 $$\begin{aligned} f(\tilde x) \geq& f(x) + \nabla f\cdot(\tilde x - x) = f(x) + z^T(\tilde x - x)\\ =& f(x) + z^T\tilde x - z^Tx = f(x)+|z_j| - z^Tx > f(x) \end{aligned} \quad (2.12)$$

证毕。$\square$

估计相对条件数 $\kappa_{CR} = \||A^{-1}|\cdot|A|\|_\infty$

$$\frac{\|\delta x\|}{\|x\|} \leq \epsilon\left\|\left|A^{-1}\right| \cdot|A|\right\|. \quad (2.8)$$

或者估计 $\left\|\left|A^{-1}\right| \cdot|r|\right\|$。

$$\|\delta x\|=\left\|A^{-1} r\right\| \leq\left\|\left|A^{-1}\right| \cdot|r|\right\|. \quad (2.9)$$

上面的两个估计问题，可以转化为相同的问题，即估计 $\||A^{-1}|\cdot g\|_\infty$， 其中 $g$ 是一个元素非负的向量。

利用引理 1.7的第 5 个结论， 可知如果矩阵 $X$ 的元素都是非负的，$e = (1, 1, \cdots, 1)^T$ 则 $$\|X\|_\infty = \|Xe\|_\infty.$$ 可推出 $$\left\|\left|A^{-1}\right| \cdot|A|\right\|_{\infty}=\left\|\left|A^{-1}\right| \cdot|A| e\right\|_{\infty}=\left\|\left|A^{-1}\right| \cdot g\right\|_{\infty}, \quad \text { where } g=|A| e.$$

引入对角矩阵$G = \text{diag}(g_1, g_2, \cdots, g_n)$, 则有 $g = Ge$. $$\begin{aligned} \left\|\left|A^{-1}\right| \cdot g\right\|_{\infty} &=\left\|\left|A^{-1}\right| \cdot G e\right\|_{\infty}=\left\|\left|A^{-1}\right| \cdot G\right\|_{\infty}=\left\|A^{-1} G |\right\|_{\infty} \\ &=\left\|A^{-1} G\right\|_{\infty} \end{aligned}$$

最后转化为 $$\left\|(A^TG)^T\right\|_1 = \left\|A^TG\right\|_\infty$$

2.4.4 实用的误差界

$$\text{error} = =\frac{\|\hat{x}-x\|_{\infty}}{\|\hat{x}\|_{\infty}} \leq\left\|A^{-1}\right\|_{\infty} \cdot \frac{\|r\|_{\infty}}{\|\hat{x}\|_{\infty}} \quad (2.13)$$

$$\text{error}=\frac{\|\hat{x}-x\|_{\infty}}{\|\hat{x}\|_{\infty}} \leq \frac{\left\|\left|A^{-1}\right| \cdot|r|\right\|_{\infty}}{\|\hat{x}\|_{\infty}}\quad (2.14)$$

课后作业

1. 第一章 Question 1.10
2. 第一章 Question 1.11