Full-Scale SILVA: Derivations, Equivalence Checks, and Training¶
This lab connects every canonical SILVA family to one execution contract. It derives the memory-aware forms used for attention, monotone graph maps, empirical measures, physics-informed derivatives, and implicit DAE stages. It then writes lazy data shards, trains a Fourier equilibrium, resumes from a checkpoint, and records the diagnostics required before a larger study.
The small tensors make the numerical identities inspectable. They are mechanism checks, not substitutes for the official benchmark splits and metrics cited for each family [4, 5, 31, 36, 43-52].
Numbered literature: [1], [4], [5], [13], [27], [29], [31], [32], [36], [43], [45], [51], [52], [54], [55], [56], [57], [58], [73], [74]. Each number opens the complete citation and its primary external source.
from pathlib import Path
import importlib.util
import subprocess
import sys
REPO_URL = "https://github.com/jseluis/silva-networks.git"
def find_local_silva_root():
candidates = [Path.cwd(), Path("/content/silva-networks")]
root = Path.cwd()
while root != root.parent:
candidates.append(root)
root = root.parent
for candidate in candidates:
if (candidate / "src" / "silva_networks").exists():
return candidate
return None
root = find_local_silva_root()
if root is not None:
sys.path.insert(0, str(root / "src"))
elif importlib.util.find_spec("silva_networks") is None:
subprocess.check_call([sys.executable, "-m", "pip", "install", f"git+{REPO_URL}"])
root = Path.cwd()
else:
root = Path.cwd()
import tempfile
import matplotlib.pyplot as plt
import torch
from silva_networks import (
SILVAImplicitDAEStep,
SILVAInjectedSelfAttention,
SILVAMonotoneGraphTransition,
SILVAShardedTensorDataset,
all_silva_family_guides,
audit_silva_family_guides,
build_scaled_silva,
distributional_discrepancy,
fit_supervised,
full_scale_solver_config,
make_periodic_elliptic_dataset,
make_silva_dataloader,
runtime_for_tier,
write_silva_tensor_shards,
)
torch.manual_seed(260)
plt.rcParams.update({"figure.dpi": 300, "savefig.dpi": 300})
1. The Contract That Every Family Preserves¶
A SILVA point solves
$$ \begin{aligned} z^\star&=F_\theta(z^\star;x),\\ z^\star&=\Phi\!\left(S_\theta(x)+H_\theta(z^\star)\right.\\ &\qquad\left.{}+L_\theta(z^\star)+G_\theta(z^\star)\right). \end{aligned} $$
The internal state may be a vector, sequence, image pyramid, graph field, Fourier field, diffusion trajectory, or empirical measure. The defining shape contract is
$$ F_\theta:\mathcal Z\times\mathcal X\longrightarrow\mathcal Z. $$
A U-Net [27], Fourier operator [31], attention block [29], graph operator [36], or user module may live inside one point when its output is projected back to the declared state space.
2. Forward and Backward Equations¶
The forward residual is
$$ r_k=\frac{\lVert F_\theta(z_k;x)-z_k\rVert} {\varepsilon+\lVert F_\theta(z_k;x)\rVert}. $$
At equilibrium, implicit differentiation solves
$$ \begin{aligned} \left(I-J_F^\top\right)u &=\frac{\partial\mathcal L}{\partial z^\star},\\ \frac{\partial\mathcal L}{\partial\theta} &=u^\top\frac{\partial F_\theta}{\partial\theta}. \end{aligned} $$
The implementation sends Jacobian-vector or vector-Jacobian products to an iterative linear solver; full Jacobians are retained only as optional low-dimensional teaching paths [13].
smoke_solver = full_scale_solver_config(tier="smoke")
full_solver = full_scale_solver_config(tier="full")
assert smoke_solver.backward_mode == full_solver.backward_mode == "implicit"
assert smoke_solver.stop_mode == full_solver.stop_mode == "relative"
print("forward budgets:", smoke_solver.max_iter, full_solver.max_iter)
print("backward budgets:", smoke_solver.backward_max_iter, full_solver.backward_max_iter)
forward budgets: 12 60 backward budgets: 20 80
3. Every Canonical Family Has an Actionable Route¶
The family guide is executable metadata used by the public factory and command line. Each row declares a tensor/data contract, primary literature, a benchmark route, scale controls, and extension points. Aliases resolve to these canonical SILVA names.
guides = all_silva_family_guides()
assert audit_silva_family_guides() == ()
assert len(guides) == 64
for index, guide in enumerate(guides, start=1):
refs = ", ".join(f"[{number}]" for number in guide.paper_refs)
print(f"{index:2d}. {guide.family:39s} {refs:18s} {guide.benchmark_tasks[0]}")
1. silva_layer [1], [4] task-defined tensor or graph prediction 2. silva_graph [1], [15], [16] node classification 3. silva_graph_preset [1], [16], [17] citation networks 4. silva_cortex [1], [4] task-defined point, image, field, or graph mapping 5. silva_cortex_network [1], [4] multistage or multimodal tasks 6. silva_image_cortex [1], [27], [29] CIFAR-10 7. compact_deq [4] WikiText-103-style sequence cores 8. message_passing_deq [4], [16] node classification 9. mdeq [5] CIFAR-10 teaching bridge 10. multiscale_vision_deq [5] ImageNet classification 11. sequence_deq [4] WikiText-103 language modeling 12. implicit_graph [36] chain tasks 13. implicit_neural_representation [37] image fitting 14. diffusion_equilibrium [38], [49] CIFAR-10 diffusion 15. scientific_operator [31], [32] Darcy flow 16. fourier_operator_equilibrium [31], [32] Darcy flow 17. implicit_time_step [7] stiff ODE steps 18. silva_deq_flow [22], [23] FlyingChairs 19. raft_deq_flow [22], [23] FlyingChairs 20. quadratic_optimization [8] synthetic QPs 21. silva_projected_qp [8], [9] b 22. silva_fno_deq [31], [43] Darcy flow 23. silva_physics_graph_deq [44] convection-diffusion sensor graphs 24. silva_homotopy_equilibrium [7], [46], [58] CIFAR-10/100 image classification 25. silva_distributional_deq [45] ModelNet40 classification 26. silva_monotone_graph_equilibrium [47] node classification 27. silva_generative_equilibrium_transformer [48] offline CIFAR-10 diffusion distillation 28. silva_poisson_mirror_equilibrium [50] Poisson inverse imaging with a declared forward operator 29. silva_physics_informed_equilibrium [51] Van der Pol oscillator 30. silva_implicit_dae_step [52] three-bus power-network DAE 31. silva_consistency_deq [59] WikiText-103 32. silva_psi_gnn [60] synthetic unstructured Poisson meshes 33. silva_ifno [61] Darcy flow 34. silva_snarf [62] 2D Stick 35. silva_mesh_inference [63] synthetic carrier chains 36. silva_physics_guided_diffusion_pde [64] Poisson 37. silva_therino [73] periodic heterogeneous elastic localization 38. silva_fixed_point_diffusion [74] ImageNet 39. silva_monotone_operator_equilibrium [75] MNIST 40. silva_positive_concave_equilibrium [76] MNIST 41. silva_non_euclidean_equilibrium [77] MNIST 42. silva_efficient_infinite_graph [78] Cora 43. silva_multiscale_graph_implicit [79] Cora 44. silva_delta_equilibrium [80] FlyingChairs 45. silva_hyper_deq [87] WikiText-103 46. silva_quantum_deq [90] MNIST-4 47. silva_bayesian_deq [94] MNIST 48. silva_joint_inference_equilibrium [95] image denoising 49. silva_implicit_spatiotemporal [96] advection-diffusion 50. silva_certified_equilibrium [97], [98] MNIST 51. silva_lipschitz_mdeq [99] CIFAR-10 52. silva_subhomogeneous_equilibrium [100] MNIST 53. silva_algorithmic_reasoner [101] CLRS-30 54. silva_hamiltonian_equilibrium [102] MD17 55. silva_inverse_imaging_equilibrium [103] compressive sensing 56. silva_snapshot_compressive_equilibrium [104] Kobe 57. silva_magnetic_particle_equilibrium [105] OpenMPIData 58. silva_sparse_hyperspectral_equilibrium [106] ICVL 59. silva_serialized_smoothing_equilibrium [107] CIFAR-10 60. silva_diffusion_restoration_equilibrium [108] ImageNet 61. silva_recurrent_equilibrium_network [109] system identification 62. silva_lipschitz_robust_equilibrium [110] MNIST 63. silva_image_matting_equilibrium [111] Adobe Composition-1k 64. silva_dynamic_economic_equilibrium [112] stochastic growth
4. Runtime Tiers Do Not Change the Mathematics¶
smoke uses small solver budgets and ordinary precision.
workstation keeps neutral precision and leaves distribution off.
full selects larger solver budgets, gradient accumulation,
bfloat16, and distributed loading. Every value can be overridden.
For $P$ processes and $K$ accumulated microbatches,
$$ B_{\mathrm{effective}}=B_{\mathrm{device}}KP. $$
for tier in ("smoke", "workstation", "full"):
runtime = runtime_for_tier(tier)
print(
tier,
"precision=", runtime.mixed_precision,
"distributed=", runtime.distributed,
"effective batch on 4 processes=", runtime.effective_batch_size(world_size=4),
)
smoke precision= none distributed= False effective batch on 4 processes= 16 workstation precision= none distributed= False effective batch on 4 processes= 32 full precision= bfloat16 distributed= True effective batch on 4 processes= 128
5. Factorized Monotone Graph Operator¶
The constrained graph channel map is
$$ \begin{aligned} D&=(1-m)I-CC^\top,\\ S&=UV^\top-VU^\top,\\ W&=D+S. \end{aligned} $$
Each factor $C$, $U$, and $V$ has shape $d\times r$.
Its symmetric certificate follows directly:
$$ \begin{aligned} W_{\mathrm{sym}}&=\frac{W+W^\top}{2},\\ I-W_{\mathrm{sym}}&=mI+CC^\top,\\ I-W_{\mathrm{sym}}&\succeq mI. \end{aligned} $$
With rank $r$, applying the factors costs $O(Ndr)$ storage-aware arithmetic instead of constructing a dense $d\times d$ matrix. The next cell verifies that both forms are numerically identical.
graph_transition = SILVAMonotoneGraphTransition(
in_dim=3,
state_dim=16,
operator_rank=4,
margin=0.15,
)
node_values = torch.randn(11, 16)
factorized = graph_transition.apply_channel_weight(node_values)
explicit = node_values @ graph_transition.channel_weight().T
graph_error = (factorized - explicit).abs().max()
assert graph_error < 1e-5
assert graph_transition.monotonicity_certificate() >= 0.15 - 1e-6
print("factorized/dense maximum error:", float(graph_error))
print("monotonicity certificate:", float(graph_transition.monotonicity_certificate()))
factorized/dense maximum error: 4.76837158203125e-07 monotonicity certificate: 0.1499999314546585
6. Injected Attention Without a Required Score Matrix¶
Generative equilibrium attention injects source information once:
$$ \begin{aligned} Q&=ZW_q+U_q,\\ K&=ZW_k+U_k,\\ V&=ZW_v+U_v,\\ A&=\operatorname{softmax}\!\left(\frac{QK^\top}{\sqrt{d_h}}\right)\\ &\qquad{}\cdot V. \end{aligned} $$
The manual path exposes the equation. Fused attention uses the backend scaled-dot-product kernel, while query chunking bounds the explicit query workspace. All three must compute the same map within floating-point tolerance.
manual_attention = SILVAInjectedSelfAttention(16, heads=4, attention_mode="manual")
fused_attention = SILVAInjectedSelfAttention(16, heads=4, attention_mode="sdpa")
chunked_attention = SILVAInjectedSelfAttention(
16,
heads=4,
attention_mode="chunked",
query_chunk_size=5,
)
fused_attention.load_state_dict(manual_attention.state_dict())
chunked_attention.load_state_dict(manual_attention.state_dict())
state = torch.randn(2, 13, 16)
injection = torch.randn(2, 13, 48)
attention_outputs = [
module(state, injection)
for module in (manual_attention, fused_attention, chunked_attention)
]
attention_errors = [
float((output - attention_outputs[0]).abs().max())
for output in attention_outputs[1:]
]
assert max(attention_errors) < 2e-5
print("fused and chunked errors:", attention_errors)
fused and chunked errors: [1.1920928955078125e-07, 1.7881393432617188e-07]
7. Distributional Equilibria and Exact Pair Chunking¶
For empirical measures $\mu_Z$ and $\mu_X$, a distributional SILVA point may minimize
$$ \mathcal E(Z)=\frac12D^2\!\left(\mu_Z, \mu_{F_\theta(Z,X)}\right). $$
Gaussian MMD and energy distance contain all particle pairs. A chunked reduction preserves the exact sum and gradient while reducing peak pair storage from $O(NM)$ to $O(CM)$ [45].
left = torch.randn(3, 17, 4, requires_grad=True)
right = torch.randn(3, 19, 4)
dense_measure = distributional_discrepancy(left, right, kernel="gaussian")
chunked_measure = distributional_discrepancy(
left,
right,
kernel="gaussian",
pairwise_chunk_size=5,
)
measure_error = (dense_measure - chunked_measure).abs()
chunked_measure.backward()
assert measure_error < 1e-6
assert left.grad is not None and torch.isfinite(left.grad).all()
print("dense/chunked discrepancy error:", float(measure_error))
dense/chunked discrepancy error:
2.2351741790771484e-08
8. Physics-Informed Equilibrium Derivative¶
If $z^\star=f_\theta(z^\star,t)$, differentiation gives
$$ \left(I-J_zf_\theta\right)\frac{dz^\star}{dt}=J_tf_\theta. $$
The matrix-free operator is
$$ v\longmapsto v-J_zf_\theta v. $$
It is evaluated with JVPs and solved by GMRES. The dense solve is retained for small-state verification and produces the comparison below [51].
physics_model = build_scaled_silva(
"pideq",
tier="smoke",
state_dim=5,
output_dim=2,
derivative_mode="matrix_free",
derivative_max_iter=30,
)
times = torch.linspace(0.0, 1.0, 4).unsqueeze(-1)
physics_output = physics_model(times, return_result=True)
dense_derivative = physics_model.implicit_time_derivative(
times,
physics_output.state,
mode="dense",
)
matrix_free_derivative = physics_model.implicit_time_derivative(
times,
physics_output.state,
mode="matrix_free",
)
derivative_error = (dense_derivative - matrix_free_derivative).abs().max()
assert derivative_error < 2e-4
print("dense/matrix-free derivative error:", float(derivative_error))
print("equilibrium residual:", physics_output.solver_result.residual)
dense/matrix-free derivative error: 1.341104507446289e-07 equilibrium residual: 5.8430846365808975e-06
9. Implicit DAE Stage as Newton-Krylov SILVA¶
For $y'=f(y,a)$ and $0=g(y,a)$, an implicit Runge-Kutta step solves
$$ \begin{aligned} s_j&=\sum_i a_{ji}f(Y_i,A_i),\\ Y_j&=y_n+h s_j,\\ g(Y_j,A_j)&=0. \end{aligned} $$
Stacking stage and endpoint equations gives $R(q)=0$. Newton's correction satisfies
$$ \begin{aligned} J_R(q_k)\delta_k&=R(q_k),\\ q_{k+1}&=q_k-\lambda\delta_k. \end{aligned} $$
The Krylov path supplies $v\mapsto J_R(q_k)v+\rho v$ through a JVP, avoiding the dense stage Jacobian [52].
dense_dae = SILVAImplicitDAEStep(
max_iter=5,
tol=1e-7,
linear_solver="dense",
)
krylov_dae = SILVAImplicitDAEStep(
max_iter=5,
tol=1e-7,
linear_solver="gmres",
linear_max_iter=20,
linear_tol=1e-7,
)
y0 = torch.tensor([[1.0], [0.5]])
a0 = 0.25 * y0
dynamics = lambda y, a: -0.4 * y + a
constraint = lambda y, a: a - 0.25 * y
dense_step = dense_dae(y0, a0, 0.1, dynamics, constraint)
krylov_step = krylov_dae(y0, a0, 0.1, dynamics, constraint)
dae_error = max(
float((dense_step.differential - krylov_step.differential).abs().max()),
float((dense_step.algebraic - krylov_step.algebraic).abs().max()),
)
assert dae_error < 2e-5
print("dense/Newton-Krylov step error:", dae_error)
print("Krylov stage residual:", krylov_step.residual)
dense/Newton-Krylov step error: 0.0 Krylov stage residual: 1.0244548320770264e-08
10. Lazy Shards for a Periodic PDE¶
The generated field satisfies
$$ (-\Delta+m)u=f $$
on a periodic grid. The example writes aligned forcing and target tensors to independently loadable shards. The dataset keeps one shard cached per process, and the same loader configuration can select a distributed sampler for larger runs.
elliptic = make_periodic_elliptic_dataset(
samples=12,
height=8,
width=8,
modes=2,
seed=26,
)
shard_workspace = tempfile.TemporaryDirectory()
manifest = write_silva_tensor_shards(
{"x": elliptic.forcing, "y": elliptic.target},
Path(shard_workspace.name) / "periodic",
shard_size=5,
)
sharded_dataset = SILVAShardedTensorDataset(manifest)
runtime = runtime_for_tier(
"smoke",
per_device_batch_size=2,
gradient_accumulation_steps=2,
checkpoint_path=Path(shard_workspace.name) / "fno-checkpoint.pt",
)
train_loader = make_silva_dataloader(
sharded_dataset,
runtime.data_config(shuffle=False),
)
assert len(sharded_dataset) == 12
first_batch = next(iter(train_loader))
print("manifest:", manifest)
print("batch shapes:", {key: tuple(value.shape) for key, value in first_batch.items()})
manifest: /var/folders/6h/6crf402d0_sd_zgs_b139tb40000gp/T/tmpl83jlbha/periodic/shard-manifest.json
batch shapes: {'x': (2, 1, 8, 8), 'y': (2, 1, 8, 8)}
11. Train a Fourier Equilibrium and Resume¶
The Fourier family lifts the forcing, solves
$$ \begin{aligned} z^\star&=B_\theta(z^\star,P_\theta f),\\ \widehat u&=Q_\theta z^\star. \end{aligned} $$
and retains the source injection at every tied transition [43]. The first call trains one epoch and writes a complete checkpoint. The second call restores model, optimizer, history, scaler, and random-number states, then continues to epoch two.
import contextlib
import io
fno_model = build_scaled_silva(
"silva_fno_deq",
tier="smoke",
in_channels=1,
state_channels=2,
out_channels=1,
modes_height=2,
modes_width=2,
block_depth=1,
state_scale=0.05,
)
library_messages = io.StringIO()
with contextlib.redirect_stdout(library_messages), contextlib.redirect_stderr(
library_messages
):
first_run = fit_supervised(
fno_model,
train_loader,
config=runtime.train_config(
task="regression",
epochs=1,
optimizer="adam",
lr=2e-3,
),
)
resumed_run = fit_supervised(
fno_model,
train_loader,
config=runtime.train_config(
task="regression",
epochs=2,
optimizer="adam",
lr=2e-3,
),
)
assert runtime.checkpoint_path is not None and runtime.checkpoint_path.exists()
assert [row.epoch for row in resumed_run.history] == [1, 2]
print("loss history:", [row.train_loss for row in resumed_run.history])
loss history: [0.2529488479097684, 0.2388739436864853]
12. Diagnostics Are Separate Questions¶
A task metric does not prove that the implicit state converged, and a small fixed-point residual does not prove that the task was learned. For PDE work, report at least:
- task error in physical units or the official normalized metric;
- forward fixed-point residual and iteration count;
- backward linear residual when using implicit gradients;
- equation or conservation residual;
- resolution, precision, memory, and wall-clock protocol.
with torch.no_grad():
diagnostic = fno_model(elliptic.forcing[:2], return_result=True)
task_mse = torch.mean((diagnostic.output - elliptic.target[:2]).square())
print("task MSE:", float(task_mse))
print("fixed-point residual:", diagnostic.solver_result.residual)
print("solver iterations:", diagnostic.solver_result.iterations)
task MSE: 0.25076115131378174 fixed-point residual: 1.4526863196806516e-06 solver iterations: 7
losses = [row.train_loss for row in resumed_run.history]
figure, axes = plt.subplots(1, 3, figsize=(8.4, 2.5))
axes[0].imshow(elliptic.forcing[0, 0], cmap="viridis")
axes[0].set_title("forcing")
axes[1].imshow(elliptic.target[0, 0], cmap="viridis")
axes[1].set_title("exact field")
axes[2].plot(range(1, len(losses) + 1), losses, marker="o")
axes[2].set(xlabel="epoch", ylabel="MSE", yscale="log", xticks=[1, 2])
for axis in axes[:2]:
axis.axis("off")
figure.tight_layout()
plt.show()
13. Multi-Process Launch Contract¶
One process owns one accelerator. Initialize the process group,
create the ordinary SILVA module, call prepare_silva_model, and
build the loader from the same runtime. The family equation and
checkpoint format remain unchanged.
import os
import torch.distributed as dist
from silva_networks import prepare_silva_model, runtime_for_tier
dist.init_process_group(backend="nccl")
local_rank = int(os.environ["LOCAL_RANK"])
runtime = runtime_for_tier("full", device=f"cuda:{local_rank}")
model = prepare_silva_model(model, runtime, local_rank=local_rank)
loader = make_silva_dataloader(dataset, runtime.data_config())
Launch with torchrun --standalone --nproc-per-node=4 train.py.
The package advances the distributed sampler epoch and avoids
redundant gradient synchronization during accumulated microbatches.
14. One Point Can Contain Multiple Operators¶
A custom point may combine any shape-preserving modules:
$$ \begin{aligned} u_{\mathrm{local}}&=R_\theta(z)+C_\theta(z),\\ u_{\mathrm{global}}&=K_\theta(z)+A_\theta(z),\\ F_\theta(z;x)&=\tanh\!\big(P_xx+u_{\mathrm{local}}\\ &\qquad{}+u_{\mathrm{global}}\big). \end{aligned} $$
where $R$ may be residual, $C$ convolutional or U-Net-like, $K$
spectral or neural-operator based, and $A$ attention or global
context. SILVACortexLayer accepts the internal module graph;
SILVACortexNetwork links heterogeneous points with explicit
projections. The only non-negotiable rule is that the transition
returns the declared state shape.
15. Reproduce, Then Go Beyond¶
For a source-paper reproduction:
- select the canonical SILVA family and open its numbered source;
- preserve the official split, preprocessing, metric, and budget;
- match the source transition before adding SILVA branches;
- record solver and physical diagnostics independently;
- compare parameters, memory, iterations, and task quality;
- add new local/global/self branches one at a time;
- retain a configuration and checkpoint for every reported run.
The compact checks in this notebook establish implementation equivalence and execution wiring. They do not claim published benchmark scores. Full benchmark readiness still requires the cited data and protocol for the selected family.
summary = {
"canonical_families": len(guides),
"guide_errors": list(audit_silva_family_guides()),
"graph_equivalence_error": float(graph_error),
"attention_equivalence_error": max(attention_errors),
"measure_equivalence_error": float(measure_error),
"physics_derivative_error": float(derivative_error),
"dae_equivalence_error": dae_error,
"checkpoint_epochs": [row.epoch for row in resumed_run.history],
}
assert summary["guide_errors"] == []
summary
{'canonical_families': 64,
'guide_errors': [],
'graph_equivalence_error': 4.76837158203125e-07,
'attention_equivalence_error': 1.7881393432617188e-07,
'measure_equivalence_error': 2.2351741790771484e-08,
'physics_derivative_error': 1.341104507446289e-07,
'dae_equivalence_error': 0.0,
'checkpoint_epochs': [1, 2]}
shard_workspace.cleanup()
print("temporary notebook data removed")
temporary notebook data removed
From 26 Full Scale Silva 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 tensor solved to equilibrium |
| Condition | the observed input or source tensor |
| Repeated computation | the state-preserving transition evaluated by the root solver |
| Required invariants | shape, device, dtype, finiteness, and differentiability |
| Replaceable components | initializer, source encoder, transition, readout, and solver |
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.
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 fixed-point residual and task error against a deterministic target. 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 width, batch size, and data volume. 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.
notebook_reproduction_record = {
"notebook": '26_full_scale_silva.ipynb',
"state": 'the tensor solved to equilibrium',
"condition": 'the observed input or source tensor',
"transition": 'the state-preserving transition evaluated by the root solver',
"invariants": 'shape, device, dtype, finiteness, and differentiability',
"compact_metric": 'fixed-point residual and task error against a deterministic target',
"scale_axis": 'state width, batch size, and data volume',
}
assert all(notebook_reproduction_record.values())
notebook_reproduction_record
{'notebook': '26_full_scale_silva.ipynb',
'state': 'the tensor solved to equilibrium',
'condition': 'the observed input or source tensor',
'transition': 'the state-preserving transition evaluated by the root solver',
'invariants': 'shape, device, dtype, finiteness, and differentiability',
'compact_metric': 'fixed-point residual and task error against a deterministic target',
'scale_axis': 'state width, batch size, and data volume'}
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 tensor solved to equilibrium, and its repeated map, the state-preserving transition evaluated by the root solver. 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.
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
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 | fixed-point residual and task error against a deterministic target |
| 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 | shape, device, dtype, finiteness, and differentiability |
| Scale sweep | Change one of state width, batch size, and data volume 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.