Skip to content

Diagnostics

SILVA layers are useful when the fixed-point dynamics can be inspected. The package exposes residual curves, Jacobian products, spectral-radius estimates, and Lyapunov-style energy traces as ordinary Python functions.

The stochastic Jacobian estimate follows Hutchinson [14], and the equilibrium regularization interpretation follows Jacobian-regularized DEQ [6].

Residual

For a transition \(f\), the fixed-point residual is

\[ r(z)=f(z)-z, \qquad \|r(z)\|_2=\|f(z)-z\|_2. \]

residual_curve records this norm during a solve:

from silva_networks import SolverConfig, residual_curve

curve = residual_curve(f, z0, SolverConfig(max_iter=20, alpha=0.5))

Damped Spectral Radius

The executed damped update is

\[ T_\alpha(z)=(1-\alpha)z+\alpha f(z). \]

Its local Jacobian is

\[ J_{T_\alpha}(z^\star) = (1-\alpha)I+\alpha J_f(z^\star). \]

damped_spectral_radius estimates

\[ \rho(J_{T_\alpha}(z^\star)) \]

with VJP-based power iteration. This is the quantity used for local fixed-point stability diagnostics.

Full Diagnostic Record

For a layer transition \(f\), a reproducible experiment record should contain:

\[ \left( z_K,\, \{\|f(z_k)-z_k\|_2\}_{k=0}^{K},\, \rho(J_{T_\alpha}(z_K)),\, \widehat{\|J_f(z_K)\|_F},\, \{E(z_k)\}_{k=0}^{K} \right). \]

In package terms:

from silva_networks import (
    SolverConfig,
    damped_spectral_radius,
    hutchinson_jacobian_norm,
    lyapunov_quadratic_energy,
    solve_with_energy,
)

config = SolverConfig(alpha=0.5, max_iter=20, tol=1e-6)
report = solve_with_energy(
    f,
    z0,
    energy_fn=lambda z: lyapunov_quadratic_energy(z, interaction(z)),
    config=config,
    include_stability=True,
)
rho = damped_spectral_radius(f, report.result.z, alpha=config.alpha)
fro = hutchinson_jacobian_norm(f, report.result.z, samples=8, squared=False)

Lyapunov-Style Energy

The package provides the quadratic alignment diagnostic

\[ E_i(z)=\|z_i\|_2^2-z_i^\top h_i, \]

where \(h_i\) is a local, global, or combined interaction evaluated at the same state. In code:

from silva_networks import lyapunov_quadratic_energy

energy = lyapunov_quadratic_energy(z, local_update + global_update)

This is a diagnostic proxy. A rigorous Lyapunov certificate requires the assumptions of the specific dynamical system being studied.

Solve With Energy

solve_with_energy runs a fixed-point solve while evaluating an energy function on each iterate:

from silva_networks import solve_with_energy

report = solve_with_energy(
    f,
    z0,
    energy_fn=lambda z: lyapunov_quadratic_energy(z, interaction(z)),
    include_stability=True,
)

report.result.z
report.energies
report.energy_deltas
report.stability.spectral_radius

descent_fraction(report.energies) returns the fraction of consecutive energy steps that did not increase.

Failure Modes

Symptom Mathematical sign First response
residual grows \(\|r_{k+1}\|>\|r_k\|\) repeatedly lower alpha, inspect operator scale
residual alternates linearized modes near unit radius increase damping, try Anderson with ridge
adjoint solve is ill-conditioned \(I-J_f^\top\) nearly singular inspect \(\rho(J_f)\), regularize or simplify
energy rises \(E_{k+1}-E_k>0\) often treat energy as warning, not proof
task metric improves but residual is high finite-depth behavior dominates report finite-solve budget honestly
graph batch leaks context global term ignores batch validate batch and use package global operators

Claim Scale

Evidence Claim scale
one residual curve this input solved under this budget
residuals over a validation split this configuration is numerically reliable on sampled data
residuals plus \(\rho(J_{T_\alpha})<1\) local stability evidence near computed states
full ablation plus diagnostics model mechanism claim
theorem-level assumptions plus diagnostics certificate or guarantee

ResidualEnergyReport dataclass

Diagnostics collected during a fixed-point solve.

Source code in src/silva_networks/diagnostics.py 
14
15
16
17
18
19
20
21
22
23
24
25
@dataclass
class ResidualEnergyReport:
    """Diagnostics collected during a fixed-point solve."""

    result: SolverResult
    energies: list[float]
    energy_deltas: list[float]
    stability: StabilityReport | None = None

    @property
    def final_energy(self) -> float:
        return self.energies[-1] if self.energies else float("nan")

damped_spectral_radius 

damped_spectral_radius(f, z_star, alpha, *, iters=20)

Estimate the spectral radius of the executed damped update.

This estimates rho((1-alpha)I + alpha J_f(z_star)) by calling the package's VJP-based spectral-radius estimator on T_alpha.

Source code in src/silva_networks/diagnostics.py 
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
def damped_spectral_radius(
    f: Callable[[Tensor], Tensor],
    z_star: Tensor,
    alpha: float,
    *,
    iters: int = 20,
) -> float:
    """Estimate the spectral radius of the executed damped update.

    This estimates ``rho((1-alpha)I + alpha J_f(z_star))`` by calling the
    package's VJP-based spectral-radius estimator on ``T_alpha``.
    """

    report = stability_report(damped_update(f, alpha), z_star, iters=iters)
    return report.spectral_radius

damped_update 

damped_update(f, alpha)

Return the executed Picard-style update T_alpha(z).

Source code in src/silva_networks/diagnostics.py 
28
29
30
31
32
33
34
35
36
37
def damped_update(
    f: Callable[[Tensor], Tensor],
    alpha: float,
) -> Callable[[Tensor], Tensor]:
    """Return the executed Picard-style update ``T_alpha(z)``."""

    def transition(z: Tensor) -> Tensor:
        return (1.0 - alpha) * z + alpha * f(z)

    return transition

descent_fraction 

descent_fraction(energies, tolerance=0.0)

Fraction of consecutive diagnostic energy changes that are nonincreasing.

Source code in src/silva_networks/diagnostics.py 
86
87
88
89
90
91
92
def descent_fraction(energies: list[float], tolerance: float = 0.0) -> float:
    """Fraction of consecutive diagnostic energy changes that are nonincreasing."""

    deltas = energy_deltas(energies)
    if not deltas:
        return float("nan")
    return sum(delta <= tolerance for delta in deltas) / len(deltas)

energy_deltas 

energy_deltas(energies)

Return consecutive energy changes E_{k+1}-E_k.

Source code in src/silva_networks/diagnostics.py 
80
81
82
83
def energy_deltas(energies: list[float]) -> list[float]:
    """Return consecutive energy changes ``E_{k+1}-E_k``."""

    return [energies[i + 1] - energies[i] for i in range(len(energies) - 1)]

lyapunov_quadratic_energy 

lyapunov_quadratic_energy(z, interaction, *, reduction='mean')

Quadratic alignment energy used in SILVA diagnostics.

Per row,

\[ E_i = \|z_i\|_2^2 - z_i^\top h_i, \]

where h is a local, global, or local+global interaction evaluated at the same state. This is a monitoring quantity; a rigorous Lyapunov certificate requires the assumptions of the specific dynamical system.

Source code in src/silva_networks/diagnostics.py 
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
def lyapunov_quadratic_energy(
    z: Tensor,
    interaction: Tensor,
    *,
    reduction: str = "mean",
) -> Tensor:
    r"""Quadratic alignment energy used in SILVA diagnostics.

    Per row,

    $$
    E_i = \|z_i\|_2^2 - z_i^\top h_i,
    $$

    where ``h`` is a local, global, or local+global interaction evaluated at
    the same state. This is a monitoring quantity; a rigorous Lyapunov
    certificate requires the assumptions of the specific dynamical system.
    """

    values = (z * z).sum(dim=-1) - (z * interaction).sum(dim=-1)
    if reduction == "none":
        return values
    if reduction == "sum":
        return values.sum()
    if reduction == "mean":
        return values.mean()
    raise ValueError("reduction must be 'none', 'sum', or 'mean'")

residual_curve 

residual_curve(f, z0, config=None)

Solve once and return ||f(z_k)-z_k||_2 values.

Source code in src/silva_networks/diagnostics.py 
40
41
42
43
44
45
46
47
48
def residual_curve(
    f: Callable[[Tensor], Tensor],
    z0: Tensor,
    config: SolverConfig | None = None,
) -> list[float]:
    """Solve once and return ``||f(z_k)-z_k||_2`` values."""

    result = fixed_point(f, z0, config)
    return result.residuals

solve_with_energy 

solve_with_energy(f, z0, energy_fn, config=None, *, include_stability=False, stability_samples=8, stability_iters=20)

Solve z=f(z) while collecting an energy trace.

energy_fn is evaluated on each iterate before applying the next solver step. This makes the function useful for custom Lyapunov-style diagnostics on a user's own SILVA module.

Source code in src/silva_networks/diagnostics.py 
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
def solve_with_energy(
    f: Callable[[Tensor], Tensor],
    z0: Tensor,
    energy_fn: Callable[[Tensor], Tensor],
    config: SolverConfig | None = None,
    *,
    include_stability: bool = False,
    stability_samples: int = 8,
    stability_iters: int = 20,
) -> ResidualEnergyReport:
    """Solve ``z=f(z)`` while collecting an energy trace.

    ``energy_fn`` is evaluated on each iterate before applying the next solver
    step. This makes the function useful for custom Lyapunov-style diagnostics
    on a user's own SILVA module.
    """

    energies: list[float] = []

    def tracked_f(z: Tensor) -> Tensor:
        energies.append(float(energy_fn(z).detach().cpu()))
        return f(z)

    result = fixed_point(tracked_f, z0, config)
    stability = (
        stability_report(f, result.z, samples=stability_samples, iters=stability_iters)
        if include_stability
        else None
    )
    return ResidualEnergyReport(
        result=result,
        energies=energies,
        energy_deltas=energy_deltas(energies),
        stability=stability,
    )

Where to Go Next

Question Page
How should residual and stability traces be interpreted? Interactive Diagnostics Lab
Which Jacobian estimates support the diagnostics? Jacobians API
Which solver result fields supply the traces? Solvers API