求解器中的预处理器

Krylov 算法中的预处理器 Preconditioner

前言

Newton-Krylov 迭代法已经能求解绝大多数非线性电路方程组了，但遇到一些刚性电路（特征值量级差异巨大的电路，雅可比矩阵病态），一个器件时间常数纳秒级，一个器件时间常数毫秒级，且带有严重非线性，收敛还是会变得非常困难。

人们为了解决上述情况，提出了预处理器（Preconditioner）这个概念。预处理的本质定义就是

把原问题 Ax=bAx=b 变成一个更容易被 Krylov 收敛的问题。

数学原理

当求解器拿到一个大型方程组 Ax=bAx=b 的时候，它通常不会直接开解，而会先对方程组进行预处理

左预处理（常见）

M−1Ax=M−1bM^{-1}Ax=M^{-1}b

右预处理

AM−1y=b,x=M−1yAM^{-1}y=b,\quad x=M^{-1}y

通过上述处理，能够使 M−1AM^{-1}A 的谱性质更好，从而更容易收敛。预处理之所以能够加速收敛，是因为

Krylov 迭代的收敛速度，近似取决于其特征值分布情况。理想情况下 M−1A≈IM^{-1}A\approx I，那 Krylov 迭代一次就能收敛。

因此，预处理器对矩阵处理的目标，就是将特征值压缩到一块儿。所谓压缩，就是通过一个权重因子，将特征值都拽倒一个量级里来。例如一个器件时间常数为 10−910^{-9}，另一个为 10−310^{-3}，通过给第一个时间常数加权 1e6，让两者都在 1e-3 的量级里。由于这种作用是可逆的，且方程两边同时作用，因此不会对结果造成影响，只会令矩阵 AA 不再雅可比矩阵病态。

ADS 中的三个预处理器

DCP

全称是直流预处理器 (DC Preconditioner)，本质是用 DC 模型近似整个系统

考虑到电路中，电容电感属于高频动态效应，电路主干是电阻静态网络，高频是扰动，DC是主体。因此构建预处理矩阵

M≈ADCM\approx A_{DC}

所以

M−1A≈I+小扰动M^{-1}A\approx I+小扰动

这种情况下 Krylov 就很容易收敛。但当强非线性/强动态电路存在时，例如振荡器、高频反馈、光电耦合系统，将电路主干认为是静态网络就有问题了，相当于 I+大扰动I+大扰动。

DCP 也就是用 DC 近似，认为 ω=0\omega=0，电容开路，电感短路，先得到 ADCA_{DC}

但问题来了，要如何构建 MM，使其能既像 AA，又容易求 M−1M^{-1}

ILU 法

对 ADCA_{DC} 做：

ADC≈LUA_{DC} \approx LU

但只保留部分非零（Incomplete LU），得到 M=LUM=LU。每次需要 M−1vM^{-1}v 时，只需先求解

Ly=vLy=v

再求解

Ux=yUx=y

其中 x=M−1vx=M^{-1}v

对角法

只取对角

M=diag(ADC),M−1v=vdiagM=\text{diag}(A_{DC}),\quad M^{-1}v=\frac{v}{\text{diag}}

BSP

分块选择预处理器 (Block Select Preconditioner)，与 DCP 相比，区别在于处理 ADCA_{DC} 时，把矩阵分块化

ADC=[A11∗∗A22]A_{DC} = \begin{bmatrix} A_{11} & * \\ * & A_{22} \end{bmatrix}

构造

M=[A1100A22]M = \begin{bmatrix} A_{11} & 0 \\ 0 & A_{22} \end{bmatrix}

然后对每个块单独 LU，或单独求逆。这种方式比对角法和 ILU 更稳定可控

SCP

舒尔补预处理器 (Schur-Complement Preconditioner)，与前两者相比，更加强大与稳定

这个预处理器是在 BSP 的基础上，也就是把系统分块后，只对难解的部分精确处理

我们有分块系统：

[ABCD][x1x2]=[b1b2]\begin{bmatrix} A & B \\ C & D \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \end{bmatrix} = \begin{bmatrix} b_1 \\ b_2 \end{bmatrix}

展开就是：

{Ax1+Bx2=b1Cx1+Dx2=b2\begin{cases} A x_1 + B x_2 = b_1 \\ C x_1 + D x_2 = b_2 \end{cases}

既然 AA 分块好解，那就先求解它，把 x1x_1 消掉

从第一行解 x1x_1

Ax1=b1−Bx2A x_1 = b_1 - B x_2

代入第二行

CA−1(b1−Bx2)+Dx2=b2C A^{-1}(b_1 - B x_2) + D x_2 = b_2

整理得到

(D−CA−1B)x2=b2−CA−1b1(D - C A^{-1} B) x_2 = b_2 - C A^{-1} b_1

由此，得到了比原系统更好的 S=D−CA−1BS=D - C A^{-1} B，规模更好，条件值也更好，整个方程组变得更容易求解了。

注意事项

后两个预处理器在做分块的时候，2×22\times 2 的分块只是一个例子，实际情况可能远不止 2×22\times 2

但因为所有多块分解，本质都可以递归成 2×22\times 2 的 SCP

假设系统是：

[ABECDFGHK][x1x2x3]=[b1b2b3]\begin{bmatrix} A & B & E \\ C & D & F \\ G & H & K \end{bmatrix} \begin{bmatrix} x_1\\x_2\\x_3 \end{bmatrix} = \begin{bmatrix} b_1\\b_2\\b_3 \end{bmatrix}

Step 1：先把前两块当一个整体

[ABCD]\begin{bmatrix} A & B \\ C & D \end{bmatrix}

Step 2：对它做 Schur 补（消掉 x1x_1）

得到：

S2=D−CA−1BS_2 = D - C A^{-1} B

Step 3：再和第三块做 Schur 补

最终得到：

S3=K−[G H][ABCD]−1[EF]S_3 = K - [G \ H] \begin{bmatrix} A & B \\ C & D \end{bmatrix}^{-1} \begin{bmatrix} E\\F \end{bmatrix}

在电光链路中，可以自然分成：

  • 电驱动（Tx）
  • 光纤/传播
  • 探测（Rx）
  • 数字处理

对应矩阵：

[ATx∗∗∗∗Afiber∗∗∗∗APD∗∗∗∗ADSP]\begin{bmatrix} A_{\text{Tx}} & * & * & * \\ * & A_{\text{fiber}} & * & * \\ * & * & A_{\text{PD}} & * \\ * & * & * & A_{\text{DSP}} \end{bmatrix}

我们可以对"容易部分"做预处理，对"难耦合部分"做 Schur 补