非线性方程组求解器算法

非线性方程组求解器算法

总结一下电路仿真器中的非线性方程组求解器算法，即 Newton-Krylov 方法

---

算法概述

1. 外层循环：非线性泰勒展开（牛顿法）

   • 动作：面对复杂的非线性方程 $F(x)=0$，牛顿法确实是在当前点进行一阶泰勒展开。

   • 结果：它丢掉了高阶项，将原本"弯曲"的非线性问题变成了"笔直"的线性方程 $J\cdot \Delta x=-F(x)$

2. 内层计算：线性方程组求解（Krylov 迭代）

   • 动作：现在我们手里有一个巨大的线性方程组。如果电路很大，我们不用"暴力"的高斯消元，而是用 Krylov 子空间迭代（如 GMRES）。

   • 结果：它不存完整的 $J$ 矩阵，而是通过不断的"向量乘法"试探，快速求出 $\Delta x$

3. 约束导向：弧长法加持（弧长约束）

   • 动作：如果此时电路处于极端非线性区（比如增益压缩的拐点），直接按牛顿法的 $\Delta x$ 跳过去可能会"脱轨"。此时弧长法介入，它在方程组里额外加了一个距离约束（源步进法）$(\Delta x)^2+(\Delta(\lambda b))^2=(\Delta s)^2$
   • 结果：它强迫牛顿法的每一次迭代都必须落在以旧解为圆心、弧长为半径的"圆弧"上。

4. 最终归位：拉回解曲线

   • 动作：通过 Krylov 算出的 $\Delta x$ 虽然是在圆弧上的，但不一定在真实的物理解曲线上。
   • 循环：牛顿法继续迭代，利用残差 $F(x)$​ 不断修正，最终将这个点垂直（或按约束方向）拉回到物理曲线上，找到真正的平衡解。

如果不考虑内层线性方程组是如何求解的，将求解一个大型 $m\times m$ 矩阵方程组看成一个已经实现了的工具黑箱，那么牛顿法已经解决了非线性电路求解的问题，前面的博文电路仿真器是如何工作的已经介绍地足够清晰了。

所以这里作为完整求解器算法的补充，主要讲如何求解这个大型矩阵方程组。如何实现在有限的内存下，更快求解出这个成千上万未知量的矩阵方程组。

Krylov 子空间

首先要引进一个概念，Krylov 子空间。为了理解这个子空间，我们可以从最简单的人类认知来讲解。拥有成千上万个未知量的方程组好比是一个解谜游戏，谜底是什么我们完全不知道，但我们可以不断猜测谜底的特征来逐渐接近它。完整的精确的解，也就是谜底，本需要与未知量相同数量的维度去描述。但描述越精确，需要的代价太大，就好比我们计算的内存和时间不够。所以有没有可能，我们将这些特征从粗略到精细排个序，先粗略猜答案的大致范围。例如我们假定谜底是个非生物，带入方程组后发现错了，因此结果告诉我们谜底是生物，这样就能确定解的大致范围。接着我们假定谜底是植物，带入后发现错了，因此缩小范围判断谜底是动物。如此反复猜测迭代，最终就能确定谜底。

同样的道理，精确的解在一个非常高维的空间中，想直接求解非常困难。于是我们想让这个解先降维下来，求解它在低维空间中的近似解，看看这个近似解是否足够精确，也就是收敛。如果在某个低维空间中，这个解已经足够精确收敛了，那我们也没必要去浪费资源硬将它的高维精确解完整求出。求解过程中的用到的低维空间，就是我们这里的 Krylov 子空间。

传统的线性方程组可以表示为

$$Ax=b$$

我们想求解 $x$，直接法例如高斯消元或者 LU 本质都在间接求 $A^{-1}$，但对于一个拥有成千上万维度的矩阵，求逆太困难了。于是我们考虑将 $A^{-1}$ 表示为 $\sum^{m-1}_{k=0}c_kA^k$，也即是

$$x=A^{-1}b\quad\Rightarrow\quad x^{(m)}=\sum^{m-1}_{k=0}c_kA^kb$$

看到这里是不是想起了泰勒级数？没错，这种思想就是矩阵形式的"泰勒级数"展开，理论上只要我展开的项足够多，近似解 $x^{(m)}$ 就足够精确！

因此，高维矩阵求解问题便转化为了求解有限低维 $m$ 级次下的展开系数。Krylov 子空间可以表示为：

$$K_m(A,b)=\text{span}\{b,Ab,A^2b,...,A^{m-1}b\}$$

假如能获取这个子空间的所有正交基 $\{v_1,v_2,...,v_m\}$，和精确解在这个子空间中各个基上的投影，即子空间坐标 $\{y_1,y_2,...,y_m\}$，我们就能得到近似解 $x^{(m)}$ 的表达式

$$x^{(m)}=y_1v_1+y_2v_2+\cdot\cdot\cdot+y_mv_m=V_m y$$

我们的目标明确了，需要求 Krylov 子空间的正交基 $\{v_1,v_2,...,v_m\}$ 和子空间坐标 $\{y_1,y_2,...,y_m\}$。

GMRES 算法

GMRES 算法大家应该在各种仿真软件（例如 COMSOL）中见过，它的全称是 Generalized Minimal Residual method，即广义最小残量方法，属于 Krylov 子空间方法。这里就以它作为例子来说明近似解是如何一步一步求出来的。

我们以解这个方程组为例：

$$\begin{cases} x_1+x_2+x_3=3\\ x_1-x_2+x_3=1\\ x_1+x_2-x_3=1 \end{cases}$$

写成矩阵形式：

$$Ax=b$$

其中

$$A= \begin{bmatrix} 1&1&1\\ 1&-1&1\\ 1&1&-1 \end{bmatrix}, \qquad b= \begin{bmatrix} 3\\1\\1 \end{bmatrix}$$

0 阶残差初始化

按 Krylov 法，先给一个初值。取最常见的

$$x^{(0)}=0$$

则初始残差

$$r_0=b-Ax^{(0)}=b= \begin{bmatrix} 3\\1\\1 \end{bmatrix}$$

其范数为

$$\beta=\|r_0\|=\sqrt{3^2+1^2+1^2}=\sqrt{11}$$

于是 Arnoldi 的第一个标准正交基向量 $v_1$ 为

$$v_1=\frac{r_0}{\|r_0\|} = \frac1{\sqrt{11}} \begin{bmatrix} 3\\1\\1 \end{bmatrix}$$

这里相当于得到了 Krylov 子空间 $K_m(A,b)=\text{span}\{b,Ab,A^2b,...,A^{m-1}b\}$ 中的归一化 $b$

Arnoldi 算法迭代第一步

继续算 $K_m(A,b)$ 中的 $Ab$

$$w=Av_1= \frac1{\sqrt{11}} \begin{bmatrix} 1&1&1\\ 1&-1&1\\ 1&1&-1 \end{bmatrix} \begin{bmatrix} 3\\1\\1 \end{bmatrix} = \frac1{\sqrt{11}} \begin{bmatrix} 5\\3\\3 \end{bmatrix}$$

对 $v_1$ 正交化，先算投影系数

$$h_{11}=v_1^\top Av_1=v_1^\top w= \frac1{\sqrt{11}} \begin{bmatrix} 3&1&1 \end{bmatrix} \frac1{\sqrt{11}} \begin{bmatrix} 5\\3\\3 \end{bmatrix}= \frac{21}{11}$$

然后消掉 $v_1$ 方向的投影：

$$w \leftarrow w-h_{11}v_1$$

代入：

$$w= \frac1{\sqrt{11}} \begin{bmatrix} 5\\3\\3 \end{bmatrix} - \frac{21}{11}\cdot \frac1{\sqrt{11}} \begin{bmatrix} 3\\1\\1 \end{bmatrix}= \frac{1}{11\sqrt{11}} \begin{bmatrix} -8\\12\\12 \end{bmatrix}$$

求第二个基向量，先算它的范数：

$$h_{21}=v^{\top}_2Av_1=v^{\top}_2v_2= \|w\|= \frac1{11\sqrt{11}} \sqrt{(-8)^2+12^2+12^2} =\frac{4\sqrt2}{11}$$

于是

$$v_2=\frac{w}{h_{21}}= \frac1{\sqrt{22}} \begin{bmatrix} -2\\3\\3 \end{bmatrix}$$

Arnoldi 算法迭代第二步

再继续算 $K_m(A,b)$ 中的 $A^2b$

$$w=A^2v_1=Av_2= \frac1{\sqrt{22}} \begin{bmatrix} 1&1&1\\ 1&-1&1\\ 1&1&-1 \end{bmatrix} \begin{bmatrix} -2\\3\\3 \end{bmatrix} = \frac1{\sqrt{22}} \begin{bmatrix} 4\\-2\\-2 \end{bmatrix}$$

先对 $v_1$ 正交化

$$h_{12}=v_1^\top Av_2=v_1^\top w=\frac1{\sqrt{11}} \begin{bmatrix} 3&1&1 \end{bmatrix}\frac1{\sqrt{22}} \begin{bmatrix} 4\\-2\\-2 \end{bmatrix}=\frac{4\sqrt2}{11}$$

再对 $v_2$ 正交化

$$h_{22}=v_2^\top Av_2=v_2^\top w= \frac1{22} \begin{bmatrix} -2&3&3 \end{bmatrix} \begin{bmatrix} 4\\-2\\-2 \end{bmatrix}=-\frac{10}{11}$$

更新：

$$w=Av_2-h_{12}v_1-h_{22}v_2=v_3^\top Av_2=0$$

所以

$$h_{32}=v_3^\top Av_2=\|w\|=0$$

这说明 Arnoldi 在第二步就停了。如果不为零，则会一直迭代到维度为 Krylov 重启长度。

这里 $h$ 是 Hessenberg 矩阵 $H_m$ 的矩阵元，它表示 Arnoldi 过程将一个大型的 $n\times n$ 的矩阵 $A$ 投影到一个小型 $m\times m$ 的 $H_m$ 矩阵上，因此有

$$H_m=V_m^\top AV_m$$

其中 $V_m$ 为正交基 $v_m$ 组成的变换矩阵，作用是将矩阵 $A$ 降维到 $m\times m$

$$V_m= \begin{bmatrix} | & | & & |\\ v_1 & v_2 &...& v_m\\ | & | & & |\\ \end{bmatrix}$$

假如这里更新 $w\ne0$，接下来要算 $A^3b$，也就是 $w=Av_3=v_4$。为了得到它在 $v_4$ 上的坐标，要依次减去它在 $v_1,v_2,v_3$ 上的分量，即 $h_{13},h_{23},h_{33}$。根据 $h$ 的定义，再得到 $h_{31},h_{43}$，再判断 $v_4$ 的坐标 $h_{43}$ 是否为零。

总结前面的 Arnoldi 算法，其正交化过程实际上就是将 $Av_j$ 拆成各低阶分量与 $j+1$ 阶分量之和，然后减去所有低价项并除以模 $h_{j+1,j}$，从而得到高维正交基 $v_{j+1}$

$$Av_j=\sum_{i=1}^jh_{ij}v_i+h_{j+1,j}v_{j+1}= \begin{bmatrix} v_1&v_2&\cdots&v_{j+1} \end{bmatrix} \begin{bmatrix} h_{1j}\\h_{2j}\\\vdots\\h_{j+1,j} \end{bmatrix}= V_{m+1} \begin{bmatrix} h_{1j}\\h_{2j}\\\vdots\\h_{j+1,j} \end{bmatrix}$$

把所有列并在一起

$$A \begin{bmatrix} v_1 & v_2 & \cdots & v_m \end{bmatrix} =AV_m= V_{m+1} \begin{bmatrix} h_{11} & h_{12} & \cdots & h_{1m} \\ h_{21} & h_{22} & \cdots & h_{2m} \\ h_{31} & h_{32} & \cdots & h_{3m} \\ \vdots & \vdots & & \vdots \\ h_{m+1,1} & h_{m+1,2} & \cdots & h_{m+1,m} \end{bmatrix} =V_{m+1}\bar H_m$$

注意这里的 $\bar H_m$ 是 $(m+1) \times m$ 的矩阵，而 $H_m$ 是 $m\times m$ 的矩阵

得到投影小矩阵

于是二维 Krylov 子空间上的 Hessenberg 小矩阵为

$$H_2= \begin{bmatrix} h_{11} & h_{12}\\ h_{21} & h_{22} \end{bmatrix} = \begin{bmatrix} \frac{21}{11} & \frac{4\sqrt2}{11}\\[4pt] \frac{4\sqrt2}{11} & -\frac{10}{11} \end{bmatrix}$$

而

$$V_2= \begin{bmatrix} | & |\\ v_1 & v_2\\ | & | \end{bmatrix}$$

其中

$$v_1=\frac1{\sqrt{11}} \begin{bmatrix} 3\\1\\1 \end{bmatrix}, \qquad v_2=\frac1{\sqrt{22}} \begin{bmatrix} -2\\3\\3 \end{bmatrix}$$

得到投影方程

$$H_2 \begin{bmatrix} y_1\\y_2 \end{bmatrix} =\beta e_1=\beta \begin{bmatrix} 1\\0 \end{bmatrix}$$

现在来解释一下这个投影方程的数学/物理意义。我们知道 $y=[y_1,y_2]$ 是精确解 $x$ 投影在 Krylov 子空间 $K_m(A,b)=\text{span}\{b,Ab\}$ 上的解（这里没有 $A^2b$ 是因为这个维度正交基的模为0）。$H_2$ 是完整系数矩阵 $A$ 投影在子空间上的系数小矩阵。$\beta$ 是第一个标准正交基 $v_1$ 的模，同样也是 $b$ 投影在子空间内第一个基上的模大小。

因此，在这个低维空间中，原 $Ax=b$ 同样对应了一个小矩阵方程组 $H_my=\beta e_1$。我们只需求解这个小矩阵方程组，得到 $y$，再空间变换回 $x=V_m y$，就能得到低维近似解 $x^{(m)}$。

由此可见，求解复杂大型矩阵方程组的核心思想，就是对其降维并求解低维下的投影方程组。如果观察到解收敛了，就认为求解成功。

最小化残差

前面提到求解降维后的投影方程组，这其实是 FOM（Full Orthogonalization Method），是 Krylov 最基础的版本。GMRES 是它的优化版，毕竟低维近似解也可能一直不收敛，导致求解的 FOM 矩阵方程组越来越大。

所以人们干脆不直接对投影方程组求解了，而是算

$$y=\arg\min_y\|\beta e_1-\bar{H}_my\|$$

注意这里相当于是找到 $m$ 阶残差的最小值

根据前面推导的 Arnoldi 公式

$$AV_m=V_{m+1}\bar H_m$$

且 $b=\beta v_1=V_{m+1}(\beta e_1)$

$$r_m=b-Ax^{(m)}=V_{m+1}(\beta e_1)-AV_my=V_{m+1}(\beta e_1)-V_{m+1}\bar H_my=V_{m+1}(\beta e_1-\bar H_my)$$

因此，所有高阶残差项都投影到了 $V_{m+1}$ 这个基里。只要我找到其坐标 $\|\beta e_1-\bar{H}_my\|$ 的最小值，就能保证残差最小。最终问题又被转化为了一个 $(m+1)\times m$ 的最小二乘问题：

找到 $y$ 使得小矩阵方程组 $H_my=\beta e_1$ 的残差的模 $\|\beta e_1-\bar{H}_my\|$ 最小

QR 分解

$QR$ 分解是将原矩阵分解为一个正交矩阵 $Q$（满足 $Q^\top Q=I$）和一个上三角矩阵 $R$ 的乘积（$A=QR$）。该方法常用于解线性最小二乘问题，核心在于利用正交矩阵的性质将复杂计算转化为简单计算。

所以进一步地，计算最小残差的时候，选择用 $QR$ 分解将 $\bar H_m$ 变成正交矩阵 $Q$ 与上三角矩阵 $R$ 的乘积。

根据前面的推导，矩阵 $\bar H_m=[h_1,h_2,\cdots ,h_{m}]$，且每个 $h_i$ 为一个 $m+1$ 维度的向量。通过 Gram-Schmidt 施密特正交化将列向量转化为一组标准正交基 $Q=[q_1,q_2,\cdots,q_{m+1}]$（就是 Arnold 算法的正交化过程）

因为在正交化的过程中，原向量 $h_i$ 可以表示为这组新基 $q_i$ 的线性组合，例如：

$$\begin{aligned} h_1&=r_{11}q_1\\ h_2&=r_{12}q_1+r_{22}q_2\\ h_3&=r_{13}q_1+r_{23}q_2+r_{33}q_3\\ \vdots\\ h_{m}&=r_{1,m}q_1+\cdots+r_{m,m}q_{m} \end{aligned} \quad\Rightarrow\quad \bar H_m= \begin{bmatrix} h_1&h_2&h_3&\cdots&h_{m} \end{bmatrix}= \begin{bmatrix} q_1&q_2&q_3&\cdots&q_{m+1} \end{bmatrix} \begin{bmatrix} r_{11}&r_{12}&r_{13}&\cdots\\ 0&r_{22}&r_{23}&\cdots\\ 0&0&r_{33}&\cdots\\ \vdots&\vdots&\vdots&r_{m,m}\\ 0&0&0&0 \end{bmatrix}=QR$$

这里的系数 $r_{ij}$ 恰好就是矩阵 $R$ 的元素，刚好形成了上三角矩阵，于是得到

• $Q=\begin{bmatrix} q_1&q_2&q_3&\cdots&q_{m+1} \end{bmatrix}，(m+1)\times(m+1)$ 维度的正交矩阵
• $R=\begin{bmatrix} r_{11}&r_{12}&r_{13}&\cdots\\ 0&r_{22}&r_{23}&\cdots\\ 0&0&r_{33}&\cdots\\ \vdots&\vdots&\vdots&r_{m,m}\\ 0&0&0&0\end{bmatrix}，(m+1)\times m$ 维度的上三角矩阵

为什么这里非要将 $\bar H_m$ 拆成 $QR$ 的形式呢？因为带入残差式中，得到

$$\|r_m\| = \|Q^\top r_m\| = \|Q^\top\beta e_1-QRy\| = \|Q^\top\beta e_1-Ry\|$$

令 $g=Q^\top\beta e_1$，我们要求残差最小值 $\min\|r_m\|$，也就是解最小二乘问题 $\min\|g-Ry\|$，找到一组 $y$ 使得 $\|g-Ry\|$ 最小。

将之展开成平方

$$\|g - R y\|^2 = \sum_{i=1}^{m+1} (g_i - (Ry)_i)^2$$

逐行展开 $Ry$

$$\begin{alignat}{3} (Ry)_1 &= r_{11}y_1 &&+ r_{12}y_2 &&+ \dots + r_{1m}y_m \\ (Ry)_2 &= &&\quad\, r_{22}y_2 &&+ \dots + r_{2m}y_m \\ \vdots\\ (Ry)_{m} &= 0&&+0&&+\cdots+r_{mm}y_m\\ (Ry)_{m+1} &= 0&&+0&&+\cdots+0 \end{alignat}$$

也就是说

$$\|g - R y\|^2 = \sum_{i=1}^{m+1} (g_i - (Ry)_i)^2= \underbrace{ \sum_{i=1}^{m} (g_i - (Ry)_i)^2 }_{\text{可以优化}} + \underbrace{ g_{m+1}^2 }_{\text{无法改变}}$$

最后这一项就是与 $y$ 无关的残差项，因此 $\min\|g-Ry\|$ 问题就变成了求解方程组

$$g_i-(Ry)_i=0\quad(y=1,2,\cdots,m)$$

总结起来，最小二乘在这里的表现形式为求解 $Ry=Q^\top\beta e_1$ 方程组，最小二乘 = 能拟合的维度全部拟合，不能拟合的维度变成残差

回到 GMRES 算法中，最终残差就是 $\|r_m\| = \|g-Ry\| = |g_{m+1}|$。因此，GMRES 的目标就是把 $g_{m+1}$ 压到最小，且 $QR$ 分解能够让最后求解 $Ry=Q^\top\beta e_1$ 方程组时利用上三角矩阵 $R$ 的优势快速回代求解。

QR 分解的算法实现

求解器实现 QR 分解的算法为旋转消元法，构造一个 $(m+1)×(m+1)$ 的矩阵 Givens：

例如作用在第 $j$ 和 $j+1$ 行：

$$G_j = \begin{bmatrix} 1 & & & \\ & \ddots & & \\ & & c & s \\ & & -s & c \\ & & & \ddots \end{bmatrix}$$

可以发现是一个 $2\times 2$ 的旋转矩阵嵌在大矩阵里（$c$ 是 $\cos$，$s$ 是 $\sin$）

假设你有一列 $h^{(m)}$：

$$h^{(m)} = \begin{bmatrix} h_1 \\ h_2 \\ \vdots \\ h_j \\ h_{j+1} \\ h_{j+2} \\ \vdots \end{bmatrix}$$

应用 Givens：

$$h^{(m)} \leftarrow G_j h^{(m)}$$

实际效果是：

$$\begin{cases} h_j \leftarrow c h_j + s h_{j+1} \\ h_{j+1} \leftarrow -s h_j + c h_{j+1} \\ \text{其它行不变} \end{cases}$$

根据 $(h_{j},h_{j+1})$ 的角度设置旋转矩阵，利用这种方式将 $h_{j+1}$ 消元为0。由于每一步都构建了一个 Givens，导致每一步都在进行旋转消元

算法流程如下

$$\boxed{ \begin{aligned} &1.\ r_0 = b - A x^{(0)},\quad v_1 = r_0/\|r_0\| \\ \\ &2.\ \text{for } j=1,\dots,m: \\ &\quad (a)\ w = A v_j \\ &\quad (b)\ \text{Arnoldi 正交化 → 得 } h_{ij} \\ &\quad (c)\ \text{形成新列 } h^{(j)} \\ &\quad (d)\ \text{应用所有已有 Givens } G_1,\dots,G_{j-1} \\ &\quad (e)\ \text{构造 } G_j \text{ 消掉 } h_{j+1,j} \\ &\quad (f)\ \text{同时更新 } g = Q^T(\beta e_1) \\ \\ &3.\ \text{解上三角： } R y = g_{1:m} \\ \\ &4.\ x^{(m)} = V_m y \\ \\ &5.\ \|r_m\| = |g_{m+1}| \\ \end{aligned} }$$

1. 通过 $x^{(0)}=0$，得到 $v_1$
2. 循环 Arnoldi 算法，计算 $v_i$ 和 $h_{ij}$，并组合成 $h^{(j)} (m\times 1向量)$
   1. 应用已有的 Givens，$h^{(j)} \leftarrow G_{j-1} \cdots G_1 \, h^{(j)}$
   2. 针对 $\begin{bmatrix} h_{j,j} \\ h_{j+1,j} \end{bmatrix}$，构造 $G_j$，使 $h_{j+1,j} \rightarrow 0$
   3. 更新 $h^{(j)}$ 完成三角化，$h^{(j)} \leftarrow G_j h^{(j)}$
   4. 同步更新 $g$，$g \leftarrow G_j \, g$，数学上 $g = Q^T (\beta e_1) = G_m \cdots G_1 (\beta e_1)$，但 $g$ 是逐步更新
3. 上三角矩阵迭代完毕，求解 $Ry=g_{1:m}$
4. 转换回原 $x$ 空间，$x^{(m)} = V_m y$
5. 计算残差 $\|r_m\|$ 是否满足收敛要求

重启长度与最大迭代次数

定义解释

• 重启长度（$m$）：这是单次构造子空间的最大维度。每当正交基达到 $m$ 个，Arnold 算法就会停下来，计算一次近似解，然后清空这 $m$ 个基，$x^{(0)} \leftarrow x^{(m)},r_0 = b - A x^{(m)},v_1=r_0/\|r_0\|$ 重新开始
• 最大迭代次数 (Max Iterations)：这是指算法累计允许计算的正交基总数。
  • 例如：设置重启长度 $m=30$，最大迭代次数为 300。这意味着算法最多可以进行 10 轮重启（$30\times 10=300$）

这里之所以要设置这两个常量，是因为

1. 残差精度不够，需要继续迭代
   我们只求出了 $m$ 个基下的最优解 $x^{(m)}$，但如果此时残差 $\|b-Ax^{(m)}\|$ 还未达到想要的精度（收敛要求），就必须继续找更多的基，继续求解。
2. 内存墙（强制重启）：
   理论上，如果不重启，一直增加基的个数（$m$ 越来越大），解会越来越准。但由于每个正交基都要存放在内存里，且新基要与之前所有基做正交化：
   • 如果 $m$ 太大，例如 1000，内存会爆
   • 计算量越来越大（$O(m^2)$），越往后越慢

具体流程为（$m=30$）：

1. 先算 30 个基，求出一个当前最牛的近似解 $x^{(30)}$，Krylov 子空间 $K_{30}(A,r_0)$
2. 重启：把这 30 个基扔掉（释放内存），以 $x^{(30)}$ 产生的残差作为新起点，Krylov 子空间 $K_{30}(A,新残差)$
3. 再算 30 个基，求出一个基于前面结果改进的 $x^{(60)}$
4. 如此反复，直到超过最大 Arnoldi 迭代次数

重启副作用

重启可能会严重降低收敛速度，因为每次 restart：

• 丢掉之前的 Krylov 信息
• 相当于"记忆被清空"

所以实际工程中策略为 GMRES(m) + 预条件（preconditioner）