JFB, SHINE, and SILVA Backward Methods¶
This lab derives exact implicit differentiation, JFB [[88]], and SHINE [[89]], then compares their gradients against an analytic fixed point.
Numbered literature: [1], [4], [10], [11], [13], [87], [88], [89]. Each number opens the complete citation and its primary external source.
In [1]:
Copied!
from pathlib import Path
import sys
root = Path.cwd()
while root != root.parent and not (root / "src" / "silva_networks").exists():
root = root.parent
if not (root / "src" / "silva_networks").exists():
root = Path("/content/silva-networks")
sys.path.insert(0, str(root / "src"))
import matplotlib.pyplot as plt
import torch
plt.rcParams.update({"figure.dpi": 300, "savefig.dpi": 300})
torch.manual_seed(91)
from pathlib import Path
import sys
root = Path.cwd()
while root != root.parent and not (root / "src" / "silva_networks").exists():
root = root.parent
if not (root / "src" / "silva_networks").exists():
root = Path("/content/silva-networks")
sys.path.insert(0, str(root / "src"))
import matplotlib.pyplot as plt
import torch
plt.rcParams.update({"figure.dpi": 300, "savefig.dpi": 300})
torch.manual_seed(91)
Out[1]:
<torch._C.Generator at 0x123f28b90>
1. One Equation, Three Adjoint Choices¶
Exact Implicit Differentiation¶
For $z^\star=T_\theta(z^\star,x)$,
$$ (I-J_T^\top)u=g, \qquad \frac{d\mathcal L}{d\theta}=u^\top\frac{\partial T_\theta}{\partial\theta}. $$
The exact path solves for $u$ and retains the complete local inverse action.
Jacobian-Free Backpropagation¶
JFB uses $u\approx g$, replacing the inverse adjoint factor by identity. This removes the backward linear solve while preserving the converged forward root.
Inverse Estimate Sharing¶
If forward Broyden retains $B\approx(J_T-I)^{-1}$, SHINE starts from
$$ u_0=-B^\top g $$
and may refine the exact adjoint residual.
Compare All Three Gradients¶
The next cell computes the analytic gradient and evaluates exact implicit, JFB, raw SHINE, and refined SHINE routes under one fixed transition.
In [2]:
Copied!
from silva_networks import SolverConfig, solve_equilibrium
def scalar_gradient(mode, refine_steps=0):
bias = torch.nn.Parameter(torch.tensor([0.4], dtype=torch.float64))
weight = torch.nn.Parameter(torch.tensor([0.2], dtype=torch.float64))
config = SolverConfig(
solver="broyden" if mode == "shine" else "anderson",
max_iter=60,
tol=1e-11,
history=8,
anderson_batch_dims=0,
backward_mode=mode,
backward_solver="gmres",
backward_max_iter=8,
backward_tol=1e-11,
shine_refine_steps=refine_steps,
)
result = solve_equilibrium(
lambda z: bias + weight * z,
torch.zeros(1, dtype=torch.float64),
config,
params=(bias, weight),
)
(0.5 * result.z.square().sum()).backward()
return result, float(bias.grad), float(weight.grad)
exact_z = 0.4 / (1.0 - 0.2)
expected_bias = exact_z / (1.0 - 0.2)
expected_weight = exact_z**2 / (1.0 - 0.2)
print("analytic:", exact_z, expected_bias, expected_weight)
for mode, refinement in (("implicit", 0), ("jfb", 0), ("shine", 0), ("shine", 2)):
result, bias_gradient, weight_gradient = scalar_gradient(mode, refinement)
print(
mode,
"refine", refinement,
"state", float(result.z),
"bias grad", bias_gradient,
"weight grad", weight_gradient,
"backward residual", result.info.get("backward_residual"),
)
from silva_networks import SolverConfig, solve_equilibrium
def scalar_gradient(mode, refine_steps=0):
bias = torch.nn.Parameter(torch.tensor([0.4], dtype=torch.float64))
weight = torch.nn.Parameter(torch.tensor([0.2], dtype=torch.float64))
config = SolverConfig(
solver="broyden" if mode == "shine" else "anderson",
max_iter=60,
tol=1e-11,
history=8,
anderson_batch_dims=0,
backward_mode=mode,
backward_solver="gmres",
backward_max_iter=8,
backward_tol=1e-11,
shine_refine_steps=refine_steps,
)
result = solve_equilibrium(
lambda z: bias + weight * z,
torch.zeros(1, dtype=torch.float64),
config,
params=(bias, weight),
)
(0.5 * result.z.square().sum()).backward()
return result, float(bias.grad), float(weight.grad)
exact_z = 0.4 / (1.0 - 0.2)
expected_bias = exact_z / (1.0 - 0.2)
expected_weight = exact_z**2 / (1.0 - 0.2)
print("analytic:", exact_z, expected_bias, expected_weight)
for mode, refinement in (("implicit", 0), ("jfb", 0), ("shine", 0), ("shine", 2)):
result, bias_gradient, weight_gradient = scalar_gradient(mode, refinement)
print(
mode,
"refine", refinement,
"state", float(result.z),
"bias grad", bias_gradient,
"weight grad", weight_gradient,
"backward residual", result.info.get("backward_residual"),
)
analytic: 0.5 0.625 0.3125
implicit refine 0 state 0.49999999999791656 bias grad 0.6249999999973956 weight grad 0.3124999999973957 backward residual 5.551115123125783e-17 jfb refine 0 state 0.49999999999958333 bias grad 0.49999999999958333 weight grad 0.24999999999874994 backward residual None shine refine 0 state 0.5 bias grad 0.625 weight grad 0.3125 backward residual 0.0 shine refine 2 state 0.5 bias grad 0.625 weight grad 0.3125 backward residual 0.0
2. Read the Difference¶
JFB intentionally differs from the exact gradient because it replaces the inverse adjoint factor by identity. SHINE should approach the exact result when the forward inverse estimate is accurate or refinement is sufficient. That is a numerical claim we can test directly.
In [3]:
Copied!
from silva_networks import fixed_point, shine_adjoint_solve
matrix = torch.tensor([[0.35, 0.12], [-0.08, 0.2]], dtype=torch.float64)
source = torch.tensor([0.3, -0.2], dtype=torch.float64)
forward = fixed_point(
lambda z: source + matrix @ z,
torch.zeros(2, dtype=torch.float64),
SolverConfig(solver="broyden", max_iter=12, tol=1e-12, history=6),
)
gradient = torch.tensor([1.0, -0.4], dtype=torch.float64)
exact = torch.linalg.solve(torch.eye(2, dtype=torch.float64) - matrix.T, gradient)
refinement_steps = list(range(6))
errors = []
residuals = []
for steps in refinement_steps:
result = shine_adjoint_solve(
lambda z: source + matrix @ z,
forward.z,
gradient,
forward.inverse_estimate,
refine_steps=steps,
tol=1e-12,
)
errors.append(float(torch.linalg.vector_norm(result.x - exact)))
residuals.append(result.residual)
print(steps, "error", errors[-1], "residual", residuals[-1])
from silva_networks import fixed_point, shine_adjoint_solve
matrix = torch.tensor([[0.35, 0.12], [-0.08, 0.2]], dtype=torch.float64)
source = torch.tensor([0.3, -0.2], dtype=torch.float64)
forward = fixed_point(
lambda z: source + matrix @ z,
torch.zeros(2, dtype=torch.float64),
SolverConfig(solver="broyden", max_iter=12, tol=1e-12, history=6),
)
gradient = torch.tensor([1.0, -0.4], dtype=torch.float64)
exact = torch.linalg.solve(torch.eye(2, dtype=torch.float64) - matrix.T, gradient)
refinement_steps = list(range(6))
errors = []
residuals = []
for steps in refinement_steps:
result = shine_adjoint_solve(
lambda z: source + matrix @ z,
forward.z,
gradient,
forward.inverse_estimate,
refine_steps=steps,
tol=1e-12,
)
errors.append(float(torch.linalg.vector_norm(result.x - exact)))
residuals.append(result.residual)
print(steps, "error", errors[-1], "residual", residuals[-1])
0 error 0.33577213397197764 residual 0.27003417243680083 1 error 5.757124627010134e-15 residual 4.440892098500626e-15 2 error 5.757124627010134e-15 residual 4.440892098500626e-15 3 error 5.757124627010134e-15 residual 4.440892098500626e-15 4 error 5.757124627010134e-15 residual 4.440892098500626e-15 5 error 5.757124627010134e-15 residual 4.440892098500626e-15
In [4]:
Copied!
fig, ax = plt.subplots(figsize=(5.2, 3.3))
ax.plot(refinement_steps, errors, marker="o", label="adjoint error")
ax.plot(refinement_steps, residuals, marker="s", label="linear residual")
ax.set(xlabel="SHINE refinement steps", ylabel="norm", title="shared inverse refinement")
ax.set_yscale("log")
ax.legend()
fig.tight_layout()
plt.show()
fig, ax = plt.subplots(figsize=(5.2, 3.3))
ax.plot(refinement_steps, errors, marker="o", label="adjoint error")
ax.plot(refinement_steps, residuals, marker="s", label="linear residual")
ax.set(xlabel="SHINE refinement steps", ylabel="norm", title="shared inverse refinement")
ax.set_yscale("log")
ax.legend()
fig.tight_layout()
plt.show()
3. Broyden Inverse Factors Are Public¶
The forward solver represents
$$ B_k=-I+\sum_{j=1}^{r}u_jv_j^\top. $$
history bounds $r$. The result can apply $B_k$, $B_k^\top$, or the
fixed-point adjoint approximation $-B_k^\top$ without a dense matrix.
In [5]:
Copied!
estimate = forward.inverse_estimate
probe = torch.tensor([0.25, -0.5], dtype=torch.float64)
print("retained rank:", estimate.rank)
print("B probe:", estimate.apply_residual_inverse(probe))
print("B^T probe:", estimate.apply_residual_inverse_transpose(probe))
print("-(B^T) probe:", estimate.apply_fixed_point_adjoint_inverse(probe))
estimate = forward.inverse_estimate
probe = torch.tensor([0.25, -0.5], dtype=torch.float64)
print("retained rank:", estimate.rank)
print("B probe:", estimate.apply_residual_inverse(probe))
print("B^T probe:", estimate.apply_residual_inverse_transpose(probe))
print("-(B^T) probe:", estimate.apply_fixed_point_adjoint_inverse(probe))
retained rank: 4 B probe: tensor([-0.4700, 0.6055], dtype=torch.float64) B^T probe: tensor([-0.4727, 0.6041], dtype=torch.float64) -(B^T) probe: tensor([ 0.4727, -0.6041], dtype=torch.float64)
4. Selection Guidance¶
| Method | Additional backward solve | Reuses forward information | Typical reason to test it |
|---|---|---|---|
| implicit | yes | no | highest local adjoint accuracy |
| JFB | no | no | lowest backward solver cost |
| SHINE | optional | Broyden inverse factors | forward/backward numerical reuse |
| phantom | short state trajectory | final state | controllable inexact gradient |
| unrolled | no separate solve | complete finite graph | finite-depth reference |
Hold the transition, forward tolerance, optimizer, data order, and random seed fixed when comparing these methods. Report task metric and gradient agreement together with runtime and memory.
From 49 Jfb Shine Backward Methods to a Custom SILVA Family¶
The construction in this notebook can be separated into the universal conditioned-equilibrium contract
$$ z_0=I_\eta(x),\qquad z^\star=T_\theta(z^\star,x),\qquad \widehat y=Q_\psi(z^\star). $$
For this topic:
| Part | Concrete interpretation |
|---|---|
| Equilibrium state | the converged latent state z_star |
| Condition | differentiable parameters and external inputs |
| Repeated computation | the map whose Jacobian defines I - J_z f |
| Required invariants | agreement of dense, JVP, VJP, and matrix-free products |
| Replaceable components | transition, Jacobian product, linear solver, regularizer, and loss |
The initializer and source path are evaluated outside or alongside the root solve. Only the state-preserving transition is repeated. Replacing an internal architecture does not change this equation, provided the transition still maps the same state space into itself.
In [6]:
Copied!
import torch as silva_extension_torch
from torch import nn as silva_extension_nn
from silva_networks import (
SILVAConditionedEquilibrium,
SILVAZeroInitializer,
SolverConfig,
validate_silva_transition,
)
class NotebookExtensionTransition(silva_extension_nn.Module):
def __init__(self, condition_dim=2, state_dim=3):
super().__init__()
self.source = silva_extension_nn.Linear(condition_dim, state_dim)
self.state_field = silva_extension_nn.Sequential(
silva_extension_nn.Linear(state_dim, 2 * state_dim),
silva_extension_nn.Tanh(),
silva_extension_nn.Linear(2 * state_dim, state_dim),
)
def forward(self, state, condition):
return silva_extension_torch.tanh(
self.source(condition) + 0.15 * self.state_field(state)
)
silva_extension_torch.manual_seed(610)
notebook_condition = silva_extension_torch.linspace(-1.0, 1.0, 8).reshape(4, 2)
notebook_state0 = silva_extension_torch.zeros(4, 3)
notebook_transition = NotebookExtensionTransition()
notebook_report = validate_silva_transition(
notebook_transition,
notebook_state0,
notebook_condition,
)
assert notebook_report.valid
with silva_extension_torch.no_grad():
notebook_reference_step = silva_extension_torch.tanh(
notebook_transition.source(notebook_condition)
+ 0.15 * notebook_transition.state_field(notebook_state0)
)
silva_extension_torch.testing.assert_close(
notebook_transition(notebook_state0, notebook_condition),
notebook_reference_step,
)
notebook_custom_model = SILVAConditionedEquilibrium(
notebook_transition,
SILVAZeroInitializer(3),
readout=silva_extension_nn.Linear(3, 1),
config=SolverConfig(
solver="picard",
max_iter=40,
tol=1e-7,
backward_mode="implicit",
backward_solver="gmres",
anderson_batch_dims=1,
),
)
notebook_custom_result = notebook_custom_model(
notebook_condition,
return_result=True,
)
assert notebook_custom_result.output.shape == (4, 1)
assert notebook_custom_result.solver_result.residual < 1e-5
notebook_custom_result.output.square().mean().backward()
assert all(
parameter.grad is not None and silva_extension_torch.isfinite(parameter.grad).all()
for parameter in notebook_custom_model.parameters()
)
print("custom transition:", notebook_report)
print("equilibrium residual:", notebook_custom_result.solver_result.residual)
import torch as silva_extension_torch
from torch import nn as silva_extension_nn
from silva_networks import (
SILVAConditionedEquilibrium,
SILVAZeroInitializer,
SolverConfig,
validate_silva_transition,
)
class NotebookExtensionTransition(silva_extension_nn.Module):
def __init__(self, condition_dim=2, state_dim=3):
super().__init__()
self.source = silva_extension_nn.Linear(condition_dim, state_dim)
self.state_field = silva_extension_nn.Sequential(
silva_extension_nn.Linear(state_dim, 2 * state_dim),
silva_extension_nn.Tanh(),
silva_extension_nn.Linear(2 * state_dim, state_dim),
)
def forward(self, state, condition):
return silva_extension_torch.tanh(
self.source(condition) + 0.15 * self.state_field(state)
)
silva_extension_torch.manual_seed(610)
notebook_condition = silva_extension_torch.linspace(-1.0, 1.0, 8).reshape(4, 2)
notebook_state0 = silva_extension_torch.zeros(4, 3)
notebook_transition = NotebookExtensionTransition()
notebook_report = validate_silva_transition(
notebook_transition,
notebook_state0,
notebook_condition,
)
assert notebook_report.valid
with silva_extension_torch.no_grad():
notebook_reference_step = silva_extension_torch.tanh(
notebook_transition.source(notebook_condition)
+ 0.15 * notebook_transition.state_field(notebook_state0)
)
silva_extension_torch.testing.assert_close(
notebook_transition(notebook_state0, notebook_condition),
notebook_reference_step,
)
notebook_custom_model = SILVAConditionedEquilibrium(
notebook_transition,
SILVAZeroInitializer(3),
readout=silva_extension_nn.Linear(3, 1),
config=SolverConfig(
solver="picard",
max_iter=40,
tol=1e-7,
backward_mode="implicit",
backward_solver="gmres",
anderson_batch_dims=1,
),
)
notebook_custom_result = notebook_custom_model(
notebook_condition,
return_result=True,
)
assert notebook_custom_result.output.shape == (4, 1)
assert notebook_custom_result.solver_result.residual < 1e-5
notebook_custom_result.output.square().mean().backward()
assert all(
parameter.grad is not None and silva_extension_torch.isfinite(parameter.grad).all()
for parameter in notebook_custom_model.parameters()
)
print("custom transition:", notebook_report)
print("equilibrium residual:", notebook_custom_result.solver_result.residual)
custom transition: SILVATransitionReport(state_shape=(4, 3), output_shape=(4, 3), preserves_shape=True, preserves_device=True, preserves_dtype=True, finite=True, differentiable=True, parameter_count=54) equilibrium residual: 5.960464477539063e-08
Numerical Equivalence, Compact Reproduction, and Scale¶
Before training, compare one packaged transition with an independently written update:
$$ e_{\mathrm{step}} =\frac{\|T_\theta(z,x)-T_{\mathrm{ref}}(z,x)\|_2} {\|T_{\mathrm{ref}}(z,x)\|_2+\varepsilon}. $$
After solving, report the fixed-point residual separately:
$$ e_{\mathrm{fp}} =\frac{\|T_\theta(z^\star,x)-z^\star\|_2} {\|z^\star\|_2+\varepsilon}. $$
For this notebook, a compact reproduction must declare and assert adjoint residual and gradient error against explicit differentiation. A full experiment must additionally record the source dataset version and split, preprocessing, architecture widths, solver and optimizer schedules, random seeds, baseline configuration, checkpoints, and every deviation from the cited protocol.
The principal scaling axes are state dimension and matrix-free linear-solver iterations. Increase one axis at a time, retain the compact deterministic case as a regression test, and record task error, domain-specific residual, forward residual, backward linear residual, memory use, and runtime independently.
Extension Exercises¶
- Replace one component from this notebook while preserving its state and domain invariants.
- Write the replacement first as an independent reference function, then as a module, and assert one-step equivalence.
- Compare two solver configurations on the identical trained transition.
- Add a compact baseline and a predeclared metric threshold.
- Create a full-scale configuration without weakening the compact tests.
The complete authoring protocol is documented in Extending SILVA.
In [7]:
Copied!
notebook_reproduction_record = {
"notebook": '49_jfb_shine_backward_methods.ipynb',
"state": 'the converged latent state z_star',
"condition": 'differentiable parameters and external inputs',
"transition": 'the map whose Jacobian defines I - J_z f',
"invariants": 'agreement of dense, JVP, VJP, and matrix-free products',
"compact_metric": 'adjoint residual and gradient error against explicit differentiation',
"scale_axis": 'state dimension and matrix-free linear-solver iterations',
}
assert all(notebook_reproduction_record.values())
notebook_reproduction_record
notebook_reproduction_record = {
"notebook": '49_jfb_shine_backward_methods.ipynb',
"state": 'the converged latent state z_star',
"condition": 'differentiable parameters and external inputs',
"transition": 'the map whose Jacobian defines I - J_z f',
"invariants": 'agreement of dense, JVP, VJP, and matrix-free products',
"compact_metric": 'adjoint residual and gradient error against explicit differentiation',
"scale_axis": 'state dimension and matrix-free linear-solver iterations',
}
assert all(notebook_reproduction_record.values())
notebook_reproduction_record
Out[7]:
{'notebook': '49_jfb_shine_backward_methods.ipynb',
'state': 'the converged latent state z_star',
'condition': 'differentiable parameters and external inputs',
'transition': 'the map whose Jacobian defines I - J_z f',
'invariants': 'agreement of dense, JVP, VJP, and matrix-free products',
'compact_metric': 'adjoint residual and gradient error against explicit differentiation',
'scale_axis': 'state dimension and matrix-free linear-solver iterations'}
Worked Convergence and Sensitivity Study¶
The preceding example demonstrates one configured solve. This additional study changes the transition feedback factor while keeping the source fixed, so solver effort and implicit sensitivity can be read separately from task behavior. Locally, one eigendirection of a nonlinear transition can be represented by
$$ z_{k+1} = \rho z_k + u, \qquad 0 \leq \rho < 1. $$
Its equilibrium is
$$ z^\star = \frac{u}{1-\rho}. $$
Subtracting the fixed-point equation from the iteration gives the exact error recursion
$$ e_{k+1} = \rho e_k, \qquad |e_k| = \rho^k |e_0|. $$
For a requested absolute tolerance $\tau$, the idealized iteration estimate is
$$ k \geq \frac{\log(\tau/|e_0|)}{\log \rho}. $$
The same factor controls sensitivity. Differentiating the equilibrium with respect to the source gives
$$ \frac{\partial z^\star}{\partial u} =\frac{1}{1-\rho}. $$
Thus a transition can remain contractive while becoming expensive and highly sensitive as $\rho$ approaches one. The table and figure below measure this effect rather than merely stating it. They provide a reference envelope for the notebook's actual state, the converged latent state z_star, and its repeated map, the map whose Jacobian defines I - J_z f. The scalar study does not replace the domain model; it supplies a result whose convergence rate and derivative are known exactly, so the same reporting code can be trusted before it is applied to the larger transition.
In [8]:
Copied!
import math as silva_deepening_math
import torch as silva_deepening_torch
silva_deepening_rates = (0.20, 0.45, 0.70, 0.85)
silva_deepening_source = 0.35
silva_deepening_tolerance = 1e-8
silva_deepening_histories = {}
silva_deepening_rows = []
for silva_deepening_rho in silva_deepening_rates:
silva_deepening_state = silva_deepening_torch.tensor(0.0)
silva_deepening_exact = silva_deepening_source / (1.0 - silva_deepening_rho)
silva_deepening_history = []
for silva_deepening_iteration in range(1, 241):
silva_deepening_next = (
silva_deepening_rho * silva_deepening_state + silva_deepening_source
)
silva_deepening_residual = abs(
float(silva_deepening_next - silva_deepening_state)
)
silva_deepening_history.append(silva_deepening_residual)
silva_deepening_state = silva_deepening_next
if silva_deepening_residual < silva_deepening_tolerance:
break
silva_deepening_u = silva_deepening_torch.tensor(
silva_deepening_source, requires_grad=True
)
silva_deepening_solution = silva_deepening_u / (1.0 - silva_deepening_rho)
silva_deepening_solution.backward()
silva_deepening_expected_sensitivity = 1.0 / (1.0 - silva_deepening_rho)
silva_deepening_gradient_error = abs(
float(silva_deepening_u.grad) - silva_deepening_expected_sensitivity
)
silva_deepening_histories[silva_deepening_rho] = silva_deepening_history
silva_deepening_rows.append(
(
silva_deepening_rho,
silva_deepening_iteration,
silva_deepening_history[-1],
abs(float(silva_deepening_state) - silva_deepening_exact),
float(silva_deepening_u.grad),
silva_deepening_gradient_error,
)
)
print('transition feedback factor')
print("rho | iterations | final residual | exact-state error | sensitivity | gradient error")
for silva_deepening_row in silva_deepening_rows:
print(
f"{silva_deepening_row[0]:.2f} | {silva_deepening_row[1]:3d} | "
f"{silva_deepening_row[2]:.3e} | {silva_deepening_row[3]:.3e} | "
f"{silva_deepening_row[4]:.4f} | {silva_deepening_row[5]:.3e}"
)
assert all(row[2] < silva_deepening_tolerance for row in silva_deepening_rows)
assert all(row[3] < 1e-6 for row in silva_deepening_rows)
assert all(row[5] < 1e-6 for row in silva_deepening_rows)
import math as silva_deepening_math
import torch as silva_deepening_torch
silva_deepening_rates = (0.20, 0.45, 0.70, 0.85)
silva_deepening_source = 0.35
silva_deepening_tolerance = 1e-8
silva_deepening_histories = {}
silva_deepening_rows = []
for silva_deepening_rho in silva_deepening_rates:
silva_deepening_state = silva_deepening_torch.tensor(0.0)
silva_deepening_exact = silva_deepening_source / (1.0 - silva_deepening_rho)
silva_deepening_history = []
for silva_deepening_iteration in range(1, 241):
silva_deepening_next = (
silva_deepening_rho * silva_deepening_state + silva_deepening_source
)
silva_deepening_residual = abs(
float(silva_deepening_next - silva_deepening_state)
)
silva_deepening_history.append(silva_deepening_residual)
silva_deepening_state = silva_deepening_next
if silva_deepening_residual < silva_deepening_tolerance:
break
silva_deepening_u = silva_deepening_torch.tensor(
silva_deepening_source, requires_grad=True
)
silva_deepening_solution = silva_deepening_u / (1.0 - silva_deepening_rho)
silva_deepening_solution.backward()
silva_deepening_expected_sensitivity = 1.0 / (1.0 - silva_deepening_rho)
silva_deepening_gradient_error = abs(
float(silva_deepening_u.grad) - silva_deepening_expected_sensitivity
)
silva_deepening_histories[silva_deepening_rho] = silva_deepening_history
silva_deepening_rows.append(
(
silva_deepening_rho,
silva_deepening_iteration,
silva_deepening_history[-1],
abs(float(silva_deepening_state) - silva_deepening_exact),
float(silva_deepening_u.grad),
silva_deepening_gradient_error,
)
)
print('transition feedback factor')
print("rho | iterations | final residual | exact-state error | sensitivity | gradient error")
for silva_deepening_row in silva_deepening_rows:
print(
f"{silva_deepening_row[0]:.2f} | {silva_deepening_row[1]:3d} | "
f"{silva_deepening_row[2]:.3e} | {silva_deepening_row[3]:.3e} | "
f"{silva_deepening_row[4]:.4f} | {silva_deepening_row[5]:.3e}"
)
assert all(row[2] < silva_deepening_tolerance for row in silva_deepening_rows)
assert all(row[3] < 1e-6 for row in silva_deepening_rows)
assert all(row[5] < 1e-6 for row in silva_deepening_rows)
transition feedback factor rho | iterations | final residual | exact-state error | sensitivity | gradient error 0.20 | 12 | 0.000e+00 | 5.551e-17 | 1.2500 | 0.000e+00 0.45 | 23 | 0.000e+00 | 1.084e-08 | 1.8182 | 6.502e-08 0.70 | 45 | 0.000e+00 | 1.589e-07 | 3.3333 | 7.947e-08 0.85 | 93 | 0.000e+00 | 5.563e-07 | 6.6667 | 1.589e-07
In [9]:
Copied!
import matplotlib.pyplot as silva_deepening_plt
silva_deepening_plt.rcParams.update({"figure.dpi": 300, "savefig.dpi": 300})
silva_deepening_figure, silva_deepening_axes = silva_deepening_plt.subplots(
1, 2, figsize=(8.6, 3.2)
)
for silva_deepening_rho, silva_deepening_history in silva_deepening_histories.items():
silva_deepening_axes[0].semilogy(
range(1, len(silva_deepening_history) + 1),
silva_deepening_history,
marker="o",
markersize=2,
linewidth=1.2,
label=f"rho={silva_deepening_rho:.2f}",
)
silva_deepening_axes[0].axhline(
silva_deepening_tolerance, color="black", linestyle="--", linewidth=0.9
)
silva_deepening_axes[0].set_xlabel("iteration")
silva_deepening_axes[0].set_ylabel("absolute residual")
silva_deepening_axes[0].set_title("Residual trajectories")
silva_deepening_axes[0].legend(fontsize=7)
silva_deepening_axes[1].plot(
[row[0] for row in silva_deepening_rows],
[row[1] for row in silva_deepening_rows],
marker="o",
label="iterations",
)
silva_deepening_sensitivity_axis = silva_deepening_axes[1].twinx()
silva_deepening_sensitivity_axis.plot(
[row[0] for row in silva_deepening_rows],
[row[4] for row in silva_deepening_rows],
color="tab:red",
marker="s",
label="sensitivity",
)
silva_deepening_axes[1].set_xlabel('transition feedback factor')
silva_deepening_axes[1].set_ylabel("iterations")
silva_deepening_sensitivity_axis.set_ylabel("implicit sensitivity", color="tab:red")
silva_deepening_axes[1].set_title("Cost and sensitivity")
silva_deepening_figure.tight_layout()
silva_deepening_plt.show()
import matplotlib.pyplot as silva_deepening_plt
silva_deepening_plt.rcParams.update({"figure.dpi": 300, "savefig.dpi": 300})
silva_deepening_figure, silva_deepening_axes = silva_deepening_plt.subplots(
1, 2, figsize=(8.6, 3.2)
)
for silva_deepening_rho, silva_deepening_history in silva_deepening_histories.items():
silva_deepening_axes[0].semilogy(
range(1, len(silva_deepening_history) + 1),
silva_deepening_history,
marker="o",
markersize=2,
linewidth=1.2,
label=f"rho={silva_deepening_rho:.2f}",
)
silva_deepening_axes[0].axhline(
silva_deepening_tolerance, color="black", linestyle="--", linewidth=0.9
)
silva_deepening_axes[0].set_xlabel("iteration")
silva_deepening_axes[0].set_ylabel("absolute residual")
silva_deepening_axes[0].set_title("Residual trajectories")
silva_deepening_axes[0].legend(fontsize=7)
silva_deepening_axes[1].plot(
[row[0] for row in silva_deepening_rows],
[row[1] for row in silva_deepening_rows],
marker="o",
label="iterations",
)
silva_deepening_sensitivity_axis = silva_deepening_axes[1].twinx()
silva_deepening_sensitivity_axis.plot(
[row[0] for row in silva_deepening_rows],
[row[4] for row in silva_deepening_rows],
color="tab:red",
marker="s",
label="sensitivity",
)
silva_deepening_axes[1].set_xlabel('transition feedback factor')
silva_deepening_axes[1].set_ylabel("iterations")
silva_deepening_sensitivity_axis.set_ylabel("implicit sensitivity", color="tab:red")
silva_deepening_axes[1].set_title("Cost and sensitivity")
silva_deepening_figure.tight_layout()
silva_deepening_plt.show()
Reading and Extending the Result¶
The measured residual curves become flatter as the transition feedback factor increases. The iteration count and the exact sensitivity rise together, but they answer different questions: iterations measure numerical work, while sensitivity describes how strongly the equilibrium reacts to the source. The gradient-error column verifies the differentiation path against the analytic derivative.
Apply the same separation to this notebook's full model:
| Report | Notebook-specific interpretation |
|---|---|
| Task evidence | adjoint residual and gradient error against explicit differentiation |
| Forward residual | Re-evaluate the complete transition at the returned state |
| Empirical rate | Compare consecutive residuals only after the transient regime |
| Backward residual | Record the linear-adjoint stopping value independently |
| Sensitivity | Perturb one declared source field while preserving all other inputs |
| Structural checks | agreement of dense, JVP, VJP, and matrix-free products |
| Scale sweep | Change one of state dimension and matrix-free linear-solver iterations at a time |
A richer experiment should now repeat the sweep with at least two forward solvers, two tolerances, and multiple seeds. Keep model parameters and data identical when comparing solvers. Then change one architecture or data-scale axis, retain the compact analytic study as a regression test, and report task quality, residuals, iterations, runtime, memory, gradient norms, and failed convergence cases together.