非线性方程(组)：高维方程解法

　　非线性方程的高维情形和一维情形既有相似处也有差异。首当其中的区别即在高维情形中不再存在介值定理，从而使得二分法不再可推广到高维。不过，仍然有许多方法可以推广。

1. 不动点迭代(高维)

　　寻找方程 $\boldsymbol{x}=\boldsymbol{g}(\boldsymbol{x})$ 的解。其中 $\boldsymbol{g}:\mathbb{R}^n\rightarrow \mathbb{R}^n$ 为一个n维空间上的变换。迭代形式：$\boldsymbol{x}_{k+1}=\boldsymbol{g}(\boldsymbol{x}_k)$ 。不动点在高维空间的表现为点，经过迭代后位置不变。

项目	一维不动点迭代	高维不动点迭代	高维解释
迭代形式	$x_{k+1}=g(x_k)$	$\boldsymbol{x}_{k+1}=\boldsymbol{g}(\boldsymbol{x}_k)$
收敛条件	$x=g(x), |g'(x)|<1$	$\boldsymbol{g}(\boldsymbol{x}), \rho(\boldsymbol{G}(\boldsymbol{x}^*))<1$	$\delta \boldsymbol{x}_{k+1}\approx \boldsymbol{G}(\boldsymbol{x}^*)\delta \boldsymbol{x}_k$ $||\delta \boldsymbol{x}_{k+1}|| \leq \rho(\boldsymbol{G}(\boldsymbol{x}^*)) \boldsymbol{x}_k$
收敛速度	线性收敛（当 $g'(x)\neq 0$ ); 至少二次收敛（当 $g'(x)=0$ );	线性收敛（当 $\rho (\boldsymbol{G}(\boldsymbol{x}^*))\neq 0$ ); 至少二次收敛（当 $\rho (\boldsymbol{G}(\boldsymbol{x}^*))=0$ );
迭代复杂度	一次标量函数求值	一次n维矢量函数求值，等价于n个标量函数求值

 

　　$\boldsymbol{G}(\boldsymbol{x})$ 是函数变换 $\boldsymbol{g}$ 在x处的Jacobi矩阵，$\rho$ 表示矩阵的谱半径，即矩阵特征值的最大模长（对于实数即绝对值）。正如 非线性方程(组)：计算基本理论 中提到的，Jacobi矩阵是高维情形下导数的等价形式，更确切地说是线性主部的系数；线性代数的知识可知，特征值的几何意义是矩阵对于特征值对应特征向量方向的向量的伸缩比，如特征值为2的特征方向的向量被矩阵左乘以后方向不变、长度翻倍。综上所述，Jacobi矩阵的谱半径可作为微分系数的大小在高维下的等价。　　

 

2. 牛顿法(高维)

　　在一维情况下，牛顿法用零点附近的微分公式，利用上一个点的值和微商值估算下一个零点的估计值。

　　在高维情况下，同样使用零点附近的微分公式：$\boldsymbol{f}(\boldsymbol{x+\Delta x})=\boldsymbol{f}(\boldsymbol{x})+\boldsymbol{J}_f(\boldsymbol{x})\boldsymbol{\Delta x}$ ，带入 $\boldsymbol{f}(\boldsymbol{x+\Delta x})=\boldsymbol{0}$ ，则得到方程：$\boldsymbol{f}(\boldsymbol{x})+\boldsymbol{J}_f(\boldsymbol{x})\boldsymbol{\Delta x}$ 。若Jacobi矩阵非奇异（对于合理的问题，在零点附近的矩阵也理应是非奇异的，除非是重根）则方程的解可以写作：$\boldsymbol{\Delta x}=-\boldsymbol{J}_f(\boldsymbol{x})^{-1}\boldsymbol{f}(\boldsymbol{x})$ 。

　　此即高维牛顿方法的迭代方程：$\boldsymbol{x}_{k+1}=\boldsymbol{x}_k-\boldsymbol{J}_f(\boldsymbol{x}_k)^{-1}\boldsymbol{f}(\boldsymbol{x}_k)$ 。注意到这里的 $\boldsymbol{f}$ 是一个向量，由Jacobi矩阵左乘仍为一个向量。数值方法中，不用显示地生成Jacobi矩阵的逆，只需要用Jacobi矩阵作为系数矩阵，$\boldsymbol{f}$ 作为右端项，用矩阵分解求解一个线性方程组即可，这样复杂度明显更低。

项目	一维牛顿法	高维牛顿法	解释
迭代形式	$x_{k+1}=x_k-\frac{f(x_k)}{f'(x_k)}$	$\boldsymbol{x}_{k+1}=\boldsymbol{x}_k-\boldsymbol{J}(\boldsymbol{x}_k)^{-1}\boldsymbol{f}(\boldsymbol{x}_k)$ 在计算中一般求出 $\boldsymbol{J}(\boldsymbol{x}_k)\boldsymbol{s}_k=-\boldsymbol{f}(\boldsymbol{x}_k)$ 的解，然后 $\boldsymbol{x}_{k+1}=\boldsymbol{x}_k+\boldsymbol{s}_k$
收敛速度	对应不动点迭代： $g'(x^*)=\frac{f(x^*)f''(x^*)}{(f'(x^*))^2}=0$ 二次收敛	对应不动点迭代： $\boldsymbol{G}(\boldsymbol{x^*})=\boldsymbol{I}-\boldsymbol{J}(\boldsymbol{x}^*)^{-1}\boldsymbol{J}(\boldsymbol{x}^*) - \sum\limits_{i=1}^n f_i(\boldsymbol{x}^*)\boldsymbol{H}_i(\boldsymbol{x}^*)=\boldsymbol{O}$  二次收敛	$\boldsymbol{I}-\boldsymbol{J}(\boldsymbol{x}^*)^{-1}-\sum\limits_{i=1}^nf_i(\boldsymbol{x}^*)\boldsymbol{H}_i=\boldsymbol{I}-\boldsymbol{I}-\boldsymbol{O}（f的各个分量均为零）=\boldsymbol{O} $
迭代步复杂度	两次函数求值（导数+原函数），固定四则运算两次； 另需要求导函数；	n维向量函数 $\boldsymbol{f}$ 求值 + Jacobi矩阵 $\boldsymbol{J} 求值，即攻击 $n^2+n$ 次标量函数求值； 求 $\boldsymbol{J}(\boldsymbol{x}_k)\boldsymbol{s}_k=\boldsymbol{f}(\boldsymbol{x}_k)$ 的解系解n维线性方程组，$T(n)\approx n^3/3, O(n)=n^3$； 另需要求雅可比矩阵关于x的函数。

  

3. 布洛伊登(Broyden)方法（高维拟牛顿法和割线更新法的一种实现）

　　在一维中，割线法通过前两个点的坐标做出过前两个点的割线，将割线与横轴的交点坐标作为下一个零点估计值，实质上是利用线性插值的方法代替牛顿法中的求导给出了切线的估计。割线法的好处是不用计算导函数。计算导函数所要求的符号求导运算可以说时间开销等价于大量的浮点运算，且随着函数的形式变复杂而越来越“昂贵”。对于非符号表达式的函数，计算导函数还需要插值，同样非常昂贵。

　　牛顿法的这一缺陷在高维情况下更为严重：求导数的运算变成了求雅可比矩阵关于x的函数，原本只有一个标量函数符号求导变成了 $n^2$ 个符号求导，更不用说非符号表达式的函数求导函数。除此之外，每一步迭代中解线性方程组要求的 $n^3/3$ 的复杂度更产生了迭代复杂度的最高次项。因此高维的牛顿法尤为昂贵。

　　同时解决这两个问题的一种方法称为割线更新(Secant Updating)法。类似于一维的割线法，它的具体思路包括利用连续几次的取点及其函数值生成Jacobi矩阵的估计值，而不需要先进行符号求导而后每迭代一次再耗费 $n^2$ 次标量函数求值生成。同时，所谓的“更新”是说每一次迭代只需要将之前Jacobi矩阵的分解(factorization)稍作更新改动即可，不用全盘重来，因此就可以每次省下 $n^3$ 量级的运算成本。由于基本框架继承牛顿法，它也称为拟牛顿法(Quasi-Newton Method)。

　　余下的问题即，如何确定这样的Jacobi矩阵的近似矩阵，使得它既可以作为一个估计值，又可以使得每次迭代矩阵分解的变化尽可能少（并且易得）呢？选取近似矩阵是有自由度的，而实现割线更新方法的一种算法就是布罗伊登(Broyden)方法。对于以上两个条件，Broyden方法要求Jacobi的近似矩阵以两个形式满足：

　　可以作为估计值 $\Leftrightarrow $ $\boldsymbol{B}_k$ 满足割线方程 $\boldsymbol{B}_{k+1}(\boldsymbol{x}_{k+1}-\boldsymbol{x}_k)=\boldsymbol{f}(\boldsymbol{x}_{k+1})-\boldsymbol{f}(\boldsymbol{x}_k)$

新矩阵分解可以用较少的计算获得 $\Leftrightarrow$ 相邻的 $\boldsymbol{B}_k$ 满足秩一变换 $\boldsymbol{B}_{k+1}=\boldsymbol{B}_k+\boldsymbol{u}\boldsymbol{v}^T$

　　于是最终选择的迭代式：$$\boldsymbol{B}_{k+1}=\boldsymbol{B}_k+\frac{\boldsymbol{f}(\boldsymbol{x}_{k+1})\boldsymbol{s}_k^T}{\boldsymbol{s}_k^T\boldsymbol{s}_k}$$　　其中，$\boldsymbol{s}_k$ 为线性方程组 $\boldsymbol{B}_k\boldsymbol{s}_k=-\boldsymbol{f}(\boldsymbol{x}_k)$ 的解，而 $\boldsymbol{x}_{k+1}=\boldsymbol{x}_k+\boldsymbol{s}_k$ 。