Download notebook

Cortex Hierarchy with the SILVA Package¶

This notebook checks the cortex-style architecture path:

$$ x \xrightarrow{R_\phi} u_0 \xrightarrow{\alpha_1} z_1^\star \xrightarrow{\tanh} u_1 \xrightarrow{\alpha_2} z_2^\star \xrightarrow{R_\psi} \hat y. $$

Each cortex point solves

$$ z_\ell^\star = F_{\theta_\ell}(z_\ell^\star,u_{\ell-1}), $$

with damped solver steps

$$ z_{\ell,k+1} = (1-\alpha_\ell)z_{\ell,k} +\alpha_\ell F_{\theta_\ell}(z_{\ell,k},u_{\ell-1}). $$

The goal is to verify package behavior: flexible internal modules, different architectures at different points, different solvers, different alphas, gradients, and an image-cortex preset.

Numbered literature: [1], [4], [5], [15], [16], [17], [18], [19], [29]. Each number opens the complete citation and its primary external source.

In [1]:

from pathlib import Path
import importlib.util
import subprocess
import sys

IN_COLAB = "google.colab" in sys.modules
REPO_URL = "https://github.com/jseluis/silva-networks.git"

def find_local_silva_root():
    candidates = [
        Path.cwd(),
        Path("/content/silva-networks"),
        Path("/content/drive/MyDrive/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 IN_COLAB and 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()

from pathlib import Path import importlib.util import subprocess import sys IN_COLAB = "google.colab" in sys.modules REPO_URL = "https://github.com/jseluis/silva-networks.git" def find_local_silva_root(): candidates = [ Path.cwd(), Path("/content/silva-networks"), Path("/content/drive/MyDrive/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 IN_COLAB and 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()

In [2]:

import torch
import torch.nn.functional as F
import matplotlib.pyplot as plt

plt.rcParams.update({"figure.dpi": 300, "savefig.dpi": 300})

from silva_networks import (
    SILVACortexLayer,
    SILVACortexNetwork,
    SILVAImageCortexClassifier,
    SolverConfig,
    resolve_device,
    silva_equilibrium_model,
)

torch.manual_seed(11)
device = resolve_device("cuda" if torch.cuda.is_available() else "cpu")
device

import torch import torch.nn.functional as F import matplotlib.pyplot as plt plt.rcParams.update({"figure.dpi": 300, "savefig.dpi": 300}) from silva_networks import ( SILVACortexLayer, SILVACortexNetwork, SILVAImageCortexClassifier, SolverConfig, resolve_device, silva_equilibrium_model, ) torch.manual_seed(11) device = resolve_device("cuda" if torch.cuda.is_available() else "cpu") device

Out[2]:

device(type='cpu')

A Deep Internal Network Inside One Equilibrium Point¶

The internal transition network $B_\theta$ can contain many trainable layers:

$$ B_\theta(a) = B_{\theta,10}\circ\cdots\circ B_{\theta,1}(a). $$

The cortex point then adds this field to the stimulus and the optional interaction branches before the solver damping is applied.

In [3]:

def deep_state_network(dim, depth):
    modules = []
    for _ in range(depth):
        modules += [torch.nn.Linear(dim, dim), torch.nn.Tanh()]
    modules.append(torch.nn.Linear(dim, dim))
    return torch.nn.Sequential(*modules)


class StimulusGate(torch.nn.Module):
    def __init__(self, dim):
        super().__init__()
        self.gate = torch.nn.Linear(dim, dim)

    def forward(self, z, stimulus):
        return torch.sigmoid(self.gate(stimulus)) * z

def deep_state_network(dim, depth): modules = [] for _ in range(depth): modules += [torch.nn.Linear(dim, dim), torch.nn.Tanh()] modules.append(torch.nn.Linear(dim, dim)) return torch.nn.Sequential(*modules) class StimulusGate(torch.nn.Module): def __init__(self, dim): super().__init__() self.gate = torch.nn.Linear(dim, dim) def forward(self, z, stimulus): return torch.sigmoid(self.gate(stimulus)) * z

In [4]:

model = SILVACortexNetwork(
    [
        SILVACortexLayer(
            input_dim=5,
            state_dim=14,
            state_network=deep_state_network(14, depth=10),
            self_terms=torch.nn.Linear(14, 14, bias=False),
            interaction_terms=[StimulusGate(14)],
            config=SolverConfig(solver="picard", max_iter=6, alpha=0.5),
        ),
        SILVACortexLayer(
            input_encoder=torch.nn.Linear(14, 10),
            state_dim=10,
            state_network=torch.nn.Sequential(
                torch.nn.Linear(10, 20),
                torch.nn.GELU(),
                torch.nn.Linear(20, 10),
            ),
            config=SolverConfig(solver="anderson", max_iter=6, alpha=0.2, history=3),
            normalize=False,
        ),
    ],
    links="tanh",
    head=torch.nn.Linear(10, 2),
).to(device)

x = torch.randn(8, 5, device=device)
y = torch.tensor([0, 1, 0, 1, 0, 1, 0, 1], device=device)

result = model(x, return_results=True)
loss = torch.nn.functional.cross_entropy(result.output, y)
loss.backward()

print("logits shape:", tuple(result.output.shape))
print("state shapes:", [tuple(state.shape) for state in result.states])
print("solvers:", [r.solver for r in result.solver_results])
print("alphas:", [layer.config.alpha for layer in model.layers])
print("gradient reaches layer 1 encoder:", model.layers[0].input_encoder.weight.grad is not None)

model = SILVACortexNetwork( [ SILVACortexLayer( input_dim=5, state_dim=14, state_network=deep_state_network(14, depth=10), self_terms=torch.nn.Linear(14, 14, bias=False), interaction_terms=[StimulusGate(14)], config=SolverConfig(solver="picard", max_iter=6, alpha=0.5), ), SILVACortexLayer( input_encoder=torch.nn.Linear(14, 10), state_dim=10, state_network=torch.nn.Sequential( torch.nn.Linear(10, 20), torch.nn.GELU(), torch.nn.Linear(20, 10), ), config=SolverConfig(solver="anderson", max_iter=6, alpha=0.2, history=3), normalize=False, ), ], links="tanh", head=torch.nn.Linear(10, 2), ).to(device) x = torch.randn(8, 5, device=device) y = torch.tensor([0, 1, 0, 1, 0, 1, 0, 1], device=device) result = model(x, return_results=True) loss = torch.nn.functional.cross_entropy(result.output, y) loss.backward() print("logits shape:", tuple(result.output.shape)) print("state shapes:", [tuple(state.shape) for state in result.states]) print("solvers:", [r.solver for r in result.solver_results]) print("alphas:", [layer.config.alpha for layer in model.layers]) print("gradient reaches layer 1 encoder:", model.layers[0].input_encoder.weight.grad is not None)

logits shape: (8, 2)
state shapes: [(8, 14), (8, 10)]
solvers: ['picard', 'anderson']
alphas: [0.5, 0.2]
gradient reaches layer 1 encoder: True

Residual Curves¶

The two equilibrium points solve different equations. Their residual curves are not expected to be identical because their transition maps, solvers, and alphas are different.

In [5]:

plt.figure(figsize=(6, 3))
for i, solver_result in enumerate(result.solver_results, start=1):
    plt.plot(solver_result.residuals, marker="o", label=f"point {i}: {solver_result.solver}")
plt.yscale("log")
plt.xlabel("solver step")
plt.ylabel("residual")
plt.title("Cortex hierarchy residuals")
plt.legend()
plt.tight_layout()

plt.figure(figsize=(6, 3)) for i, solver_result in enumerate(result.solver_results, start=1): plt.plot(solver_result.residuals, marker="o", label=f"point {i}: {solver_result.solver}") plt.yscale("log") plt.xlabel("solver step") plt.ylabel("residual") plt.title("Cortex hierarchy residuals") plt.legend() plt.tight_layout()

Spatial Architecture Inside a SILVA Point¶

The equilibrium state may be an image tensor rather than a feature vector. In this example, the first point has state shape (batch, 4, 8, 8) and evaluates a residual convolutional block plus a U-Net-shaped transition during every solver iteration. The U-Net may downsample internally, but it restores the equilibrium-state shape before returning:

$$ F_{\theta_1}:\mathbb R^{4\times8\times8} \rightarrow\mathbb R^{4\times8\times8}. $$

The solved spatial state is flattened and passed to a second SILVA point with a different vector architecture and solver.

In [6]:

class SILVAResidualConvTransition(torch.nn.Module):
    def __init__(self, channels):
        super().__init__()
        self.conv1 = torch.nn.Conv2d(channels, channels, 3, padding=1)
        self.conv2 = torch.nn.Conv2d(channels, channels, 3, padding=1)
        self.norm1 = torch.nn.GroupNorm(1, channels)
        self.norm2 = torch.nn.GroupNorm(1, channels)

    def forward(self, z):
        update = F.gelu(self.norm1(self.conv1(z)))
        return z + 0.25 * self.norm2(self.conv2(update))


class SILVATinyUNetTransition(torch.nn.Module):
    def __init__(self, channels):
        super().__init__()
        expanded = 2 * channels
        self.encoder = SILVAResidualConvTransition(channels)
        self.down = torch.nn.Conv2d(channels, expanded, 3, stride=2, padding=1)
        self.bottleneck = SILVAResidualConvTransition(expanded)
        self.up = torch.nn.ConvTranspose2d(expanded, channels, 2, stride=2)
        self.decoder = torch.nn.Conv2d(2 * channels, channels, 3, padding=1)

    def forward(self, z):
        skip = self.encoder(z)
        low = self.bottleneck(F.gelu(self.down(skip)))
        up = self.up(low)
        if up.shape[-2:] != skip.shape[-2:]:
            up = F.interpolate(up, size=skip.shape[-2:], mode="bilinear", align_corners=False)
        return 0.25 * torch.tanh(self.decoder(torch.cat([skip, up], dim=1)))


class SILVASpatialToVectorLink(torch.nn.Module):
    def forward(self, z):
        return z.flatten(start_dim=1)


def make_tiny_pattern_dataset(samples=24, image_size=8):
    generator = torch.Generator().manual_seed(57)
    images = 0.04 * torch.randn(samples, 1, image_size, image_size, generator=generator)
    labels = torch.arange(samples) % 2
    center = image_size // 2
    for index, label in enumerate(labels):
        if int(label) == 0:
            images[index, 0, :, center - 1:center + 1] += 1.0
        else:
            images[index, 0, center - 1:center + 1, :] += 1.0
    return images, labels

class SILVAResidualConvTransition(torch.nn.Module): def __init__(self, channels): super().__init__() self.conv1 = torch.nn.Conv2d(channels, channels, 3, padding=1) self.conv2 = torch.nn.Conv2d(channels, channels, 3, padding=1) self.norm1 = torch.nn.GroupNorm(1, channels) self.norm2 = torch.nn.GroupNorm(1, channels) def forward(self, z): update = F.gelu(self.norm1(self.conv1(z))) return z + 0.25 * self.norm2(self.conv2(update)) class SILVATinyUNetTransition(torch.nn.Module): def __init__(self, channels): super().__init__() expanded = 2 * channels self.encoder = SILVAResidualConvTransition(channels) self.down = torch.nn.Conv2d(channels, expanded, 3, stride=2, padding=1) self.bottleneck = SILVAResidualConvTransition(expanded) self.up = torch.nn.ConvTranspose2d(expanded, channels, 2, stride=2) self.decoder = torch.nn.Conv2d(2 * channels, channels, 3, padding=1) def forward(self, z): skip = self.encoder(z) low = self.bottleneck(F.gelu(self.down(skip))) up = self.up(low) if up.shape[-2:] != skip.shape[-2:]: up = F.interpolate(up, size=skip.shape[-2:], mode="bilinear", align_corners=False) return 0.25 * torch.tanh(self.decoder(torch.cat([skip, up], dim=1))) class SILVASpatialToVectorLink(torch.nn.Module): def forward(self, z): return z.flatten(start_dim=1) def make_tiny_pattern_dataset(samples=24, image_size=8): generator = torch.Generator().manual_seed(57) images = 0.04 * torch.randn(samples, 1, image_size, image_size, generator=generator) labels = torch.arange(samples) % 2 center = image_size // 2 for index, label in enumerate(labels): if int(label) == 0: images[index, 0, :, center - 1:center + 1] += 1.0 else: images[index, 0, center - 1:center + 1, :] += 1.0 return images, labels

In [7]:

channels = 4
spatial_point = SILVACortexLayer(
    input_encoder=torch.nn.Conv2d(1, channels, 3, padding=1),
    state_network=torch.nn.Sequential(
        SILVAResidualConvTransition(channels),
        SILVATinyUNetTransition(channels),
    ),
    normalizer=torch.nn.GroupNorm(1, channels),
    config=SolverConfig(solver="picard", max_iter=3, alpha=0.35),
)
vector_point = SILVACortexLayer(
    input_dim=channels * 8 * 8,
    state_dim=12,
    state_network=torch.nn.Sequential(
        torch.nn.Linear(12, 24),
        torch.nn.GELU(),
        torch.nn.Linear(24, 12),
    ),
    config=SolverConfig(solver="anderson", max_iter=3, alpha=0.2, history=2),
)
spatial_model = silva_equilibrium_model(
    "silva_cortex_network",
    layers=[spatial_point, vector_point],
    links=[SILVASpatialToVectorLink()],
    head=torch.nn.Linear(12, 2),
).to(device)

pattern_images, pattern_labels = make_tiny_pattern_dataset()
pattern_images = pattern_images.to(device)
pattern_labels = pattern_labels.to(device)
optimizer = torch.optim.Adam(spatial_model.parameters(), lr=2e-2)

for _ in range(4):
    spatial_result = spatial_model(pattern_images, return_results=True)
    spatial_loss = F.cross_entropy(spatial_result.output, pattern_labels)
    optimizer.zero_grad()
    spatial_loss.backward()
    optimizer.step()

spatial_accuracy = (spatial_result.output.argmax(dim=1) == pattern_labels).float().mean()
spatial_gradients = [
    spatial_model.layers[0].input_encoder.weight.grad is not None,
    spatial_model.layers[1].input_encoder.weight.grad is not None,
]
if not all(spatial_gradients):
    raise RuntimeError("gradients did not reach both SILVA equilibrium points")
print("state shapes:", [tuple(state.shape) for state in spatial_result.states])
print("solvers:", [item.solver for item in spatial_result.solver_results])
print("loss:", float(spatial_loss.detach().cpu()))
print("accuracy:", float(spatial_accuracy.detach().cpu()))
print("point gradients:", spatial_gradients)

channels = 4 spatial_point = SILVACortexLayer( input_encoder=torch.nn.Conv2d(1, channels, 3, padding=1), state_network=torch.nn.Sequential( SILVAResidualConvTransition(channels), SILVATinyUNetTransition(channels), ), normalizer=torch.nn.GroupNorm(1, channels), config=SolverConfig(solver="picard", max_iter=3, alpha=0.35), ) vector_point = SILVACortexLayer( input_dim=channels * 8 * 8, state_dim=12, state_network=torch.nn.Sequential( torch.nn.Linear(12, 24), torch.nn.GELU(), torch.nn.Linear(24, 12), ), config=SolverConfig(solver="anderson", max_iter=3, alpha=0.2, history=2), ) spatial_model = silva_equilibrium_model( "silva_cortex_network", layers=[spatial_point, vector_point], links=[SILVASpatialToVectorLink()], head=torch.nn.Linear(12, 2), ).to(device) pattern_images, pattern_labels = make_tiny_pattern_dataset() pattern_images = pattern_images.to(device) pattern_labels = pattern_labels.to(device) optimizer = torch.optim.Adam(spatial_model.parameters(), lr=2e-2) for _ in range(4): spatial_result = spatial_model(pattern_images, return_results=True) spatial_loss = F.cross_entropy(spatial_result.output, pattern_labels) optimizer.zero_grad() spatial_loss.backward() optimizer.step() spatial_accuracy = (spatial_result.output.argmax(dim=1) == pattern_labels).float().mean() spatial_gradients = [ spatial_model.layers[0].input_encoder.weight.grad is not None, spatial_model.layers[1].input_encoder.weight.grad is not None, ] if not all(spatial_gradients): raise RuntimeError("gradients did not reach both SILVA equilibrium points") print("state shapes:", [tuple(state.shape) for state in spatial_result.states]) print("solvers:", [item.solver for item in spatial_result.solver_results]) print("loss:", float(spatial_loss.detach().cpu())) print("accuracy:", float(spatial_accuracy.detach().cpu())) print("point gradients:", spatial_gradients)

state shapes: [(24, 4, 8, 8), (24, 12)]
solvers: ['picard', 'anderson']
loss: 0.127569779753685
accuracy: 1.0
point gradients: [True, True]

The module contract for one SILVA point is:

1. the completed transition returns the equilibrium-state shape;
2. the transition is deterministic during one solve;
3. tensors remain on the same device and dtype;
4. every operation supports the selected backward mode.

For spatial states, GroupNorm avoids mutable batch statistics. Random masks must remain consistent during repeated transition evaluations. Larger internal architectures may require stronger damping, residual scaling, or spectral normalization; residual curves remain the direct convergence check.

Alpha Sweep¶

For a fixed transition, smaller $\alpha$ mixes less of the new field value at each iteration:

$$ z_{k+1}-z_k = \alpha(F_\theta(z_k,x)-z_k). $$

The step size is proportional to $\alpha$. This is why a fast first point and a slower second point can be encoded directly by separate SolverConfig objects.

In [8]:

base_x = torch.randn(4, 3, device=device)
alpha_summaries = []
for alpha in [0.2, 0.5, 0.8]:
    layer = SILVACortexLayer(
        input_dim=3,
        state_dim=6,
        state_network=torch.nn.Linear(6, 6),
        config=SolverConfig(solver="picard", max_iter=5, alpha=alpha),
        normalize=False,
    ).to(device)
    out = layer(base_x, return_result=True)
    alpha_summaries.append((alpha, out.residuals))

for alpha, residuals in alpha_summaries:
    print(alpha, [round(v, 6) for v in residuals])

base_x = torch.randn(4, 3, device=device) alpha_summaries = [] for alpha in [0.2, 0.5, 0.8]: layer = SILVACortexLayer( input_dim=3, state_dim=6, state_network=torch.nn.Linear(6, 6), config=SolverConfig(solver="picard", max_iter=5, alpha=alpha), normalize=False, ).to(device) out = layer(base_x, return_result=True) alpha_summaries.append((alpha, out.residuals)) for alpha, residuals in alpha_summaries: print(alpha, [round(v, 6) for v in residuals])

0.2 [2.634781, 2.085235, 1.663175, 1.335796, 1.078927]
0.5 [2.859459, 1.547449, 0.841202, 0.46107, 0.256476]
0.8 [2.555252, 1.062598, 0.390488, 0.142308, 0.054648]

In [9]:

plt.figure(figsize=(6, 3))
for alpha, residuals in alpha_summaries:
    plt.plot(residuals, marker="o", label=f"alpha={alpha}")
plt.yscale("log")
plt.xlabel("solver step")
plt.ylabel("residual")
plt.title("Effect of solver damping")
plt.legend()
plt.tight_layout()

plt.figure(figsize=(6, 3)) for alpha, residuals in alpha_summaries: plt.plot(residuals, marker="o", label=f"alpha={alpha}") plt.yscale("log") plt.xlabel("solver step") plt.ylabel("residual") plt.title("Effect of solver damping") plt.legend() plt.tight_layout()

Image Cortex Preset¶

SILVAImageCortexClassifier packages the convolutional-retina path:

$$ u_0=C_\psi(x), \qquad z_1^\star \xrightarrow{\tanh} z_2^\star. $$

The local branch can be the dynamic hidden-channel kNN term, and the global branch can be per-sample channel attention. Both remain configurable.

In [10]:

image_model = SILVAImageCortexClassifier(
    in_channels=3,
    hidden_dim=[8, 6],
    num_classes=2,
    image_size=8,
    attention_mode="simple",
    graph_mode="GAT",
    k_neighbors=2,
    alphas=(0.5, 0.2),
    max_iter=3,
    internal_depth=2,
    self_interaction=True,
    dropout=0.0,
).to(device)

images = torch.randn(4, 3, 8, 8, device=device)
labels = torch.tensor([0, 1, 0, 1], device=device)
image_result = image_model(images, return_results=True)
image_loss = torch.nn.functional.cross_entropy(image_result.output, labels)
image_loss.backward()

print("image logits:", tuple(image_result.output.shape))
print("image states:", [tuple(state.shape) for state in image_result.states])
print("image solvers:", [r.solver for r in image_result.solver_results])
print("retina gradient:", image_model.retina.conv1.weight.grad is not None)

image_model = SILVAImageCortexClassifier( in_channels=3, hidden_dim=[8, 6], num_classes=2, image_size=8, attention_mode="simple", graph_mode="GAT", k_neighbors=2, alphas=(0.5, 0.2), max_iter=3, internal_depth=2, self_interaction=True, dropout=0.0, ).to(device) images = torch.randn(4, 3, 8, 8, device=device) labels = torch.tensor([0, 1, 0, 1], device=device) image_result = image_model(images, return_results=True) image_loss = torch.nn.functional.cross_entropy(image_result.output, labels) image_loss.backward() print("image logits:", tuple(image_result.output.shape)) print("image states:", [tuple(state.shape) for state in image_result.states]) print("image solvers:", [r.solver for r in image_result.solver_results]) print("retina gradient:", image_model.retina.conv1.weight.grad is not None)

image logits: (4, 2)
image states: [(4, 8), (4, 6)]
image solvers: ['picard', 'picard']
retina gradient: True

What This Enables¶

The same grammar covers the article configurations and extensions:

Need	Package control
convolutional front end	SILVAImageCortexClassifier.retina or custom input_encoder
many layers inside one point	state_network=nn.Sequential(...) or a module list
different point architectures	one SILVACortexLayer per point
different alphas	one SolverConfig(alpha=...) per point
different solvers	Picard, Anderson, or Broyden per point
local/global/self ablations	pass modules, omit modules, or use zero modules
user datasets	adapt to tensors, then call the same PyTorch modules

Citation¶

Dr. Jose Luis Silva. SILVA Networks. Version 1.2.2. MIT License. https://github.com/jseluis/silva-networks https://doi.org/10.5281/zenodo.21770098

From 11 Cortex Hierarchy 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	one image tensor per resolution or linked SILVA point
Condition	image features and per-scale source injections
Repeated computation	shape-preserving convolutional, U-Net, attention, or multiscale fusion blocks
Required invariants	channel/spatial shape at every scale and deterministic fusion
Replaceable components	stem, per-scale injections, transition blocks, links, task head, and solvers

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 [11]:

silva-extension-curriculum

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 task error, per-scale residuals, and gradient agreement. 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 image resolution, channels, scales, internal depth, and batch size. 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¶

1. Replace one component from this notebook while preserving its state and domain invariants.
2. Write the replacement first as an independent reference function, then as a module, and assert one-step equivalence.
3. Compare two solver configurations on the identical trained transition.
4. Add a compact baseline and a predeclared metric threshold.
5. Create a full-scale configuration without weakening the compact tests.

The complete authoring protocol is documented in Extending SILVA.

In [12]:

silva-extension-curriculum

notebook_reproduction_record = {
    "notebook": '11_cortex_hierarchy.ipynb',
    "state": 'one image tensor per resolution or linked SILVA point',
    "condition": 'image features and per-scale source injections',
    "transition": 'shape-preserving convolutional, U-Net, attention, or multiscale fusion blocks',
    "invariants": 'channel/spatial shape at every scale and deterministic fusion',
    "compact_metric": 'task error, per-scale residuals, and gradient agreement',
    "scale_axis": 'image resolution, channels, scales, internal depth, and batch size',
}
assert all(notebook_reproduction_record.values())
notebook_reproduction_record

notebook_reproduction_record = { "notebook": '11_cortex_hierarchy.ipynb', "state": 'one image tensor per resolution or linked SILVA point', "condition": 'image features and per-scale source injections', "transition": 'shape-preserving convolutional, U-Net, attention, or multiscale fusion blocks', "invariants": 'channel/spatial shape at every scale and deterministic fusion', "compact_metric": 'task error, per-scale residuals, and gradient agreement', "scale_axis": 'image resolution, channels, scales, internal depth, and batch size', } assert all(notebook_reproduction_record.values()) notebook_reproduction_record

Out[12]:

{'notebook': '11_cortex_hierarchy.ipynb',
 'state': 'one image tensor per resolution or linked SILVA point',
 'condition': 'image features and per-scale source injections',
 'transition': 'shape-preserving convolutional, U-Net, attention, or multiscale fusion blocks',
 'invariants': 'channel/spatial shape at every scale and deterministic fusion',
 'compact_metric': 'task error, per-scale residuals, and gradient agreement',
 'scale_axis': 'image resolution, channels, scales, internal depth, and batch size'}

Worked Convergence and Sensitivity Study¶

The preceding example demonstrates one configured solve. This additional study changes the spatial 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, one image tensor per resolution or linked SILVA point, and its repeated map, shape-preserving convolutional, U-Net, attention, or multiscale fusion blocks. 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 [13]:

silva-extension-curriculum

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('spatial 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('spatial 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)

spatial 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 [14]:

silva-extension-curriculum

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('spatial 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('spatial 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 spatial 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	task error, per-scale residuals, and gradient agreement
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	channel/spatial shape at every scale and deterministic fusion
Scale sweep	Change one of image resolution, channels, scales, internal depth, and batch size 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.

Where to Go Next¶

Question	Page
How are linked equilibrium points derived?	Cortex Hierarchies
Where is a smaller stacked architecture executed?	Stacked Architecture Example
Which hierarchy containers are public?	Architectures API