如何保证矩阵计算的精确度：backward stable & well conditioned

大纲：介绍范数；介绍条件数;介绍条件数的几何意义；介绍有条件数得出的误差估计；介绍backward stable；介绍backward stable + well conditioned 可以得出求解结果精确的结论；介绍如何判断是否backward stable。
介绍实际操作中该怎么做：1）判断是否well conditioned，如果不是该怎么调整；2）通过residual判断是否backward stable。这也就是为什么迭代方法可以以residual小于某一阈值作为结束条件的理论基础

1. 范数

1.1 向量的范数

1.1.1 向量范数定义

向量$x\in R^n$的范数（vector norm）是一个$R^n\to R$的函数，记为$||x||$，且满足以下三个条件：
- $||x||>0\quad if\quad x\ne 0, \quad and \quad ||0||\ne 0$
- $||\alpha x||=|\alpha|||x||$
- $||x+y||\le ||x|| + ||y||$

1.1.2 $p$范数

对于任何实数$p\ge 1$，定义$p$-范数为：

$$||x||_p = \Big ( \sum_{i=1}^n{|x_i|^p}\Big )^{1/p}$$

易证，所有$p$-范数都是满足上述向量范数的三个性质。
常用的p范数有：

• 0-范数：非零项的个数。
• 1-范数：各项的绝对值之和。
• 2-范数：各项平方和再开根，又称Euclidean-范数，就是我们常用的向量的模。
• $\infty -范数$：绝对值最大的值。

  1.2 矩阵范数

  1.2.1 矩阵范数定义

  矩阵$A\in R^{n\times n}$的范数（matrix norm）是一个$R^{n\times n}\to R$的函数，记为$||A||$，且满足以下四个条件：
  - $||A||\quad>\quad 0\quad if \quad A\ne 0$
  - $||\alpha A||\quad=\quad |\alpha|||A||$
  - $||A+B||\quad \le \quad ||A||+||B||$
  - $||AB||\quad \le \quad ||A||||B||$

  相较于向量范数，矩阵范数只是多了第四个性质。注意：根据上面的定义，只有方正才有范数。

  Frobenius 范数定义为各项平方和再开根：
  

  $$||A||_F=\Big( \sum_{i=1}^n{\sum_{j=1}^n{|a_{ij}|^2}}\Big)^{1/2}$$

  这是从直观上我们最想定义的范数，注意与下面要介绍的诱导2-范数区分。

  1.2.2 矩阵的诱导范数

  给定一个向量范数$||\cdot||_v$（例如1-范数），矩阵的诱导范数（induced matrix norm）定义如下：
  

  $$||A||_M=max_{x\ne 0}\frac{||Ax||_v}{||x||_v}$$

  如果把矩阵$A$看成$R^n \to R^n$的线性变换，那么$||Ax||_v/||x||_v$就是从$x$变换到$Ax$时的放大率，因此，矩阵的诱导范数就是最大的放大率。

  可证明，矩阵的诱导范数满足上面定义的矩阵范数的四个性质。

  常用的矩阵范数通常由向量的$p$-范数诱导，称为矩阵 $p$-范数：
  

  $$||A||_p=max_{x\ne 0}\frac{||Ax||_p}{||x||_p}$$
  • 矩阵 2-范数又称谱范数（spectral norm），具有重要的理论意义（以后会介绍），但计算十分耗时。注意矩阵 2-范数并不是上面介绍的Frobenius 范数。
  • 矩阵 1-范数为“列绝对值和”的最大值，又称“列和”-范数（column-sum norm）：
    $$||A||_1 = max_{1\le j \le n}\sum_{i=1}^{n}|a_{ij}|$$
  • 矩阵 $\infty$-范数为“行绝对值和”的最大值，又称“行和”-范数（row-sum norm）：
    $$||A||_{\infty} = max_{1\le i \le n}\sum_{i=1}^{n}|a_{ij}|$$

  2. 条件数

  下文的介绍都基于以下线性系统：$Ax=b$，其中$A$非奇异，$b$非零，$x$是该方程的精确解。

  2.1 条件数定义

  矩阵$A$的条件数（condition number），$k(A)$定义为：
  

  $$k(A)=||A||\cdot||A^{-1}||$$
  由上式可知：
  - 对于同一个矩阵，使用不同的矩阵范数，得到的条件数也是不同的。
  - $k(A)\ge 1$，因此最好的条件是就是1。

  证明：由1.2.1矩阵范数定义可知，$||AB||\ \le ||A||||B||$，那么，$k(A)= ||A||||A^{-1}|| \ge ||AA^{-1}|| = 1$

  2.2 条件数与矩阵计算精度

  2.2.1 $b$有波动对结果的影响

  假设由于某些原因，$b$有些许误差$\delta b$. 令，$\tilde x=x+\delta x$是$A\tilde x=b+\delta b$的解 ，则（我们不加推导的给出结论），
  

  $$\frac{||\delta x||}{||x||}\le k(A)\frac{||\delta b||}{||b||}$$

  2.2.2 $A$有波动对结果的影响

  假设由于某些原因，矩阵$A$有些许误差$\delta A$，令，$\tilde x=x+\delta x$是$(A+\delta A)\tilde x=b$的解，则（我们不加推导的给出结论），
  

  $$\frac{||\delta x||}{||\tilde x||}\le k(A)\frac{||\delta A||}{||A||}$$

  2.2.3 总结

  只要系数矩阵的条件数$k(A)$足够小，误差$\delta x$、$\delta A$足够小，那么，求解结果就是精确的。

  那么我们在求解线性系统时就要做两个工作：

  • 确保$k(A)$足够小：如果不是，那么应从物理模型分析如何调整。
  • 判断$\delta x$、$\delta A$是否足够小：使用后验方法检验。

  2.1 如何估计条件数

  • 条件数等于最大投影和最小投影之比
  • 条件数大于任意两列向量的范数之比
  • $k_2(A)=k_2(A^T)$

  3. backward stable：将算法的误差转化为操作数的波动

  3.1 Backward Error Analysis

  • Forward Error Analysis
    与 Backward Error Analysis相对应的是Forward Error Analysis，后者试图通过对参与运算的操作数（例如，矩阵A和向量b）以及求解算法（例如，LU分解）进行分析，从而在计算前预测最终的计算结果的误差，这通常是很困难的。
  • Backward Error Analysis
    聪明的科学家最终发明了backward error analysis，所谓backward，是指，我已经计算出结果了，然后我想知道我的误差大概有多大。方法是：将未知的最终误差“推给”操作数，转换成操作数误差（在求解线性系统时，即，转换成$\delta A$或$\delta b$）。如果我们能证明，操作数误差是“很小的”，那么就称算法是backward stable的。

  3.2 Backward Stable + Well Conditioned = Accurate

  • backward stable 是用来描述算法的，同一个问题（$Ax=b$）我们可以用不同的算法求解（LU分解，Cholesky分解等）
  • sensitivity analysis是用来分析问题本身的，如果对于一个线性系统$Ax=b$，$k(A)$过大（即，ill-conditioned），那么由于计算机浮点误差和舍入误差（round-off error）的存在，我们不能指望任何一个算法能够精确求解该问题。
  • 如果一个算法是backward stable ，且问题是well conditioned，那么我们可以断定求解是精确的。

  3.3 Small Residual Implies Backward Stability

  那么如何验证算法是否backward stable呢？这里有个后验方法（posterior method）：检测余项（residual）。如果余项较小，那么算法就是backward stable. 证明如下：

  假设求解$Ax=b$，我们得到结果$\tilde x$，计算余项$\tilde r = b-A\tilde x$，易知，$\tilde x$是方程$A\tilde x=b+\delta b,\quad \delta b=\tilde r$的精确解（将$\delta b 带入即可证明$）。那么只要$\tilde r$足够小，那么$||\delta b||/||b||$也就足够小，即，算法backward stable. （当然，也可以将误差“推给”$\delta A$，这里不再证明。）

  4. 实际中如何操作

  问题的根源：计算机的浮点付出和舍入误差

  上面已经提到，要保证求解的精确一要保证系数矩阵的条件数足够小；二要保证求解结果的residual要小。

  4.1 保证小的条件数

  可以使用matlab或其它库直接计算系数矩阵 A A <script type="math/tex" id="MathJax-Element-2470">A</script>的条件数，如果过大，可用尝试如下方法：

  4.1.1 保证列（行）向量范数基本相同

  如果矩阵列范数相差太多，那么应该考虑对某些列进行缩放。

  • If a system is ill conditioned because its rows or columns are badly out of scale, one must refer back to the underlying physical problem in order to determine whether the ill conditioning is inherent in the problem or simply a consequence of poor choices of measurement units.
  • Any decision about whether to rescale or not, and how to rescale, should be guided by the underlying physical problem.

  4.1.2 给对角线元素加正值

  在实践中我发现这是可行的，对于有些问题，还可以证明这么做的数学依据。

  4.1.x ……

  4.2 检测residual是否足够小

  这也是为什么迭代法求解时，可以使用residual作为结束条件的依据。

  5. 总结

  算了，写得跟屎一样，还总结个毛。