伴随法

伴随法（adjoint method）通过求解一个与状态方程转置相关的方程，计算标量目标对许多参数的局部导数。它把“对每个参数求一次状态变化”，改成“对一个目标求一次伴随变量”。当参数很多、关心的目标较少时，这种重排可以大幅减少线性方程求解次数。

状态先由方程决定

设状态 $x \in ℝ^{n}$ 、参数 $p \in ℝ^{m}$ 在所讨论点的一个邻域内满足连续可微的残差方程 $R (x, p) = 0.$ 目标是在该邻域内连续可微的标量函数 $J (x, p)$ 。在所讨论解处，假设状态雅可比矩阵 $R_{x}$ 可逆；隐函数定理保证该点附近可把状态看成参数的可微函数。对第 $j$ 个参数微分： $R_{x} \frac{\partial x}{\partial p_{j}} + R_{p_{j}} = 0, \frac{d J}{d p_{j}} = J_{p_{j}} + J_{x} \frac{\partial x}{\partial p_{j}} .$ 直接灵敏度法为每个参数解一次状态导数方程。^[1]

定义列向量 $λ$ 为伴随方程的解： $R_{x}^{𝖳} λ = J_{x}^{𝖳} .$ 代入全导数并利用状态导数方程，即得 $\frac{d J}{d p_{j}} = J_{p_{j}} - λ^{𝖳} R_{p_{j}} .$ 负号来自本文把状态方程写成 $R = 0$ 的约定。若残差符号或拉格朗日函数符号改变，伴随变量的符号约定也可能改变；最终的目标导数必须一致。

二元线性方程的手算

令 $p > - 1 / 2$ ，状态由 $\begin{matrix} A (p) x & = b, \\ A (p) & = (\begin{matrix} 2 + p & 1 \\ 1 & 2 \end{matrix}), \\ b & = (1, 0)^{𝖳}, J = x_{2} . \end{matrix}$ 决定。行列式为 $3 + 2 p > 0$ ，矩阵可逆。直接解可得 $x_{1} = \frac{2}{3 + 2 p}, x_{2} = - \frac{1}{3 + 2 p}, \frac{d J}{d p} = \frac{2}{(3 + 2 p)^{2}} .$

在 $p = 0$ ，先由原方程求出 $x = (2 / 3, - 1 / 3)^{𝖳}$ 。目标梯度为 $J_{x} = (0, 1)$ ，伴随方程为 $(\begin{matrix} 2 & 1 \\ 1 & 2 \end{matrix}) λ = (\begin{matrix} 0 \\ 1 \end{matrix}), λ = (- 1 / 3, 2 / 3)^{𝖳} .$ 本例 $R = A (p) x - b$ ，所以 $R_{p} = A_{p} x = (x_{1}, 0)^{𝖳}$ ， $J_{p} = 0$ 。因此 $\frac{d J}{d p} = - λ_{1} x_{1} = - (- \frac{1}{3}) \frac{2}{3} = \frac{2}{9},$ 与直接微分一致。

目标J等于负一除以3加2p的曲线，以及p为0处斜率九分之二的切线 — 切线斜率为九分之二。它描述当前点的局部变化，离开该点后曲线一般不与切线重合。

多参数中的计算节省

线性方程 $A (p) x = b (p)$ 对应的公式为 $A^{𝖳} λ = J_{x}^{𝖳}, \frac{d J}{d p_{j}} = J_{p_{j}} + λ^{𝖳} (b_{p_{j}} - A_{p_{j}} x) .$ 一次状态求解和一次伴随求解之后，每个参数的导数由内积计算。若有很多不同的标量目标，通常每个目标都要有相应伴随求解；伴随法的优势取决于参数数与目标数的相对大小。

变式把参数放进载荷，令 $A = (\begin{matrix} 2 & 1 \\ 1 & 2 \end{matrix})$ 固定、 $b (p) = (1, p)^{𝖳}$ 。仍取 $J = x_{2}$ ，伴随不变，而 $b_{p} = (0, 1)^{𝖳}$ ，故导数为 $λ_{2} = 2 / 3$ 。直接解 $x_{2} = (- 1 + 2 p) / 3$ 也给出相同结果。这检查了参数进入矩阵与进入右端时的符号差别。

可运行的线性伴随程序

下面的目标为 $J = c^{𝖳} x$ ，其中 $c$ 固定。输入给出当前参数点的 $A, b$ 及每个参数的 $A_{p_{j}}, b_{p_{j}}$ ，返回目标值与梯度。函数不通过有限差分近似这些导数，不改写输入；奇异矩阵由线性求解器报错。

Python 3

import numpy as np

def linear_adjoint(A, b, c, dA, db):
    A, b, c = (np.asarray(v, dtype=float) for v in (A, b, c))
    dA = np.asarray(dA, dtype=float)
    db = np.asarray(db, dtype=float)
    n = b.size
    if b.shape != (n,) or c.shape != (n,) or A.shape != (n, n) or n == 0:
        raise ValueError("inconsistent state dimensions")
    if dA.ndim != 3 or dA.shape[1:] != (n, n):
        raise ValueError("dA must have shape (parameters,n,n)")
    if db.shape != (dA.shape[0], n):
        raise ValueError("db must have shape (parameters,n)")
    if not all(np.isfinite(v).all() for v in (A, b, c, dA, db)):
        raise ValueError("inputs must be finite")
    x = np.linalg.solve(A, b)
    adjoint = np.linalg.solve(A.T, c)
    gradient = (db - np.einsum('kij,j->ki', dA, x)) @ adjoint
    return float(c @ x), gradient

A = np.array([[2., 1.], [1., 2.]])
dA = np.array([[[1., 0.], [0., 0.]]])
value, gradient = linear_adjoint(A, [1., 0.], [0., 1.], dA, [[0., 0.]])
print(value, gradient)  # -0.333333... [0.222222...]
assert np.allclose(gradient, [2 / 9])
_, rhs_gradient = linear_adjoint(A, [1., 0.], [0., 1.], np.zeros((1, 2, 2)), [[0., 1.]])
assert np.allclose(rhs_gradient, [2 / 3])

对稠密 $n$ 阶系统和 $m$ 个参数，这个实现的时间为 $O (n^{3} + m n^{2})$ ，存储全部导数张量需要 $O (m n^{2})$ 空间。大型实现通常只计算所需的导数作用与内积，并复用矩阵分解，避免显式存储导数张量。

可微性与检验

伴随公式给的是当前模型、当前解分支上的局部导数。若状态雅可比奇异、求解分支发生切换，或目标含不可微操作，上述可微推导可能不成立。迭代方程未充分收敛时，残差也会影响导数质量。

可选若干非极小步长，用中心差分 $[J (p + h) - J (p - h)] / (2 h)$ 与梯度交叉核对；二者应在截断误差与舍入误差之间出现一致区间。不能仅把 $h$ 取得极小，因为相近数相减会放大浮点误差。对于时变微分方程，伴随还涉及终端条件与反向时间积分，不能把这里的静态转置求解公式直接当作完整时变算法。

参考文献

↑ Steven G. Johnson：Notes on Adjoint Methods for 18.335，2006年建立、2021年更新，§1—§3，线性与非线性方程的伴随灵敏度推导。

[1] Steven G. Johnson：Notes on Adjoint Methods for 18.335，2006年建立、2021年更新，§1—§3，线性与非线性方程的伴随灵敏度推导。

[1]