Educational NumPy API

The educational module mirrors the PyTorch API with small NumPy functions. Use it when you want to see the algebra without autograd, modules, batching, or GPU concerns.

Why This Module Exists

The PyTorch package is the production path. The NumPy helpers are for hand-sized derivations:

Helper	Mathematical object
np_picard	Damped fixed-point iteration
np_finite_difference_jacobian	Central-difference Jacobian
np_exact_tanh_affine_jacobian	Closed-form Jacobian for tanh(Wz + s)
np_power_iteration	Dominant mode estimate for a materialized matrix
np_implicit_gradient	Explicit adjoint solve for a small DEQ

Minimal Fixed Point

For

\[ f(z)=\tanh(Wz+s), \]

Picard iteration computes

\[ z_{k+1}=(1-\alpha)z_k+\alpha f(z_k). \]

import numpy as np
from silva_networks import np_picard

W = np.array([[0.2, 0.1], [-0.1, 0.25]])
s = np.array([0.5, -0.2])
trace = np_picard(lambda z: np.tanh(W @ z + s), np.zeros(2), alpha=0.8)

The result stores the final z and the residual curve.

Exact Jacobian for tanh(Wz + s)

Let

\[ u=Wz+s, \qquad f(z)=\tanh(u). \]

Because

\[ \frac{d}{du}\tanh(u)=1-\tanh^2(u), \]

the state Jacobian is

\[ J_f(z) = \operatorname{diag}\left(1-\tanh^2(Wz+s)\right)W. \]

from silva_networks import np_exact_tanh_affine_jacobian

J = np_exact_tanh_affine_jacobian(W, trace.z, s)

Small Implicit Gradient

At a solved equilibrium, the total derivative obeys

\[ (I-J_f)\frac{dz^\star}{d\theta} = \frac{\partial f}{\partial \theta}. \]

For a loss gradient \(g=\partial \mathcal L/\partial z^\star\), solve

\[ (I-J_f^\top)\lambda=g, \]

then compute

\[ \frac{\partial \mathcal L}{\partial \theta} = \lambda^\top \frac{\partial f}{\partial \theta}. \]

from silva_networks import np_implicit_gradient

grad_theta = np_implicit_gradient(J, grad_z, df_dtheta)

This is the explicit small-matrix version of the adjoint system used by the PyTorch diagnostics.

The corresponding equilibrium, implicit-function, and numerical-method sources are listed in Paper and References. Use Mathematical Foundations for the full derivations and tensor notation.

NumpySolverTrace dataclass

Transparent NumPy trace for hand-sized fixed-point examples.

Source code in src/silva_networks/educational.py

@dataclass(frozen=True)
class NumpySolverTrace:
    """Transparent NumPy trace for hand-sized fixed-point examples."""

    z: np.ndarray
    residuals: list[float]

    @property
    def converged(self) -> bool:
        return bool(self.residuals and self.residuals[-1] < 1e-8)

np_exact_tanh_affine_jacobian

np_exact_tanh_affine_jacobian(W, z, s)

Jacobian of f(z) = tanh(W z + s).

Source code in src/silva_networks/educational.py

def np_exact_tanh_affine_jacobian(W: np.ndarray, z: np.ndarray, s: np.ndarray) -> np.ndarray:
    """Jacobian of ``f(z) = tanh(W z + s)``."""

    pre = W @ z + s
    D = np.diag(1.0 - np.tanh(pre) ** 2)
    return D @ W

np_finite_difference_jacobian

np_finite_difference_jacobian(f, z, eps=1e-05)

Central-difference Jacobian with columns J[:, i] = d f / d z_i.

Source code in src/silva_networks/educational.py

def np_finite_difference_jacobian(
    f: Callable[[np.ndarray], np.ndarray],
    z: np.ndarray,
    eps: float = 1e-5,
) -> np.ndarray:
    """Central-difference Jacobian with columns ``J[:, i] = d f / d z_i``."""

    z = np.array(z, dtype=float)
    y0 = np.asarray(f(z), dtype=float).reshape(-1)
    J = np.zeros((y0.size, z.size), dtype=float)
    for i in range(z.size):
        step = np.zeros_like(z).reshape(-1)
        step[i] = eps
        step = step.reshape(z.shape)
        yp = np.asarray(f(z + step), dtype=float).reshape(-1)
        ym = np.asarray(f(z - step), dtype=float).reshape(-1)
        J[:, i] = (yp - ym) / (2.0 * eps)
    return J

np_implicit_gradient

np_implicit_gradient(J, grad_z, df_dtheta)

Compute grad_theta L = lambda^T df/dtheta from the DEQ adjoint solve.

Source code in src/silva_networks/educational.py

def np_implicit_gradient(J: np.ndarray, grad_z: np.ndarray, df_dtheta: np.ndarray) -> np.ndarray:
    """Compute ``grad_theta L = lambda^T df/dtheta`` from the DEQ adjoint solve."""

    lam = np.linalg.solve(np.eye(J.shape[0]) - J.T, grad_z.reshape(-1))
    return lam @ df_dtheta.reshape(J.shape[0], -1)

np_picard

np_picard(f, z0, max_iter=50, tol=1e-08, alpha=1.0)

Damped Picard iteration written as visible NumPy linear algebra.

Source code in src/silva_networks/educational.py

def np_picard(
    f: Callable[[np.ndarray], np.ndarray],
    z0: np.ndarray,
    max_iter: int = 50,
    tol: float = 1e-8,
    alpha: float = 1.0,
) -> NumpySolverTrace:
    """Damped Picard iteration written as visible NumPy linear algebra."""

    z = np.array(z0, dtype=float)
    residuals: list[float] = []
    for _ in range(max_iter):
        fz = np.asarray(f(z), dtype=float)
        residuals.append(float(np.linalg.norm(fz - z)))
        z = (1.0 - alpha) * z + alpha * fz
        if residuals[-1] < tol:
            break
    return NumpySolverTrace(z=z, residuals=residuals)

np_power_iteration

np_power_iteration(A, iters=50)

Dominant singular/eigenmode magnitude estimate for a materialized matrix.

Source code in src/silva_networks/educational.py

def np_power_iteration(A: np.ndarray, iters: int = 50) -> tuple[float, np.ndarray]:
    """Dominant singular/eigenmode magnitude estimate for a materialized matrix."""

    rng = np.random.default_rng(7)
    v = rng.normal(size=(A.shape[1],))
    v = v / (np.linalg.norm(v) + 1e-12)
    rho = 0.0
    for _ in range(iters):
        Av = A @ v
        rho = float(np.linalg.norm(Av))
        v = Av / (rho + 1e-12)
    return rho, v

Where to Go Next

Question	Page
Where is the underlying fixed-point mathematics derived?	Mathematical Foundations
Where is a scalar equilibrium checked exactly?	Scalar Equilibrium Example
Which tensor solvers implement the same ideas?	Solvers API