数值分析笔记

这学期选了一个数值分析。说老实话,是实在找不到特别有意思的课,稍微有点凑学分拔将军的意味。不过这课教的知识的确非常有用,所以最好还是写一下笔记,以备以后不时之需。这课用的教材是heath的scientific computing。这笔记基本上就是抄书。

数值分析笔记

观念比原理更重要,原理比方法更重要。

第一章

第一章是对于科学计算中重要概念的介绍,没有提到具体的算法。从中学到了一些基本的考虑问题的抓手。

  • 科学计算的核心就是误差和误差分析。误差大致有两类,一种是“软性”起源的,比如说解析函数需要做taylor展开,在我们需要的范围内留前几项,或者把非线性方程近似成线性来解,或者解稀疏矩阵,等等。另一种是“硬性”起源的——我这么形容的时候,想的是软件和硬件的区别——它指的是计算机进行一般运算时的精度问题。

    比如说,求导。对于一个小量h,$f^{\prime}(x) \approx \frac{f(x+h)-f(x)}{h}$就是不考虑二阶导的结果。考虑了呢?会多一项$h\cdot\frac{|f’’(z)|}{2}\propto h$。这就是“软性”起源的truncation error。另一方面,实际计算$\frac{f(x+h)-f(x)}{h}$的时候,由于计算机的精度最多只能到$\epsilon>0$,这会导致这个计算多了$2\epsilon/h\propto h^{-1}$的误差,这是“硬性”起源的rounding error。

    有一些时候,计算过程中有许多步。我们担心误差会被步步放大。定量描述这种现象的概念是conditioning number $\frac{|\Delta y / y|}{|\Delta x / x|}$。还有forward/backward error的概念,分别指的是因变量的误差和其对应的自变量的误差,我们也可以把conditioning number理解为relative fwd error和relative bkwd error之比:$|(\Delta xf’(x)/f(x))/(\Delta x/x)|$,或$|xf’(x)/f(x)|$。

    计算机喜欢处理连续而不过分陡峭的数学问题。如果不连续,就斥之为ill-posed;连续但过分陡峭的,斥之为ill-conditioned。

  • 理解Rounding error,或计算机精度导致的误差,需要先理解计算机的浮点数机制。浮点数的表示公式如下:

    $x=\pm\left(d_{0}+\frac{d_{1}}{\beta}+\frac{d_{2}}{\beta^{2}}+\cdots+\frac{d_{p-1}}{\beta^{p-1}}\right) \beta^{E}$

    其中,beta是进制,E是$[L, U]$范围内的一个数,描述这个浮点数的数量级;p是精度,所以括号内的部分代表有效数字。一般浮点数都是要normalize的,意思就是说,除非x=0,否则$d_0$不能是0。这不难理解,我们只要考虑科学计数法,并与上面的公式作比较,就能明白。

    Remark: 我居然才意识到科学计数法里面的“科学”居然指的是科学计算。我小时候一直以为是表示大数方便才这么规定,现在才发现,原来是计算上的需求。

    会有一个underflow level $\mathrm{UFL}=\beta^{L}$ 和一个overflow level $\mathrm{OFL}=\beta^{U+1}\left(1-\beta^{-p}\right)$。

    机器处理那些精度过高的数字的时候会给round到一定精度去,这会导致一个误差。比如,如果是四舍五入,那么可能出现的(在有效数字部分的)最大误差为$\epsilon_{\mathrm{mach}}=\frac{1}{2} \beta^{1-p}$。

    浮点数运算会造成一定的误差。一种情况是,数量级不同的数相加,较小的数会被损失掉很多信息。例如,对于上述$\epsilon_{mach}$,有$(1+\epsilon)+\epsilon=1$,而$1+(\epsilon+\epsilon)>1$。一般的,对于加减乘除来说,有$\mathrm{fl}(x$ op $y)=(x$ op $y)(1+\delta)$(fl指float,op指operation)。其中$|\delta| \leq \epsilon_{\mathrm{mach}}$限制了这部分rounding error的程度。这很好理解,并且可以直接用于误差分析。

    另一种情况是,首位数字及数量级相同的数相减,会放大这些数本来自带的误差,因为误差的绝对大小不变的情况下,产生误差的数却小了很多。这叫cancellation error。

第二章

第二章粗浅介绍了一下解线性系统的事情。

  • 引入了矩阵的范数来讨论误差分析的事情,因为矩阵的范数有性质$|\boldsymbol{A} \boldsymbol{x}| \leq|\boldsymbol{A}| \cdot|\boldsymbol{x}|$。例如,对于线性系统Ax=b,如果b出了一些偏差,那么x的误差将会被控制在$\frac{|\Delta \boldsymbol{x}|}{|\boldsymbol{x}|} \leq\left|\boldsymbol{A}^{-1}\right| \cdot|\Delta \boldsymbol{b}| \frac{|\boldsymbol{A}|}{|\boldsymbol{b}|}$的范围,我们说$|\boldsymbol{A}| \cdot\left|\boldsymbol{A}^{-1}\right|$是系统的conditioning number。

  • 引入了残差的概念:对于有偏差的解$\hat{\boldsymbol{x}}$,残差是$r=b-A \hat{x}$,但更合理的应该是scale过的“相对”残差,即$|r| /(|A| \cdot|\hat{x}|)$。同样,$\frac{|\Delta \boldsymbol{x}|}{|\hat{\boldsymbol{x}}|} \leq \operatorname{cond}(\boldsymbol{A}) \frac{|\boldsymbol{r}|}{|\boldsymbol{A}| \cdot|\hat{\boldsymbol{x}}|}$,所以残差与误差之间还隔着一个condition。

  • 高斯消元法和L-U分解,这是复习大一的内容了。回忆一下,我们把一个矩阵乘以许多的下三角矩阵($L^{-1}$),它们各处理一个列,最终把待分解的矩阵变成一个上三角矩阵($U$)。

    在做三角矩阵的时候,我们如果想把矩阵的某一列中,某一元素以上/下的元素都变成0,常规的办法是用这一元素与那些变零的元素来消去,这意味着这一元素(pivot)不能是0。实际计算中,应当避免太小的pivot,因为它的倒数太大,会造成严重的rounding误差。

    这一方法(LU分解)需要的时间复杂度是$\mathcal{O}(n^{3})$,但分解完之后求解LUx=b的时间复杂度则是$\mathcal{O}(n^{2})$,比求逆矩阵快,而且精度高。分解完了之后,对于不同的b都可以解;对于稍微变动(变了一个秩一矩阵)的A也可以用$\mathcal{O}(n^{2})$的Sherman-Morrison公式来处理,仍然可以解得快。

    对于正定对称阵,可以用Cholesky分解,这是说$\boldsymbol{A}=\boldsymbol{L} \boldsymbol{L}^{T}$。这种办法的时间和空间复杂度都约只有普通LU分解的一半。

第三章

第三章介绍了最小二乘法。

  • 在解决$A x=b$问题时,如果A不是方阵而使得问题overdetermined了,就会成为一个最小二乘法问题$A x \cong b$。人们这时希望求解使得残差$r=b-A x$在2-norm之下最小的x。这种情况是很常见的,一个例子就是多项式拟合,因为通常数据点的数目远多于多项式的次数。

  • 可以由$\phi(\boldsymbol{x})=|\boldsymbol{r}|_{2}^{2}=\boldsymbol{r}^{T} \boldsymbol{r}=(\boldsymbol{b}-\boldsymbol{A} \boldsymbol{x})^{T}(\boldsymbol{b}-\boldsymbol{A} \boldsymbol{x})=\boldsymbol{b}^{T} \boldsymbol{b}-2 \boldsymbol{x}^{T} \boldsymbol{A}^{T} \boldsymbol{b}+\boldsymbol{x}^{T} \boldsymbol{A}^{T} \boldsymbol{A} \boldsymbol{x}$求导解出最优的$\boldsymbol{A}^{T} \boldsymbol{A} \boldsymbol{x}=\boldsymbol{A}^{T} \boldsymbol{b}$。另有一种几何的办法,即注意到最优的情形下事实上对应将b投影到Ax张成的空间,$A x=P b$,这时在两边同乘$\boldsymbol{A}^{T}$可得到相同结果。

  • 在condition的时候,按理说可以使用伪逆来计算conditioning number,但又注意到,如果b几乎垂直于Ax张成的平面,那么b的微扰会极大干扰x的求解。为将这一事实纳入考虑,可用$\frac{|\boldsymbol{A} \boldsymbol{x}|_{2}}{|\boldsymbol{b}|_{2}}=\frac{|\boldsymbol{y}|_{2}}{|\boldsymbol{b}|_{2}}=\cos (\theta)$去bound这一因素,经计算可得

    $\frac{|\Delta \boldsymbol{x}|_{2}}{|\boldsymbol{x}|_{2}} \leq\operatorname{cond}_{previous}(\boldsymbol{A}) \frac{1}{\cos (\theta)} \frac{|\Delta \boldsymbol{b}|_{2}}{|\boldsymbol{b}|_{2}}$

    对A的微扰造成的影响复杂些,但精神是一样的。

  • 求解$\boldsymbol{A}^{T} \boldsymbol{A} \boldsymbol{x}=\boldsymbol{A}^{T} \boldsymbol{b}$有多种办法。如A满秩,那么$A^TA$可做Cholesky分解。但数值上还有诸多tricky之处,例如矩阵乘法带来的cancellation等。故人们也去探索别的方法,而非仅从此式求解最小二乘问题。

    • 例如可将问题embed到更大的问题中,使得解线性问题时的pivot更加友好。如,注意到残差垂直于span(A),故可以写下

      $\left[\begin{array}{ll}\boldsymbol{I} & \boldsymbol{A} \\ \boldsymbol{A}^{T} & \boldsymbol{O}\end{array}\right]\left[\begin{array}{l}\boldsymbol{r} \\ \boldsymbol{x}\end{array}\right]=\left[\begin{array}{l}\boldsymbol{b} \\ \boldsymbol{o}\end{array}\right]$

      当然还得properly scale r和x,于是实际上是$\left[\begin{array}{ll}\alpha \boldsymbol{I} & \boldsymbol{A} \\ \boldsymbol{A}^{T} & \boldsymbol{O}\end{array}\right]\left[\begin{array}{c}\boldsymbol{r} / \alpha \\ \boldsymbol{x}\end{array}\right]=\left[\begin{array}{l}\boldsymbol{b} \\ \boldsymbol{o}\end{array}\right]$的形式,实验上一般凭经验取$\alpha=\max _{i, j}\left|a_{i j}\right| / 1000$即可。

    • 另外还有一种办法。若我们的A可以变成上三角矩阵,那么可以将问题划分成$\left[\begin{array}{l}\boldsymbol{R} \\ \boldsymbol{O}\end{array}\right] \boldsymbol{x} \cong\left[\begin{array}{l}c_{1} \\ c_{2}\end{array}\right]$,这个时候实际上可以求得$|\boldsymbol{r}|_{2}^{2}=\left|\boldsymbol{c}_{2}\right|_{2}^{2}$。于是只需要考虑问题的上半部分就可以上一章的办法来解了。为达到这个情况所采用的的手段叫做QR分解。就是把A写成$\boldsymbol{A}=\boldsymbol{Q}\left[\begin{array}{l}\boldsymbol{R} \\ \boldsymbol{O}\end{array}\right]$的形式。

  • QR分解当然可以使用前述LU分解的类似手段。但是我们希望能够保A的2-范数,这个时候应当做一些正交变换。有许多种这样的变换。

    • householder(这是个人名)变换。首先,设想存在一个正交矩阵叫做householder矩阵,它作用在一个向量a上,得到的a’只有前n个元素非零。那么我们就可以用一连串的householder矩阵将A造成一个上三角矩阵。

      实际上,将变换的头尾连起来,若它们的长度相等,那么过此连线v的中点做一垂直平面,则可认为变换前后的向量是以此平面为镜面做的反射变换的结果。故我们算出a的模长,得到其连线v的方向$a-\alpha a’$。

      用$\boldsymbol{P}=\boldsymbol{v}\left(\boldsymbol{v}^{T} \boldsymbol{v}\right)^{-1} \boldsymbol{v}^{T}=\left(\boldsymbol{v} \boldsymbol{v}^{T}\right) /\left(\boldsymbol{v}^{T} \boldsymbol{v}\right)$可以将a投影到v的方向。故用$I-P$可投影到上述垂直平面,而$H=I-2 P=I-2\left(v v^{T}\right) /\left(v^{T} v\right)$则意味着投影到a’。H当然也确实是正交的,这就是householder矩阵。

    • givens(这也是个人名)旋转。非常简单,就是取两个分量,拉出来做一个平面旋转, 然后把其中一个转到0。例如

      $\left[\begin{array}{rrrrr}1 & 0 & 0 & 0 & 0 \\ 0 & c & 0 & s & 0 \\ 0 & 0 & 1 & 0 & 0 \\ 0 & -s & 0 & c & 0 \\ 0 & 0 & 0 & 0 & 1\end{array}\right]\left[\begin{array}{c}a_{1} \\ a_{2} \\ a_{3} \\ a_{4} \\ a_{5}\end{array}\right]=\left[\begin{array}{c}a_{1} \\ \alpha \\ a_{3} \\ 0 \\ a_{5}\end{array}\right]$

    • Gram-Schimidt正交化。这个也是很久以前学过的,把下一个向量投影到已经处理好的部分张成的空间,然后把这个投影砍掉,只留下垂直的部分,然后scale。

    • 如果A不满秩,我们在正交化的每一步应该用permutation矩阵把范数大的列挪到前面来。这叫column pivoting。这样最后就得到$Q^{T} A P=\left[\begin{array}{ll}R & S \\ O & O\end{array}\right]$,其中R是一个上三角的方阵,并且满秩。

  • SVD也是一个经典方法。

  • 不加说明地列举一下各个方法的性能。

    若A为m*n矩阵,那么在$m \approx n$时,householder QR不如最普通的办法($\boldsymbol{A}^{T} \boldsymbol{A} \boldsymbol{x}=\boldsymbol{A}^{T} \boldsymbol{b}$)。$m \gg n$时,普通法更优。但householder法的误差基本永远更小。

    cholesky分解在$\operatorname{cond}(\boldsymbol{A}) \approx 1 / \sqrt{\epsilon_{\mathrm{mach}}}$时不可用,而householder法仅在$\operatorname{cond}(\boldsymbol{A}) \approx 1 / \epsilon_{\mathrm{mach}}$时才不可用。

    SVD是最稳定,同时也最慢的。

第四章

特征值和与之有关的性质刻画了线性算子的几何性质。这一章就来讲一讲特征值问题的数值方法。

  • 一些评论

    • 数学上,需要意识到,由于五次以上的方程没有根式解,大矩阵的特征值也很难去精确求解。不过,我们有收敛很快的近似方法。
    • 我们可以大致划定特征值的范围。Gershgorin定理说明本征值被限制在一系列圆盘里面,其圆心为对角元,而半径则是这一行其他元素的模长之和。
  • 问题的敏感性

    • 取方阵A,其有本征值$\lambda_i$,本征向量组X。对A做微扰,变为$A+E$,有本征值$\mu$,则(可以证明):

      $\left|\mu-\lambda_{k}\right|\le\operatorname{cond}_{2}(\boldsymbol{X})|\boldsymbol{E}|_{2}$。

      因此,本征向量的性质决定conditioning的好坏。本征向量越是倾向于两两垂直,则conditioning就越好。正交矩阵也因此总是well-conditioned的。对于接近不满秩的矩阵,则本征向量接近于线性相关,而conditioning不好。

  • 问题的清洗

    进行预处理,尽量不改变本征值,而将问题转化为可解的形式。一种最普遍的办法是利用相似形$B=T^{-1} A T$,其中T满秩。则:

    $B y=\lambda y \quad \Rightarrow \quad T^{-1} A T y=\lambda y \quad \Rightarrow \quad A T y=\lambda T y$。

    例如,你可以(通过一些相似变换)把矩阵变成General Schur form

    $\boldsymbol{A}=\left[\begin{array}{cccc}\boldsymbol{A}_{11} & \boldsymbol{A}_{12} & \cdots & \boldsymbol{A}_{1 p} \\ & \boldsymbol{A}_{22} & \cdots & \boldsymbol{A}_{2 p} \\ & & \ddots & \vdots \\ & & & \boldsymbol{A}_{p p}\end{array}\right]$

    其中$\lambda(\boldsymbol{A})=\bigcup_{j=1}^{p} \lambda\left(\boldsymbol{A}_{j j}\right)$。

    这基本上是最坏的情况。如果好一些,A是实的,那么该form的对角元要么是实的本征值,要么是2*2的block,对应一组复的本征值。如果还要好一些(但不是同一个意义上的好),就是说,如果A是满秩的,那么对于相似变换$\boldsymbol{B}=\boldsymbol{T}^{-1} \boldsymbol{A T}$,可以让B是diag$(\lambda_i)$,T是本征向量。如果再更好,那么B和T的性质也可以更好。

    对于下面的一些办法,例如QR迭代,应该尽量把矩阵洗干净,弄成尽量tridiagonal(原矩阵对称的情况)或hessenberg(原矩阵不对称的情况)的形式。这样做,可以把每次QR的花销从$\mathcal{O}\left(n^{3}\right)$ 降到 $\mathcal{O}\left(n^{2}\right)$,并且减少QR的次数。总的来说,QR迭代是一个$\mathcal{O}\left(n^{3}\right)$的算法。

  • 计算方法。

    • Power iteration。任取向量$x_{0}=\sum_{j=1}^{n} \alpha_{j} v_{j}$,那么乘以A(并重新归一化)许多次,肯定会迭代到对应最大本征值的那个本征向量。当然,有可能会失败,比如$\alpha_0$恰巧为0,或者最大的本征值是简并的,等等。
      • 换成A的逆,这么做就可以知道最小的本征值。
      • 平移一下,弄成$\boldsymbol{A}-\sigma \boldsymbol{I}$(的逆),就可以知道离$\sigma$最近的本征值。(这就是inverse iteration强的地方)恰当的平移可以在任何情况下改进最大和次大的本征值的比例。
  • 瑞利商。取一个x,当成是近似的本征向量,可以像求解最小二乘问题那样求解本征值:$\boldsymbol{x} \lambda \cong \boldsymbol{A} \boldsymbol{x}$。结果是$\lambda=\frac{x^{T} A x}{x^{T} x}$,这被称为瑞利商。在用A的逆去迭代求解时,你可以先得到一个迭代的结果x,算出瑞利商$\sigma$,然后平移。

  • 算出一个本征值之后,就像算出多项式的根那样,你可以除以(x-a),然后去求解别的根。

    • 一种做法是,用Householder变换对第一个本征向量得出$H x_{1}=\alpha e_{1}$,然后以$H A H^{-1}=\left[\begin{array}{cc}\lambda_{1} & \boldsymbol{b}^{T} \\ \mathbf{0} & \boldsymbol{B}\end{array}\right]$得到B,B的本征值就是接下来要求的了。但它的本征向量和A不一样。求出B的本征向量$y_2$后,求出$\alpha=\frac{\boldsymbol{b}^{T} \boldsymbol{y}_{2}}{\lambda_{2}-\lambda_{1}}$,再求出A的本征向量$x_{2}=H^{-1}\left[\begin{array}{l}\alpha \\ y_{2}\end{array}\right]$。
    • 或者,取$\boldsymbol{u}_{1}^{T} \boldsymbol{x}_{1}=\lambda_{1}$,然后$\boldsymbol{A}-\boldsymbol{x}_{1} \boldsymbol{u}_{1}^{T}$。我觉得这样更自然。
  • 还有一些迭代法。比如,从$n \times p$的$\boldsymbol{X}_{0}$开始,一次迭代p个向量,$\boldsymbol{X}_{k}=\boldsymbol{A} \boldsymbol{X}_{k-1}$,你将会得到最大的p个本征值所对应的那些本征向量张成的空间。但是,你会总是需要归一化你的各个向量。一个简单的办法是每轮都做QR分解,然后从单位矩阵$X_{0}=I$开始,用

    $\begin{aligned} \hat{\boldsymbol{Q}}_{k} \boldsymbol{R}_{k} &=\boldsymbol{X}_{k-1} \\ \boldsymbol{X}_{k} &=\boldsymbol{A} \hat{\boldsymbol{Q}}_{k} \end{aligned}$

    来取代之前的过程。我们可以把这个认为是每轮更新Q,使得$\boldsymbol{A}_{k}=\hat{\boldsymbol{Q}}_{k}^{H} \boldsymbol{A} \hat{\boldsymbol{Q}}_{k}$是一个三角之类容易算本征值的形式($H\equiv\dagger$)。同样,也可以在这个里面引入平移来加速整个计算。从$\boldsymbol{A}_{0}=\boldsymbol{A}$开始,用

    $\boldsymbol{Q}_{k} \boldsymbol{R}_{k}=\boldsymbol{A}_{k-1}-\sigma_{k} \boldsymbol{I}$
    $\boldsymbol{A}_{k}=\boldsymbol{R}_{k} \boldsymbol{Q}_{k}+\sigma_{k} \boldsymbol{I}$

    来迭代。这么做能行,是因为在迭代$\boldsymbol{A}^{k}=\hat{\boldsymbol{Q}}_{k} \hat{\boldsymbol{R}}_{k}$的同时,你实际上也迭代了$A^\dagger$的逆:$\left(\boldsymbol{A}^{-H}\right)^{k}=\hat{\boldsymbol{Q}}_{k} \hat{\boldsymbol{R}}_{k}^{-H}$(因为Q是厄密的)。实用上说,$\boldsymbol{A}_{k-1}$的右下角取一个2*2矩阵,它的本征值可以作为对平移距离$\sigma_k$(或对最大本征值)的估计。或者更粗略,直接取右下角的元素$a_{n n}^{(k-1)}$。

  • Krylov子空间方法。首先,将Power iteration的各步作为列向量排成一个矩阵:

    $\boldsymbol{K}_{k}=\left[\begin{array}{llllll}\boldsymbol{x}_{0} & \boldsymbol{x}_{1} & \cdots & \boldsymbol{x}_{k-1}\end{array}\right]=\left[\begin{array}{lllll}\boldsymbol{x}_{0} & \boldsymbol{A} \boldsymbol{x}_{0} & \cdots & \boldsymbol{A}^{k-1} \boldsymbol{x}_{0}\end{array}\right]$。

    该序列的元素张成Krylov子空间$\mathcal{K}_{k}=\operatorname{span}\left(\boldsymbol{K}_{k}\right)$。此时,我们发现

    $\begin{array}{llllll}\boldsymbol{A} \boldsymbol{K}_{n} & = & {\left[\begin{array}{llll}\boldsymbol{A} x_{0} & \cdots & \boldsymbol{A} \boldsymbol{x}_{n-2} & \boldsymbol{A} \boldsymbol{x}_{n-1}\end{array}\right]}\end{array}=\boldsymbol{K}_{n}\left[\boldsymbol{e}_{2} \quad \cdots \quad \boldsymbol{e}_{n} \quad \boldsymbol{a}\right] \equiv \boldsymbol{K}_{n} \boldsymbol{C}_{n}$

    而$\boldsymbol{a}=\boldsymbol{K}_{n}^{-1} \boldsymbol{x}_{n}$。注意到这里C是一个Hessenberg矩阵。这样,实际上就是$K_{n}^{-1} A K_{n}=C_{n}$。进一步,由于Kn不是一组好的基,我们应当做QR分解,$Q_{n} R_{n}=K_{n}$,这样仍然是$\boldsymbol{Q}_{n}^{H} \boldsymbol{A} \boldsymbol{Q}_{n}=\boldsymbol{R}_{n} \boldsymbol{C}_{n} \boldsymbol{R}_{n}^{-1} \equiv \boldsymbol{H}$,其中H是个hessenberg矩阵。

    我们使用Arnoldi迭代来计算Q的列向量q,并用q得到H,用H得到A的本征值和本征向量。H的本征值叫做Ritz values,它们近似等于A的本征值。用Qk乘以H的本征值(Ritz向量),就得到近似的A的本征值。计算q们所使用的办法就是一边乘以A,一边Gram-Schimidt正交化。这种办法每次需要动用此前算出的所有q,因此时间和内存上都比较贵。

    如果待分解的矩阵是对称或厄密的,那么计算起来会容易很多。我们使用Lanczos迭代(这就是个名字)来计算。为什么?因为这个时候H是tridiagonal的,我们并不需要每次正交化q的时候都让它和之前的所有q都垂直,而只需要让它和它之前的两个q垂直,剩下的就自动垂直。所以我们大量减少了计算量和内存。

    • 雅克比方法。这是一种计算对阵矩阵本征值的办法。它通过平面旋转来一次消去待求解矩阵的对称的两个元素。但是这种方法比较慢,而且也很难推广到非对称的矩阵。

第五章

这一章的主题是非线性问题。这样的问题可以归结为求一个非线性函数或函数组的零点。

  • 解的存在性和唯一性。一般来说,这是比较复杂的。一般我们在一个区间或区域内求解,这叫做bracket,而在bracket上使用中值定理、反函数定理可以得到一些关于解的数学性质(比如存在性)。但是总的来说,似乎没有什么好的办法。
  • 敏感性。absolute condition number,可以想见,是雅可比矩阵决定的。在某根$x*$处的此数值为$\left|\boldsymbol{J}_{f}^{-1}\left(\boldsymbol{x}^*\right)\right|$。
  • 迭代方法。
    • 人们用指数记录迭代方法收敛的速度,例如,线性收敛的convergence rate是1。
    • 二分法。它当然是线性收敛的。还有不动点法,一般来说,它的收敛速度是和泰勒展开有关的,如果在不动点处$g^{\prime}\left(x^{*}\right)=0$,那么就肯定至少是二次以上的快速收敛。
    • 牛顿法。这就是说,$x_{k+1}=x_{k}-\frac{f\left(x_{k}\right)}{f^{\prime}\left(x_{k}\right)}$,或$x_{k+1}=x_{k}-J\left(x_{k}\right)^{-1} f\left(x_{k}\right)$。这相当于解$g(x)=x-f(x) / f^{\prime}(x)$的不动点,可以知道,一般情况下,牛顿法至少是二次的收敛(但要离根够近,否则不收敛)。计算雅克比矩阵的成本是$n^{2}$,n是维数,而解线性系统的成本(回忆两章前)是$\mathcal{O}\left(n^{3}\right)$。
    • 割线法。牛顿法有个问题,就是每次都要在下一个点重新算一遍导数的值。割线法就是用$f^{\prime}\left(x_{k}\right) \approx \frac{f\left(x_{k}\right)-f\left(x_{k-1}\right)}{x_{k}-x_{k-1}}$来算这个导数。在二维以上,如果最初的Jacobian是$B_{0}$,那么每一轮的jacobian是$\boldsymbol{B}_{k+1}=\boldsymbol{B}_{k}+\left(\left(\boldsymbol{y}_{k}-\boldsymbol{B}_{k} \boldsymbol{s}_{k}\right) \boldsymbol{s}_{k}^{T}\right) /\left(\boldsymbol{s}_{k}^{T} \boldsymbol{s}_{k}\right)$(这个时候叫Broyden方法)。它的conv rate是1.618。
    • 反向内插。实际上,可以认为割线法里面用了一次函数去内插下一轮的解,而我们也可以用更高次的函数。例如,用二次函数,用前三轮去(lagrange法)逼近一个二次函数,找它的零点。它的conv rate是1.839。
    • 一些安全方法。安全是说,其他的一些方法虽然更快,但如果起点离实际值太远,则可能不收敛,而安全的方法不会有这个情况。特别是,如果某一轮使用快速方法的结果跑到bracket外面去了,那么这时应该做一轮安全方法。这些办法的核心基本上就是不断搜索适当的步长和方向,使得这一步迭代真的能提高精度。