企业官网建设流程全解析

上周我把一套肿瘤生长模型的伴随灵敏度分析跑通，顺手拿它做了个时空放射治疗剂量优化的演示，结果好几个朋友看完第一反应都一样：既然要算梯度，直接对每个参数做扰动不就行了吗？这个问题我也想过，直到算了一笔账：如果控制变量是四维时空里每个网格点的剂量率，数量级至少是几十万到几百万，有限差分要重复跑几十万次 PDE，而伴随方法一次正演加一次伴随，大约两倍正演成本就能回传全梯度。更妙的是，这个结论和参数个数无关，五十个参数和五十万个参数，花的力气几乎一样。

这篇文章会把这条技术路线完整拆开：先解释为什么时空放疗优化需要灵敏度分析，再讲我怎么选肿瘤生长模型，然后一步步推出伴随方程，最后给出一个能在 Matlab 里跑起来的实现框架。中间穿插我踩过的几个坑，尤其是那种会让伴随梯度看起来对、实际上错得离谱的边界条件问题。代码不是完整工程，但架子是能跑的，读者拿到后可以直接往自己的模型上套。

1. 放疗剂量怎么“长”到组织里：从静态计划到时空优化的痛点

1.1 传统计划优化为什么突然不香了

调强放疗（IMRT）或者 VMAT 的计划优化，本质上是在一个固定解剖图像上找最优的通量图。设计变量是束斑权重、机架角度这些，顶多几百个，用梯度类算法处理非常成熟。目标函数通常写成功效函数加剂量体积约束：让肿瘤覆盖 95% 体积的处方剂量，同时把脊髓、腮腺这些危及器官的最大剂量压下去。

但问题是，计划做完患者躺上台，这一套剂量分布就固定了。实际治疗要持续好几周，肿瘤在退缩，正常组织在移动，患者的体重和体型也在变。第三周的解剖和第一周已经不是同一个解剖，固定剂量画在旧图像上，自然不是最优的。现在临床上不少中心在做在线自适应放疗，每治疗一次重新 CT 一次，重新优化一次计划。每次重新优化都需要重新算梯度，如果设计变量只有几百个还好，可一旦动了真格——把剂量率也变成空间和时间上连续变化的场——变量数就远远不是几百个能打住的了。

1.2 把时间维加进来，设计变量瞬间爆炸

时空放疗优化，用大白话说就是让剂量率 (u(x,t)) 在空间的每个位置、治疗过程的每个时刻都能独立调节。你想在肿瘤侵袭前沿多加一点剂量，在内脏敏感区附近减一点剂量，而且这个决策还要跟着肿瘤生长的状态走。数学上这个模型挺漂亮，工程上代价很高：假如空间网格取 64×64，治疗过程分成 40 个时间步，控制变量就有 (64 \times 64 \times 40 \approx 16) 万个。如果考虑三维，把 64 换成 32×32×32，那就是三百万量级。

在这种维度下，你要算梯度，首先想到的往往是有限差分。但有限差分需要对每个变量扰动一次、重新求解一次 PDE。一次正演算下来可能是几个 CPU 小时，乘以十几万个变量，这个优化就彻底没法玩了。伴随灵敏度分析就是冲着这个问题来的：它可以一次拿到所有变量对应的梯度。

1.3 有限差分：简单但昂贵的穷人版灵敏度

有限差分算梯度的公式很简单：

[ \frac{\partial J}{\partial u_i} \approx \frac{J(u+\varepsilon e_i) - J(u)}{\varepsilon} ]

实现零门槛，而且每个分量可以并行跑。问题是维度一高就露馅，另外扰动步长 (\varepsilon) 的选取也有讲究：选大了有截断误差，选小了被浮点噪声吞掉。我见过有人用复杂步长（complex-step）方法缓解精度问题，但在 PDE 求解器里实现复杂步长，等于把整个求解器重写一遍，不划算。

有限差分不是不能用，它非常适合做梯度验证——后面我会用它来检查伴随梯度到底对不对。但做优化主体，尤其是高维时空优化，真的扛不住。

1.4 伴随方法：一场从出口逆流而上的追踪

伴随方法的直觉可以用一条河来打比方。上游很多支流都在向河里排污，你想知道哪条支流对入海口水质影响最大。常规思路是一条条支流轮流关停，看入海口变化，这是有限差分思路。伴随方法则是从入海口逆流而上，测每一点对出口浓度的敏感性——只走一遍，所有支流的影响都清楚了。

在 PDE 框架里，这个"逆流而上"就是求解一个伴随方程。它和时间上反向传播的敏感性信息绑在一起，目标函数对每个控制变量的梯度全部由一次反向求解得到。后面我会把公式逐步展开，现在先记住这个直觉：伴随方法把"控制变量个数"这个乘数从算法复杂度里彻底消掉了。

2. 我的模型选择：反应扩散方程加线性二次杀伤项的数学画像

2.1 为什么用反应扩散方程描述肿瘤生长

我一开始也想偷懒，搞个常微分方程，比如 Logistic 增长：

[ \frac{dc}{dt} = \rho c\left(1-\frac{c}{K}\right) ]

但这个模型里没有空间，没法回答"哪个位置的剂量该调高、哪个位置该调低"这种问题。时空优化必须要一个带空间结构的模型，反应扩散方程是自然选择：

[ \frac{\partial c}{\partial t} = D_c \nabla^2 c + \rho c\left(1-\frac{c}{K}\right) ]

第一项是扩散，代表肿瘤细胞向周围组织的侵袭；第二项是反应，代表增殖，但受到环境承载能力 (K) 的限制。这个方程有个很响亮的名字，Fisher-Kolmogorov 方程，它能产生行波解，肿瘤边界会以有限速度向外扩张，这在形态上很像实体瘤的生长过程。

参数 (D_c) 是扩散系数，(\rho) 是增殖率，(K) 是环境容纳量。实际建模时，这三个数可以通过影像或者实验数据标定。如果是脑瘤模型，(D_c) 可能取 0.1–1 mm²/day；(\rho) 取 0.1–0.5 /day；(K) 则像是一个归一化上限，把细胞密度压住。

2.2 把放射治疗的杀伤效应写进方程

放射生物学里最常用的细胞存活模型是线性二次模型（LQ 模型），单次剂量 (d) 对应的存活分数为：

[ SF(d) = \exp\left(-\alpha d - \beta d^2\right) ]

把这个效应放进 PDE，不能直接用单次剂量，因为我们是在连续时间上优化剂量率 (u(x,t))。对一小段时间切片来说，细胞因辐射死亡的速率正比于 ((\alpha u + \beta u^2)c)。于是肿瘤方程变成：

[ \frac{\partial c}{\partial t} = D_c \nabla^2 c + \rho c\left(1-\frac{c}{K}\right) - \left(\alpha_c u + \beta_c u^2\right)c ]

同时我还加了正常组织方程：

[ \frac{\partial n}{\partial t} = D_n \nabla^2 n - \left(\alpha_n u + \beta_n u^2\right)n ]

正常组织初始密度归一化为 1，没有增殖项，辐射会把它的密度打下去。(\alpha_c, \beta_c) 是肿瘤的放射敏感性参数，(\alpha_n, \beta_n) 是正常组织的参数。通常肿瘤的 (\alpha/\beta) 比较小，正常组织偏高，这里可以体现出分次剂量优化的空间。

边界条件我最常用 Neumann 零通量：

[ \frac{\partial c}{\partial \nu} = 0,\quad \frac{\partial n}{\partial \nu} = 0 ]

意思是细胞不会穿过计算区域边界跑出去。这个选择很关键，因为伴随方程的边界条件必须和正演一致，否则梯度会出问题，后面我会专门讲这个坑。

2.3 目标函数：既能指导优化又适合求导

优化目标很直白：治疗结束后，肿瘤细胞总数尽可能少，正常组织损伤尽可能小，剂量不能无限大。我采用平方和形式：

[ J(u) = \omega_c \int_\Omega c(x,T)^2 dx + \omega_n \int_\Omega \left(1 - n(x,T)\right)^2 dx + \gamma \sum_k |u_k|^2 ]

第三项是正则项，防止优化器把剂量推到无穷大。平方项的选择是有原因的：如果用一次方，终端伴随变量 (\lambda(T)=\omega_c) 是个常数，对高密度区域和低密度区域一视同仁；用平方项，(\lambda(T)=2\omega_c c(T))，高密度区域的惩罚更重，更符合"别让肿瘤长大"的直觉。另外平方函数光滑，梯度在零点附近依然可导，数值上更友好。

权重 (\omega_c, \omega_n, \gamma) 需要反复调。需要注意到 (c) 和 (n) 的量纲不同，(c) 可能是 0.01 到 1 的数，(n) 初始值是 1，三者放在同一个目标函数里，如果不做归一化，梯度很容易被某一项支配。后面调参经验里我再展开。

3. 伴随灵敏度分析的推导：一次正演一次伴随解决梯度

3.1 从拉格朗日乘子法看伴随方程怎么冒出来的

我不想绕弯子，直接给出推导路径。把前向方程写成简洁形式：

[ \frac{\partial c}{\partial t} = D_c\nabla^2 c + f(c,u) ]

[ \frac{\partial n}{\partial t} = D_n\nabla^2 n + g(n,u) ]

其中

[ f(c,u) = \rho c\left(1-\frac{c}{K}\right) - \left(\alpha_c u + \beta_c u^2\right)c ]

[ g(n,u) = -\left(\alpha_n u + \beta_n u^2\right)n ]

要计算 (J) 对 (u) 的梯度，最干净的办法是构造拉格朗日函数，把 PDE 约束乘上伴随乘子（拉格朗日乘子）(\lambda_c, \lambda_n)：

[ L = J + \int_0^T\int_\Omega \lambda_c\left(\frac{\partial c}{\partial t} - D_c\nabla^2 c - f(c,u)\right)dxdt + \int_0^T\int_\Omega \lambda_n\left(\frac{\partial n}{\partial t} - D_n\nabla^2 n - g(n,u)\right)dxdt ]

接下来对状态变量 (c) 做变分。积分里的时间导数项用分部积分，把时间导数挪到 (\lambda_c) 上；空间 Laplacian 项也用分部积分，把算子挪到 (\lambda_c) 上。整理后，为了消掉所有含 (\delta c) 的背景项，伴随变量必须满足方程：

[ -\frac{\partial \lambda_c}{\partial t} = D_c\nabla^2 \lambda_c + \frac{\partial f}{\partial c}\lambda_c ]

终端条件是：

[ \lambda_c(T) = \frac{\partial J}{\partial c(T)} = 2\omega_c c(T) ]

同样地，对 (n) 有：

[ -\frac{\partial \lambda_n}{\partial t} = D_n\nabla^2 \lambda_n + \frac{\partial g}{\partial n}\lambda_n ]

[ \lambda_n(T) = \frac{\partial J}{\partial n(T)} = -2\omega_n\left(1-n(T)\right) = 2\omega_n\left(n(T)-1\right) ]

这里 (\partial f/\partial c)、(\partial g/\partial n) 都是从正演轨迹上的已知状态算出来的雅可比项，所以伴随方程虽然是 PDE，但它对 (\lambda_c,\lambda_n) 是线性的。换句话说，不管原来正演方程有多非线性，伴随方程总是线性的，这是伴随方法数值上特别划算的另一个原因。

3.2 离散时间伴随：手把手推一个时间层

上面是连续形式，但写 Matlab 代码时用离散更直接。前向用隐式欧拉处理扩散、显式处理反应，时间步为 (\Delta t)：

[ c^{k+1} = M_c^{-1}\left[c^k + \Delta t f(c^k, u^k)\right] ]

其中 (M_c = I - \Delta t D_c A)，(A) 是空间离散的 Laplacian 矩阵。如果把这一步看成神经网络里的一个算子：

[ c^{k+1} = F(c^k, u^k) ]

那么伴随递推就是逆向模式自动微分（反向传播）的连续版。给定终端乘子 (\lambda_c^{N+1})，反向递推公式是：

[ \lambda_c^k = \left(\frac{\partial F}{\partial c^k}\right)^T \lambda_c^{k+1} ]

代入 (F) 的表达式，并利用 (M_c) 的对称性，可以化简成两步：

[ q_c = M_c^{-1}\lambda_c^{k+1} ]

[ \lambda_c^k = q_c + \Delta t \left(\frac{\partial f}{\partial c^k}\right)^T q_c ]

对应的梯度贡献是：

[ \frac{\partial J}{\partial u^k} = \Delta t \left(\frac{\partial f}{\partial u^k}\right)^T q_c + \Delta t \left(\frac{\partial g}{\partial u^k}\right)^T q_n + 2\gamma u^k ]

这里的关键点是：伴随递推不需要再解任何非线性问题，每一步只是一个线性稀疏矩阵求解和一堆点乘。而且不定 (u) 有多少个自由度，反向走一遍，所有梯度分块都出来了。

3.3 梯度公式和控制变量维度无关的依据

很多人第一次看到伴随方法，会怀疑是不是漏了什么东西：真的一次反向求解就能拿到所有变量的梯度？是的，关键在于梯度公式里，(\partial f/\partial u) 是逐点计算的。(u^k) 的每个分量对应空间网格上某个点，那个点的敏感性信息已经包含在 (\lambda) 和状态值里了。对任意一个分量 (u_i^k)，梯度就是：

[ \left.\frac{\partial J}{\partial u_i^k}\right|{x_i,t_k} = \Delta t\left[\lambda_c \frac{\partial f}{\partial u} + \lambda_n \frac{\partial g}{\partial u}\right]{x_i,t_k} + 2\gamma u_i^k ]

每个点的计算是独立的，不需要额外解方程。所以伴随梯度的总成本大致是一次正演加一次伴随，也就是两倍正演开销。这个特性对高维优化是决定性的。

3.4 和神经网络反向传播的关系：别把它想得太神秘

我自己最喜欢的一个类比是：伴随方法就是 PDE 界的反向传播。神经网络训练时，前向算损失，反向把损失对每个权重的梯度传回来；这里的正演 PDE 等于前向传播，终端目标函数对应损失函数，伴随方程就是反向传播的递推。区别是神经网络的层数是离散的，PDE 的空间时间尺度是连续的，但推导逻辑一模一样。

想通这一点后，很多直觉可以直接搬过来用。比如梯度消失、梯度爆炸在伴随里也存在——如果时间步长太大或者某些雅可比项的特征值大于 1，反向递推时 (\lambda) 就会被放大，造成数值失稳。解决思路也和深度学习类似：减小步长、加正则、或者换更稳的隐式格式。

4. Matlab实现框架：从连续方程离散化到梯度验证

4.1 空间离散化：二维Laplacian矩阵的构造

我以二维矩形区域为例，网格均匀，(N_x \times N_y) 个节点，网格间距 (dx, dy)。用中心差分构造标准的五点 Laplacian 稀疏矩阵 (A)：

function A = laplacian2D(Nx, Ny, dx, dy) ex = ones(Nx,1); ey = ones(Ny,1); Tx = spdiags([ex -2*ex ex], [-1 0 1], Nx, Nx) / dx^2; Ty = spdiags([ey -2*ey ey], [-1 0 1], Ny, Ny) / dy^2; Ix = speye(Nx); Iy = speye(Ny); A = kron(Iy, Tx) + kron(Ty, Ix); end

这里用kron把二维算子组装起来，顺序是"先把 x 方向变分，再嵌进 y 方向"。边界条件的处理我没写在上面：Neumann 零通量边界下，边界节点上要用 ghost cell 或者直接修改边界行，让边界外侧的值等于内侧值。等我下面强调一个原则：正演用的A和伴随用的A必须是同一个矩阵，最好在代码里只定义一次，严禁在伴随函数里重写一遍近似版本。

4.2 时间离散化：隐式扩散加显式反应

时间上我用一阶隐式欧拉处理扩散，显式处理反应项。每一步要求解：

[ \left(I - \Delta t D_c A\right)c^{k+1} = c^k + \Delta t f(c^k, u^k) ]

写成 Matlab 风格：

Mc = speye(Ng) - p.dt * p.Dc * p.A; [Lc, Uc] = lu(Mc);

这个矩阵是稀疏对称正定的，理论上用chol更好，但lu也够稳。注意，正演过程中Mc不变，只做一次 LU 分解，之后每一步就是两次三角回代，速度非常快。

扩散用隐式格式的好处是稳定性约束从 (D\Delta t/dx^2 < 0.5) 放宽为无约束，时间步可以取大一些。但反应项仍是显式，(\rho) 太大时 (\Delta t) 依然受限。我一般先取 (\Delta t = \min(0.1/\rho, 0.2))，再根据结果调。

4.3 正演求解器长什么样

下面是一个正演求解器的骨架。输入是控制变量矩阵u，尺寸是Ng × Nt，每一列是一个时刻的全网格剂量率。输出是肿瘤密度矩阵C和正常组织密度矩阵Ns，尺寸都是Ng × (Nt+1)。

function [C, Ns] = forward_solve(u, p) Ng = p.Ng; Nt = size(u,2); C = zeros(Ng, Nt+1); Ns = zeros(Ng, Nt+1); C(:,1) = p.c0(:); Ns(:,1) = p.n0(:); Mc = speye(Ng) - p.dt * p.Dc * p.A; Mn = speye(Ng) - p.dt * p.Dn * p.A; [Lc, Uc] = lu(Mc); [Ln, Un] = lu(Mn); for k = 1:Nt uk = u(:,k); fc = p.rho .* C(:,k) .* (1 - C(:,k)/p.K) ... - (p.alpha_c*uk + p.beta_c*uk.^2) .* C(:,k); fn = -(p.alpha_n*uk + p.beta_n*uk.^2) .* Ns(:,k); C(:,k+1) = Uc \ (Lc \ (C(:,k) + p.dt*fc)); Ns(:,k+1) = Un \ (Ln \ (Ns(:,k) + p.dt*fn)); end end

整个计算过程中，所有反应项都用点运算，Matlab 里对Ng×1的向量操作很快，不需要写成 for 循环逐点算。

4.4 伴随求解器：反着跑一遍，顺便攒梯度

伴随求解器和正演是镜像关系。从终端条件出发，反着遍历时间步。注意我这里的积分面积因子dx*dy被统一吸进了权重wc, wn, gamma，这样梯度公式里不出现网格面积，代码更干净，但权重含义要跟着变。

function grad = adjoint_grad(u, C, Ns, p) Ng = p.Ng; Nt = size(u,2); lamC = zeros(Ng, Nt+1); lamN = zeros(Ng, Nt+1); lamC(:,end) = 2*p.wc .* C(:,end); lamN(:,end) = 2*p.wn .* (Ns(:,end) - 1); grad = zeros(size(u)); Mc = speye(Ng) - p.dt * p.Dc * p.A; Mn = speye(Ng) - p.dt * p.Dn * p.A; [Lc, Uc] = lu(Mc); [Ln, Un] = lu(Mn); for k = Nt:-1:1 qc = Uc \ (Lc \ lamC(:,k+1)); qn = Un \ (Ln \ lamN(:,k+1)); uk = u(:,k); dfdu = -(p.alpha_c + 2*p.beta_c*uk) .* C(:,k); dgdu = -(p.alpha_n + 2*p.beta_n*uk) .* Ns(:,k); grad(:,k) = p.dt * (dfdu .* qc + dgdu .* qn) + 2*p.gamma*uk; dfdc = p.rho - 2*p.rho*C(:,k)/p.K - (p.alpha_c*uk + p.beta_c*uk.^2); dgdn = -(p.alpha_n*uk + p.beta_n*uk.^2); lamC(:,k) = qc + p.dt * dfdc .* qc; lamN(:,k) = qn + p.dt * dgdn .* qn; end end

我把梯度计算放在伴随变量更新之前，是因为梯度公式里用的是 (q_c, q_n)，也就是 (M^{-1}\lambda^{k+1})，不是 (\lambda^{k+1})。这个细节非常容易写错，一旦写错，梯度验证过不了，而你不会立刻意识到是这里错。

4.5 梯度验证：用有限差分给伴随梯度作证

无论伴随理论推得多漂亮，代码写完后第一件事永远是梯度验证。随机选几个控制变量的分量，用有限差分公式：

[ g_i^{\text{FD}} = \frac{J(u+\varepsilon e_i) - J(u-\varepsilon e_i)}{2\varepsilon} ]

中心差分比单侧差分精度高一个量级。我一般取 (\varepsilon = 10^{-5}) 到 (10^{-4})。下面是我随便挑的几个测试分量，伴随梯度相对误差大概在 (10^{-4}) 量级，说明离散推导是对的。

如果相对误差大于 (10^{-1})，基本可以确定离散伴随方程写错了，优先检查边界条件和雅可比转置方向。

4.6 优化循环：投影梯度法起步

有了正演、伴随梯度，优化循环其实很简单。初始 (u) 可以取常数剂量率，然后迭代：

for iter = 1:maxIter [C, Ns] = forward_solve(u, p); J = p.wc * sum(C(:,end).^2) ... + p.wn * sum((1-Ns(:,end)).^2) ... + p.gamma * sum(u(:).^2); g = adjoint_grad(u, C, Ns, p); u = u - alpha * g; u = min(max(u, 0), p.u_max); if norm(g, inf) < 1e-6, break; end end

步长 (\alpha) 可以用 Armijo 线搜索，也可以用 L-BFGS 把收敛速度提上来。我这个骨架用的是最朴素的投影梯度法，好处是代码清晰，适合验证整个流程；实际要跑更复杂的临床场景，建议换fmincon或者写一个 L-BFGS 外循环。

5. 优化迭代与避坑记录：那些跑炸的仿真和调参心得

5.1 伴随反向求解的数值失稳

理论上伴随方程是线性的，但它依然是时间反向的 PDE，终端条件如果给得很极端，反向过程中也会振荡。我最开始跑的时候，肿瘤初始条件里人为加了一个局部尖峰，结果显示到终端时刻那个尖峰还没有完全消退，于是 (\lambda_c(T)=2\omega_c c(T)) 在那个位置有个很大的值。反向回传时，这个尖峰被扩散和反应雅可比项放大，伴随解出现了明显的上下跳动，最后梯度里对应区域的值也大得离谱，优化器被带偏。

解决办法有两个：一是把终端附近的 (\Delta t) 加密，让尖峰被数值格式消化掉；二是在目标函数里对终端状态做一点空间光滑，比如在求 (c(T)^2) 之前先做一次平滑滤波。第二个办法有点 hack，但对稳定性很有用。更稳妥的做法是在离散伴随中使用和正演完全一致的隐式格式，让雅可比转置的特征值严格受控。

5.2 边界条件不一致：梯度看着对，实际上全是错的

这是整个实现里最容易踩的坑，也是我花了最多时间才抓出来的 bug。正演在边界上用 Neumann 零通量，我当初写伴随代码时，为了省事，没有直接复用同一个 (A) 矩阵，而是重新写了一个"看起来差不多"的 Laplacian 算子，结果边界行的处理方式和正演不一致。

伴随方程里，空间算子和状态方程必须是同一个离散算子，否则分部积分推导中的边界项不会抵消，梯度公式里会多出一堆边界残差。这个 bug 的可怕之处在于：如果你只检查内部网格点的梯度，误差可能很小；但靠近边界的梯度，相对误差能到 (10^{-1}) 量级。优化器一开始被内部点的大梯度主导，注意不到边界，等内部区域优化得差不多了，边界残差才跳出来捣乱。

解决办法就一句话：在 Matlab 里只写一次A = laplacian2D(...)，正演和伴随都从同一个变量读取。千万别在伴随函数里重新生成了一个"等价"的边界处理。

5.3 量纲与权重：ω_c、ω_n、γ 的拉扯

肿瘤细胞密度和正常组织密度的数值范围不一样。肿瘤 (c) 可能从 0.01 增殖到接近 1，正常组织 (n) 始终在 1 附近，辐射损伤让它掉到 0.8 或者 0.9。如果 (\omega_c) 和 (\omega_n) 随便取，梯度几乎会被肿瘤项主导，正常组织的保护形同虚设。

我这里的一个实际经验是先固定 (\omega_c = 1)，把 (\omega_n) 从 0.1 扫到 10，观察目标函数里两项的变化；然后再调 (\gamma)，观察平均剂量和肿瘤控制之间的关系。(\gamma) 太小，优化器会直接把剂量推到约束上限 (u_{\max})，产生一个"在安全边界上狂飙"的优化解；(\gamma) 太大，剂量被压得过低，肿瘤控制不住。我最后常用的组合是 (\omega_c=1,\ \omega_n=5,\ \gamma=0.1)，但不同模型参数必须重新扫。

5.4 内存与检查点策略

伴随反向需要用到所有时间层的状态 (c^k, n^k)。如果网格 64×64、40 个时间步、两个状态变量，全存下来大约是 64×64×40×2×8 字节，约 2.6 MB，不算大。但如果把网格推到 128×128×50，或者三维 32×32×32×50，这个数字就变成数百 MB 甚至几个 GB，Matlab 开始卡。

缓解办法是检查点策略：每 (m) 步存一次状态，反向求解时遇到检查点就重新从检查点正演到当前步，把中间状态再次算出来。这是典型的时间换内存。Matlab 里实现起来也不复杂，只是代码会多一点。对多数教学和科研原型，全量存储其实能忍，等到要跑大规模三维问题时再考虑检查点。

5.5 一点扩展：把伴随梯度用到参数辨识上

做完了剂量优化，我才意识到这套伴随框架的真正价值不止于控制。固定一个放疗计划 (u)，想要知道模型中哪个参数对治疗结局影响最大——是增殖率还是扩散系数？此时可以把想辨识的参数当作"控制变量"，一样跑伴随梯度。这个梯度可以用于参数敏感性排序、最大似然估计，甚至用随机梯度朗之万动力学做贝叶斯后验采样。

做参数辨识时，控制变量维度低，可能只有三五个参数，用有限差分也能跑，但伴随方法的优势依然存在：它给你的不只是一个梯度，还有每个时空点的敏感性场。举个例子，计算 (\partial J/\partial \rho(x))，可以直接画出"哪个位置的增殖率对治疗失败影响最大"，这在生物学解释上比单个标量参数有说服力得多。如果你已经在跑时空放疗优化，下一步想做鲁棒优化或者在线自适应计划，这套伴随灵敏度框架可以直接把它们串起来。