线性求解的离散伴随

形式陈述

设参数 $\theta\in\mathbb R^p$ ，实矩阵 $A(\theta)\in\mathbb R^{n\times n}$ 、右端 $b(\theta)\in\mathbb R^n$ 与标量目标 $J(\theta,u)$ 都为 $C^1$ 。在所考察的 $\theta$ 处，假设 $A(\theta)$ 可逆。定义线性方程组的精确解和约化目标

A(\theta)u(\theta)=b(\theta),\qquad \widehat J(\theta)=J(\theta,u(\theta)).

行列式的连续性保证该点的某个邻域中 $A$ 仍可逆， $u=A^{-1}b$ 在此邻域为 $C^1$ 。本页计算的是这个解映射的导数；算法实现使用求解而不显式形成逆矩阵。

切向方程与伴随恒等式

固定参数方向 $h$ ，用点号表示沿 $h$ 的导数。微分原方程得到

\dot A\,u+A\dot u=\dot b,\qquad A\dot u=\dot b-\dot A\,u.

由链式法则，

D\widehat J(\theta)[h] =\partial_\theta J(\theta,u)[h]+g^{\mathsf T}\dot u, \qquad g=\nabla_uJ(\theta,u).

这里 $\partial_\theta J$ 固定 $u$ ，不能漏掉目标对参数的直接依赖。定义伴随向量为转置系统的解

A^{\mathsf T}\lambda=g.

于是无需显式算出 $\dot u$ 就有

g^{\mathsf T}\dot u =\lambda^{\mathsf T}A\dot u =\lambda^{\mathsf T}(\dot b-\dot A\,u),

从而

\boxed{ D\widehat J(\theta)[h] =\partial_\theta J(\theta,u)[h] +\lambda^{\mathsf T}\bigl(Db(\theta)[h]-DA(\theta)[h]u\bigr). }

取各坐标方向，梯度分量为

\frac{\partial\widehat J}{\partial\theta_j} =\frac{\partial J}{\partial\theta_j} +\lambda^{\mathsf T} \left(\frac{\partial b}{\partial\theta_j} -\frac{\partial A}{\partial\theta_j}u\right).

这就是离散伴随：对已经给定的有限维代数方程求导，再解其转置系统。它不需要先对连续 PDE 求伴随，也没有自动证明连续伴随离散化与这套方程完全相同。

有限求解过程与因子复用

输入包括 $A,b,J$ 及其所需导数接口、参数值，以及数值实现使用的残差容差。执行顺序是：先解 $Au=b$ ；在所得状态计算 $J,g,\partial_\theta J$ ；再解 $A^{\mathsf T}\lambda=g$ ；最后计算上述方向收缩或全部参数分量。输出目标值、导数、两次求解的残差及状态。奇异、不可用主元或非有限值应返回失败；采用有迭代预算的求解器时，预算耗尽而未满足标准也应明确返回未收敛状态。只有各步骤完成才返回成功，数值误差仍需结合条件性解释。

对稠密可逆矩阵，用带置换 LU得到 $PA=LU$ 。正向求解是

Ly=Pb,\qquad Uu=y.

因为 $A^{\mathsf T}=U^{\mathsf T}L^{\mathsf T}P$ ，转置求解可以复用同一组因子：

U^{\mathsf T}z=g,\qquad L^{\mathsf T}w=z,\qquad \lambda=P^{\mathsf T}w.

置换在转置求解的末端施加，不能照抄正向的 $Pb$ 顺序。一次稠密因子分解需 $O(n^3)$ 工作和 $O(n^2)$ 存储，每个普通或转置右端需 $O(n^2)$ 工作；不必重新分解 $A^{\mathsf T}$ 。

对一个标量目标，伴随右端数不随参数个数 $p$ 增长。这只是在数求解次数：导数接口、参数收缩和输出 $p$ 个分量仍有成本。若只要一个方向导数，解一次切向方程同样自然；若有多个独立目标，每个目标通常有自己的伴随右端。

直觉

切向方法先问每个参数扰动如何改变整份状态 $\dot u$ ，再把状态变化投到目标上。伴随方法先把目标敏感度 $g$ 通过转置系统传回来，得到 $\lambda$ ，随后任何参数方向只需与方程变化 $\dot b-\dot A\,u$ 做内积。

转置来自内积恒等式 $g^{\mathsf T}A^{-1}q=(A^{-\mathsf T}g)^{\mathsf T}q$ 。它并不意味着原方程是对称的；即便 $A$ 对称， $u$ 和 $\lambda$ 的右端通常不同，也不能把两个向量混为一谈。

图中两条求解路径共用同一组因子。状态 $u$ 既用于计算目标的状态梯度 $g$ ，也进入右侧参数收缩；直接参数项必须另外相加。

例子与边界

二维系统的三种核对

取

A(\theta)=\begin{pmatrix}\theta+2&1\\1&2\end{pmatrix}, \qquad b=\begin{pmatrix}1\\0\end{pmatrix}, \qquad J(\theta,u)=\tfrac12u^{\mathsf T}u.

在 $\theta=0$ ，原方程、切向方程与伴随方程分别给出：

对象	方程右端	解
原状态 $u$	$b=(1,0)^{\mathsf T}$	$(2/3,-1/3)^{\mathsf T}$
切向状态 $\dot u$	$-\dot A u=(-2/3,0)^{\mathsf T}$	$(-4/9,2/9)^{\mathsf T}$
伴随 $\lambda$	$g=u=(2/3,-1/3)^{\mathsf T}$	$(5/9,-4/9)^{\mathsf T}$

取参数方向 $h=1$ 。切向计算得到

u^{\mathsf T}\dot u =\frac23\left(-\frac49\right) -\frac13\frac29=-\frac{10}{27}.

伴随计算只取 $A$ 左上角的参数变化，因此

-\lambda^{\mathsf T}\dot A\,u =-\frac59\frac23=-\frac{10}{27}.

第三种核对直接解出

u(\theta)=\frac1{3+2\theta}\begin{pmatrix}2\\-1\end{pmatrix}, \qquad \widehat J(\theta)=\frac5{2(3+2\theta)^2}, \qquad \widehat J'(0)=-\frac{10}{27}.

这里 $A$ 对称而 $\lambda\ne u$ ；两次求解是同一系数矩阵、不同右端。

有限次迭代是另一个映射

考虑标量系统 $(\theta+2)u=1$ ，目标仍为 $J=\tfrac12u^2$ 。精确解目标为

\widehat J(\theta)=\frac1{2(\theta+2)^2}, \qquad \widehat J'(0)=-\frac18.

从 $u_0=0$ 运行一步固定步长 Richardson 迭代，

u_1=u_0+\tfrac13\bigl(1-(\theta+2)u_0\bigr)=\tfrac13.

这一步程序输出与 $\theta$ 无关，所以对实际程序目标 $\tfrac12u_1^2$ 求导，结果为 $0$ 。若把尚未收敛的 $u_1$ 填进隐式公式，再精确解 $(\theta+2)\widetilde\lambda=u_1$ ，则在零点得到

-\widetilde\lambda u_1=-\frac1{18}.

$0$ 、 $-1/18$ 和 $-1/8$ 分别对应有限程序导数、代入近似状态的隐式估计、精确解导数。此迭代在零点邻域内是收敛的；差异来自只运行了一步，不能归因于选了一个发散方法。对有限迭代反向传播可以正确计算该程序的导数，但它的目标并不是精确解映射。

两个残差与条件性

对近似状态 $\widetilde u$ ，令 $r=b-A\widetilde u$ ，则

u-\widetilde u=A^{-1}r.

对固定的精确右端 $g$ ，若 $r_\lambda=g-A^{\mathsf T}\widetilde\lambda$ ，则

\lambda-\widetilde\lambda=A^{-\mathsf T}r_\lambda.

这些误差通过乘积进入导数，因此小残差必须结合线性系统条件性解释。实际 $g=\nabla_uJ(\theta,u)$ 还可能随状态变化；若在 $\widetilde u$ 上计算它，伴随右端也有误差，不能只报告转置求解器相对于这个近似右端的残差。伴随公式减少了右端数，并不自动改善问题条件性。

推论与应用

把线性求解看作 $S(A,b)=A^{-1}b$ 原语，给定输出敏感度 $g$ 后，其反向自动微分规则可写成

\bar b=\lambda,\qquad \bar A=-\lambda u^{\mathsf T},\qquad A^{\mathsf T}\lambda=g.

这是因为 $\lambda^{\mathsf T}(db-dA\,u)$ 中矩阵项恰为 Frobenius 配对 $\langle-\lambda u^{\mathsf T},dA\rangle_F$ 。它只描述通过求解输出的贡献；目标若还直接依赖 $A$ 或 $b$ ，相应偏导要另外累加。参数化稀疏矩阵时可以直接计算 $\lambda^{\mathsf T}(\partial_j A)u$ ，不一定显式生成稠密 $\bar A$ 。

该接口可用于参数估计、离散约束优化与可微仿真。若后续还需二阶方向信息，Hessian–向量积提供组合求导的组织方式，但须增加所需的二阶光滑性，并让各次原方程和转置求解达到相应精度。本页只证明实数、非奇异线性系统的一阶导数；奇异系统的选解规则和复数求导需另行规定。

参考资料

Mike B. Giles, An Extended Collection of Matrix Derivative Results for Forward and Reverse Mode Algorithmic Differentiation, Oxford Numerical Analysis Group Report NA-08/01, 2008，§2.1 的矩阵内积约定、§2.3.1（p.5）的 matrix inverse product 前向与反向公式。
Mike B. Giles and Niles A. Pierce, “An Introduction to the Adjoint Approach to Design”, Flow, Turbulence and Combustion 65, 2000, pp. 393–415，§2.2 的线性对偶与求解次数、§2.4 的显式参数依赖；引言区分连续与离散伴随。

形式陈述 ​

切向方程与伴随恒等式 ​

有限求解过程与因子复用 ​

二维系统的三种核对 ​

有限次迭代是另一个映射 ​

两个残差与条件性 ​

反向链接