SILVA Networks as Structured Implicit Layers and Vector Attractors via Dynamic Interaction Fields

# SILVA Networks as Structured Implicit Layers and Vector Attractors via Dynamic Interaction Fields
Source: [https://arxiv.org/html/2607.28989](https://arxiv.org/html/2607.28989)
Jose Luis Lima de Jesus Silva\*1,2 1 Federal University of Bahia, Department of Geophysics, Salvador, BA 40170\-115, Brazil 2 Grupo de Estudos e Aplicação de Inteligência Artificial em Geofísica \(GAIA\), Federal University of Bahia, Salvador, BA 40170\-115, Brazil \* Corresponding author:[jseluis\.silva@gmail\.com](https://arxiv.org/html/2607.28989v1/mailto:[email protected])

###### Abstract

Many learning problems require representations that reconcile direct input, nearby structure, and broader context\. In implicit neural layers, these influences are usually absorbed into a single fixed\-point update, making it hard to identify what enters from the stimulus, what propagates locally, what comes from global context, and what is produced by solver dynamics\. Here we introduce SILVA Networks, Structured Implicit Layers and Vector Attractors via Dynamic Interaction Fields\. SILVA separates stimulus, local interaction, global interaction, damping, and readout inside one fixed\-point architecture\. The same template is instantiated for images, molecules, citation networks, and long\-range graph benchmarks through domain\-specific definitions of nodes, neighborhoods, and global summaries\. Experiments and ablations show task\-dependent roles for these terms: local interactions are load\-bearing in the graph tasks, MNIST gains little from recurrence at the tested capacity, and the clearest global benefit appears in a long\-range node\-classification benchmark\. SILVA therefore provides an implicit representation whose internal interaction dynamics can be trained, ablated, visualized, and diagnosed\.

## 1Introduction

Latent representations often have to reconcile several simultaneous sources of influence\. A pixel region, graph node, or molecular atom is driven by its own observed features, by nearby structure, and by broader context from the system in which it appears\[[11](https://arxiv.org/html/2607.28989#bib.bib30),[23](https://arxiv.org/html/2607.28989#bib.bib34)\]\. Modern architectures combine these effects through depth, message passing, attention, or recurrence\[[23](https://arxiv.org/html/2607.28989#bib.bib34),[52](https://arxiv.org/html/2607.28989#bib.bib11)\]\. In an implicit architecture, the same ingredients can be viewed as a driven dynamical system: an input stimulus enters a state space, local and global interactions modify that state, and a numerical solver returns the state reached at convergence\[[6](https://arxiv.org/html/2607.28989#bib.bib5),[24](https://arxiv.org/html/2607.28989#bib.bib8)\]\. The proposed Structured Implicit Layers and Vector Attractors via Dynamic Interaction Fields \(SILVA Networks\) architecture makes this decomposition explicit by separating the stimulus, interaction operators, solver update, and readout state\.

This view continues a long progression from one\-pass transformations toward iterative construction of latent states\. Backpropagation made explicitly layered networks trainable\[[46](https://arxiv.org/html/2607.28989#bib.bib1)\]\. Hopfield networks and energy\-based models showed that a representation can also be the state reached by a dynamical system\[[27](https://arxiv.org/html/2607.28989#bib.bib23),[38](https://arxiv.org/html/2607.28989#bib.bib42)\]\. Residual networks made deep iterative refinement trainable by learning increments around an identity map\[[26](https://arxiv.org/html/2607.28989#bib.bib3)\], while neural ordinary differential equations recast depth as continuous\-time dynamics solved numerically\[[16](https://arxiv.org/html/2607.28989#bib.bib4)\]\. Graph neural networks \(GNNs\) extended iterative state updates to relational data by recursively updating node states over edges\[[48](https://arxiv.org/html/2607.28989#bib.bib2)\]\.

Deep equilibrium models \(DEQs\)\[[6](https://arxiv.org/html/2607.28989#bib.bib5)\]make the fixed\-point view explicit\. Instead of specifying a finite stack of distinct layers, they define a transformation and solve forz∗=fθ​\(z∗,x\)z^\{\\ast\}=f\_\{\\theta\}\(z^\{\\ast\},x\)\. In the original DEQ formulation, implicit differentiation decouples training memory from the number of solver iterations\[[6](https://arxiv.org/html/2607.28989#bib.bib5)\]\. Graph\-specific implicit neural networks extend this idea to node\-state matrices and long\-range graph dependencies\[[24](https://arxiv.org/html/2607.28989#bib.bib8)\], with later operator\-theoretic formulations emphasizing well\-posedness, stability, and accelerated fixed\-point solvers\[[9](https://arxiv.org/html/2607.28989#bib.bib9)\]\. Jacobian\-regularized DEQs further connect the conditioning of the equilibrium Jacobian to forward and backward stability\[[8](https://arxiv.org/html/2607.28989#bib.bib7)\]\. Multiscale DEQs \(MDEQs\) further show that equilibria can couple several feature resolutions inside one joint fixed point\[[7](https://arxiv.org/html/2607.28989#bib.bib6)\]\. The fixed\-point abstraction establishes a powerful representation principle\. SILVA focuses on the internal composition of that fixed\-point update, specifying how direct input drive, local interaction, and system\-wide context enter the same equilibrium computation\.

Several neural architectures already implement interaction mechanisms\. Geometric deep learning treats grids, graphs, and manifolds as one family by generalizing convolution beyond Euclidean arrays\[[15](https://arxiv.org/html/2607.28989#bib.bib35)\]\. Graph convolution, message passing, and interaction networks provide efficient local aggregation over edges or objects\[[11](https://arxiv.org/html/2607.28989#bib.bib30),[34](https://arxiv.org/html/2607.28989#bib.bib18),[23](https://arxiv.org/html/2607.28989#bib.bib34)\]\. Graph attention networks \(GATs\) and transformer attention learn data\-dependent interaction weights\[[53](https://arxiv.org/html/2607.28989#bib.bib10),[52](https://arxiv.org/html/2607.28989#bib.bib11)\]\. Purely local propagation, however, can compress distant signals into fixed\-size node states\[[1](https://arxiv.org/html/2607.28989#bib.bib45)\]\. Dense global attention and latent bottleneck architectures can also entangle local computation, global context, and repeated solver application inside a single learned transformation\[[30](https://arxiv.org/html/2607.28989#bib.bib38),[5](https://arxiv.org/html/2607.28989#bib.bib39)\]\. SILVA separates four computational objects within one fixed\-point layer: the external stimulus, the local interaction operator, the global interaction operator, and the damped solver that accumulates their effects\.

The terminology of stimulus, interaction, stability, and attractor has a long history in physics\-inspired and neurodynamic models\. Hopfield networks framed computation as collective settling in a physical system of interacting units\[[27](https://arxiv.org/html/2607.28989#bib.bib23)\]\. Wilson\-Cowan population dynamics and Amari neural fields model neural activity as a driven interacting population over time\[[56](https://arxiv.org/html/2607.28989#bib.bib24),[2](https://arxiv.org/html/2607.28989#bib.bib25)\]\. Modern neuroscience uses attractor networks to explain persistent activity, cue integration, and error correction in neural populations\[[32](https://arxiv.org/html/2607.28989#bib.bib26)\]\. Visual\-cortical models further show that tuning and perception can depend on recurrent, lateral, predictive, and top\-down interactions rather than feedforward drive alone\[[12](https://arxiv.org/html/2607.28989#bib.bib27),[44](https://arxiv.org/html/2607.28989#bib.bib28),[22](https://arxiv.org/html/2607.28989#bib.bib29)\]\. These works motivate the dynamical vocabulary used here, while SILVA remains a machine\-learning architecture rather than a biological circuit model\.

The local–global tradeoff appears differently across benchmarks\. The ZINC molecular\-property benchmark often rewards chemically local bond structure\[[29](https://arxiv.org/html/2607.28989#bib.bib21)\]\. CLUSTER, from the Benchmarking GNNs suite, is an inductive node\-classification benchmark over synthetic graph instances in which long\-range aggregation is central\[[20](https://arxiv.org/html/2607.28989#bib.bib17)\]\. Cora, Citeseer, and Pubmed citation\-network node classification sit at another boundary\[[49](https://arxiv.org/html/2607.28989#bib.bib41),[58](https://arxiv.org/html/2607.28989#bib.bib13)\]\. Their labelled budgets are small and strong local baselines often suffice, making these datasets a direct test of whether a global term contributes measurable signal once local structure already explains much of the task\.

SILVA Networksnames the mechanism’s three components\.*Structured implicit layers*are fixed\-point layers whose update is split into explicit computational roles\.*Vector attractors*are the settled vector\-valued statesz∗z^\{\\ast\}reached by the solver\.*Dynamic interaction fields*are learned fieldsfθ​\(z,x\)−zf\_\{\\theta\}\(z,x\)\-zthat combine stimulus, local interaction, and global interaction\. At the layer level, SILVA uses a damped updatezk\+1=\(1−α\)​zk\+α​fθ​\(zk,x\)z\_\{k\+1\}=\(1\-\\alpha\)z\_\{k\}\+\\alpha f\_\{\\theta\}\(z\_\{k\},x\), wherefθf\_\{\\theta\}is organized as stimulus plus local and global terms, with domain\-specific adapters defining the state space and operators\.

This formulation generalizes a conventional neural layer in a specific sense\. A standard feedforward transformation computes a representation in one pass\. SILVA instead injects a stimulus and lets it propagate through an interacting state until an approximate vector attractor is reached\. In the locally linearized analysis developed below, a path is an ordered product of self, local, and global operators carrying stimulus to a state coordinate\. These stimulus\-to\-state paths make precise how direct input, local neighborhoods, global context, damping, and finite solver depth combine\. The same outer fixed\-point solver is used across domains, while each domain defines what counts as a node, a local neighborhood, and a global summary\. The shared object is an inspectable interaction template instantiated through domain\-specific state spaces and operators\.

The term equilibrium is used in two precise senses\. Architecturally, all models run a damped fixed\-point solver and read out the terminal statezKz\_\{K\}as an approximate fixed point\. In the strict DEQ training sense, the adjoint\-trained MNIST diagnostic uses a GMRES adjoint solve\[[47](https://arxiv.org/html/2607.28989#bib.bib37)\], whereas the remaining experiments use truncated backpropagation through the unrolled fixed\-point solver\.

The empirical results are informative because the global term is not uniformly beneficial\. On MNIST\[[37](https://arxiv.org/html/2607.28989#bib.bib52)\]and citation\-network node classification, the global term does not automatically improve accuracy over local or no\-interaction variants, so the added operators are tested directly rather than assumed to help\. On the long\-range CLUSTER benchmark, where broader aggregation is expected to matter, the full local\+global SILVA model reaches73\.04±0\.60%73\.04\\pm 0\.60\\%validation accuracy across five seeds\. This exceeds the best matched local graph baseline by14\.2914\.29percentage points and exceeds a no\-global SILVA ablation by5\.495\.49percentage points\. The same experiments also expose solver residuals, spectral\-radius behavior, learned interaction structure, and equilibrium\-state geometry, providing measurements of the dynamics as well as task scores\.

Contributions\.

1. 1\.SILVA Networks define structured implicit layers in which a stimulus term, a local interaction term, a global interaction term, and a damped solver appear as separate parts of one dynamic interaction field \(Section[2](https://arxiv.org/html/2607.28989#S2)\)\.
2. 2\.The same template is instantiated for vision, molecular regression, citation\-network node classification, and long\-range graph benchmarks by specifying the state nodes, local graph, global summary, solver, and backward mode used in each domain\.
3. 3\.Stacked SILVA models compose separately solved approximate vector attractor states with distinct parameters and timescales, rather than increasing the iteration count of a single equilibrium solve\.
4. 4\.A finite path\-sum interpretation of the locally linearized solver connects solver depth, damping, spectral radius, and effective stimulus\-to\-state influence through the same Jacobian\-level operator \(Sections[2\.10](https://arxiv.org/html/2607.28989#S2.SS10)and[2\.11](https://arxiv.org/html/2607.28989#S2.SS11)\)\.
5. 5\.The architecture is evaluated across vision, chemistry, citation\-network, and long\-range graph settings, including local/global ablations that identify when the global interaction contributes measurable signal and when local structure already dominates\.
6. 6\.The empirical characterization includes seed sensitivity, finite solver budgets, truncated backpropagation, stability diagnostics, and direct benchmark overlap with general\-purpose multi\-domain architectures\.

## 2Method

### 2\.1The equilibrium layer

Letz∈ℝN×dz\\in\\mathbb\{R\}^\{N\\times d\}be a layer’s state overNNentities \(pixels/units for vision, atoms for chemistry, nodes for citation graphs\) with hidden widthdd, andxxthe layer’s input embedding\. All experiments share one outer template\. A per\-domain update mapfθ​\(z,x\)f\_\{\\theta\}\(z,x\)reaches a fixed pointz∗=fθ​\(z∗,x\)z^\{\\ast\}=f\_\{\\theta\}\(z^\{\\ast\},x\)via damped fixed\-point \(Picard\) iteration\[[10](https://arxiv.org/html/2607.28989#bib.bib22)\]fork=0,…,K−1k=0,\\dots,K\-1,

zk\+1=\(1−α\)​zk\+α⋅fθ​\(zk,x\),z\_\{k\+1\}=\(1\-\\alpha\)\\,z\_\{k\}\+\\alpha\\cdot f\_\{\\theta\}\(z\_\{k\},x\),\(1\)Herezkz\_\{k\}is the solver iterate afterkkupdates,KKis the total iteration budget,α∈\(0,1\]\\alpha\\in\(0,1\]is the damping coefficient,θ\\thetadenotes all learned parameters in the update map, andzKz\_\{K\}is the approximate fixed point used by the downstream readout\. The inputxxis already encoded by the domain\-specific input encoder defined below\. The update map always combines stimulus, local, and global terms\. Its exact composition varies by domain, so the equations below state the evaluated compositions separately\. The local and global maps act on a recurrent signalz~=χ𝒟​\(z\)\\tilde\{z\}=\\chi\_\{\\mathcal\{D\}\}\(z\)rather than necessarily on the raw solver state, where𝒟\\mathcal\{D\}indexes the domain or model family\. For ZINC,χ𝒟\\chi\_\{\\mathcal\{D\}\}is the identity\. For the Cora, Citeseer, and Pubmed citation networks\[[49](https://arxiv.org/html/2607.28989#bib.bib41),[58](https://arxiv.org/html/2607.28989#bib.bib13)\], CLUSTER\[[20](https://arxiv.org/html/2607.28989#bib.bib17)\], MNIST\[[37](https://arxiv.org/html/2607.28989#bib.bib52)\], and CIFAR\-10\[[36](https://arxiv.org/html/2607.28989#bib.bib57)\],χ𝒟​\(z\)=tanh⁡\(z\)\\chi\_\{\\mathcal\{D\}\}\(z\)=\\tanh\(z\), as in the evaluated models\.

The chemistry experiments use the ZINC molecular\-property benchmark\[[29](https://arxiv.org/html/2607.28989#bib.bib21)\]\. The node\-classification experiments use the Cora, Citeseer, and Pubmed citation graphs\[[49](https://arxiv.org/html/2607.28989#bib.bib41),[58](https://arxiv.org/html/2607.28989#bib.bib13)\]\. The long\-range node\-classification experiments use the CLUSTER benchmark\[[20](https://arxiv.org/html/2607.28989#bib.bib17)\]\. These non\-vision experiments use

z~\\displaystyle\\tilde\{z\}=χ𝒟​\(z\),\\displaystyle=\\chi\_\{\\mathcal\{D\}\}\(z\),\(2\)fθ​\(z,x\)\\displaystyle f\_\{\\theta\}\(z,x\)=LayerNorm​\(ReLU​\(Wstim​x\+L​\(z~\)\+G​\(z~\)\)\)\.\\displaystyle=\\mathrm\{LayerNorm\}\\Big\(\\mathrm\{ReLU\}\\big\(W\_\{\\mathrm\{stim\}\}x\+L\(\\tilde\{z\}\)\+G\(\\tilde\{z\}\)\\big\)\\Big\)\.\(3\)HereWstim​xW\_\{\\mathrm\{stim\}\}x,L​\(z~\)L\(\\tilde\{z\}\), andG​\(z~\)G\(\\tilde\{z\}\)all have the same shape aszz\. The matrixWstimW\_\{\\mathrm\{stim\}\}maps the encoded input into the state space\. The normalization is standard layer normalization\[[4](https://arxiv.org/html/2607.28989#bib.bib47)\]\. The mapLLis Graph Attention Network \(GAT\)\-style attention\[[53](https://arxiv.org/html/2607.28989#bib.bib10)\]\(Section[2\.3](https://arxiv.org/html/2607.28989#S2.SS3)\), andGGis a mean\-field broadcast over graph nodes \(Section[2\.6](https://arxiv.org/html/2607.28989#S2.SS6)\)\. The vision experiments use MNIST\[[37](https://arxiv.org/html/2607.28989#bib.bib52)\]and CIFAR\-10\[[36](https://arxiv.org/html/2607.28989#bib.bib57)\]\. They use a distinct, simpler composition without the outer LayerNorm/ReLU,

z~\\displaystyle\\tilde\{z\}=tanh⁡\(z\),\\displaystyle=\\tanh\(z\),\(4\)fθ​\(z,x\)\\displaystyle f\_\{\\theta\}\(z,x\)=Wstim​x\+Lvis​\(z~\)\+Gvis​\(z~\),\\displaystyle=W\_\{\\mathrm\{stim\}\}x\+L\_\{\\mathrm\{vis\}\}\(\\tilde\{z\}\)\+G\_\{\\mathrm\{vis\}\}\(\\tilde\{z\}\),\(5\)whereLvisL\_\{\\mathrm\{vis\}\}andGvisG\_\{\\mathrm\{vis\}\}are vision\-specific rather than instances ofLL/GGabove,

Lvis​\(z~\)i\\displaystyle L\_\{\\mathrm\{vis\}\}\(\\tilde\{z\}\)\_\{i\}=1\|𝒩​\(i\)\|​∑j∈𝒩​\(i\)WL​z~j,\\displaystyle=\\frac\{1\}\{\|\\mathcal\{N\}\(i\)\|\}\\sum\_\{j\\in\\mathcal\{N\}\(i\)\}W\_\{L\}\\tilde\{z\}\_\{j\},\(6\)Gvis​\(z~\)\\displaystyle G\_\{\\mathrm\{vis\}\}\(\\tilde\{z\}\)=z~​softmax​\(\(Wq​z~\)​\(Wk​z~\)⊤/d\)\.\\displaystyle=\\tilde\{z\}\\,\\mathrm\{softmax\}\\\!\\big\(\(W\_\{q\}\\tilde\{z\}\)\(W\_\{k\}\\tilde\{z\}\)^\{\\top\}/\\sqrt\{d\}\\big\)\.\(7\)In Eq\. \([6](https://arxiv.org/html/2607.28989#S2.E6)\),iiandjjindex hidden\-channel nodes within one sample,\|𝒩​\(i\)\|\|\\mathcal\{N\}\(i\)\|is the number of selected neighbors, andWLW\_\{L\}is shared across all such neighbors\. In Eq\. \([7](https://arxiv.org/html/2607.28989#S2.E7)\),WqW\_\{q\}andWkW\_\{k\}are learned query and key projections applied toz~\\tilde\{z\}for the channel self\-attention map\.LvisL\_\{\\mathrm\{vis\}\}is a degree\-normalized average over akk\-nearest\-neighbor graph in*channel space*, meaning that𝒩​\(i\)\\mathcal\{N\}\(i\)ranges over theddhidden channels of one sample, not over other samples or graph nodes\. It uses one shared linear mapWLW\_\{L\}and no per\-edge learned attention weights\.GvisG\_\{\\mathrm\{vis\}\}is genuine self\-attention, but over that one sample’s ownddchannels, not a mean\-field vector pooled across nodes or samples\. Both compositions are damped by the same Eq\. \([1](https://arxiv.org/html/2607.28989#S2.E1)\), so the outer fixed\-point template is shared across every domain\. Because the exact local and global maps vary, the theoretical results below \(spectral radius, energy diagnostic, local descent surrogate, and path\-sum\) are stated generically in terms offθf\_\{\\theta\}and its Jacobian\. They apply to either composition\. The graph\-attention and mean\-field\-broadcast descriptions used in the rest of this section refer specifically to Eq\. \([3](https://arxiv.org/html/2607.28989#S2.E3)\), not Eq\. \([5](https://arxiv.org/html/2607.28989#S2.E5)\)\.

In the equations below, the coordinateziz\_\{i\}denotes the state of node, atom, channel, or hidden unitii, and an interaction coefficient indexed by\(i,j\)\(i,j\)describes influence from sourcejjto receiverii\. The stimulus issθ​\(x\)=Wstim​xs\_\{\\theta\}\(x\)=W\_\{\\mathrm\{stim\}\}x, with any bias absorbed into the encoded input or intosθs\_\{\\theta\}\. The symbolz~\\tilde\{z\}denotes the recurrent signal consumed by local and global maps\. The symbolsLθL\_\{\\theta\}andGθG\_\{\\theta\}denote the local and global update maps executed inside the solver\. Their local linearizations at the evaluated state areMLM\_\{L\}andMGM\_\{G\}, andM=ML\+MGM=M\_\{L\}\+M\_\{G\}denotes the one\-step interaction operator used only in the linearized analysis\. The damped step operatorsPselfP\_\{\\mathrm\{self\}\},PlocalP\_\{\\mathrm\{local\}\}, andPglobalP\_\{\\mathrm\{global\}\}, the one\-step propagatorΠi←j\\Pi\_\{i\\leftarrow j\}, and the response operatorRK,αR\_\{K,\\alpha\}are derived objects, not additional learned parameters\. Indicesi,ji,jdenote nodes or state coordinates,k,tk,tdenote solver iteration, andrrdenotes path length or matrix power\. Batched computations may store states as row tensors, but the displayed equations use this receiver\-first coordinate convention\.

###### Definition 1\(SILVA layer\)\.

A SILVA layer is a damped implicit layer with statez∈ℝN×dz\\in\\mathbb\{R\}^\{N\\times d\}, stimulussθ​\(x\)s\_\{\\theta\}\(x\), domain\-specific interaction mapsLθL\_\{\\theta\}andGθG\_\{\\theta\}, and dynamic interaction field

ℱθ​\(z,x\)\\displaystyle\\mathcal\{F\}\_\{\\theta\}\(z,x\)=fθ​\(z,x\)−z,\\displaystyle=f\_\{\\theta\}\(z,x\)\-z,\(8\)zk\+1\\displaystyle z\_\{k\+1\}=zk\+α​ℱθ​\(zk,x\)\.\\displaystyle=z\_\{k\}\+\\alpha\\,\\mathcal\{F\}\_\{\\theta\}\(z\_\{k\},x\)\.\(9\)A*vector attractor*is a locally stable zero of this field,ℱθ​\(z∗,x\)=0\\mathcal\{F\}\_\{\\theta\}\(z^\{\\ast\},x\)=0, equivalentlyz∗=fθ​\(z∗,x\)z^\{\\ast\}=f\_\{\\theta\}\(z^\{\\ast\},x\)\. Under the local linearization used for the path\-sum diagnostics, writings=Wstim​xs=W\_\{\\mathrm\{stim\}\}xandM=ML\+MGM=M\_\{L\}\+M\_\{G\}gives

zk\+1≈α​s\+\(\(1−α\)​I\+α​ML\+α​MG\)​zk\.z\_\{k\+1\}\\approx\\alpha s\+\\big\(\(1\-\\alpha\)I\+\\alpha M\_\{L\}\+\\alpha M\_\{G\}\\big\)z\_\{k\}\.\(10\)HereIIis the identity operator on the flattenedN​dNd\-dimensional state,s∈ℝN×ds\\in\\mathbb\{R\}^\{N\\times d\}is the stimulus reshaped as a state\-sized tensor, andML,MG∈ℝN​d×N​dM\_\{L\},M\_\{G\}\\in\\mathbb\{R\}^\{Nd\\times Nd\}are the local and global Jacobian\-level operators at the evaluated state\. The finite solver therefore propagates an external stimulus through three operator types\. Self\-persistence contributes\(1−α\)​I\(1\-\\alpha\)I, local interactions contributeα​ML\\alpha M\_\{L\}, and global interactions contributeα​MG\\alpha M\_\{G\}\. The layer uses damped identity persistence for the self term, with local and global maps providing the learned or state\-dependent interactions\. Figure[1](https://arxiv.org/html/2607.28989#S2.F1)a gives the graphical reading of Eqs\. \([1](https://arxiv.org/html/2607.28989#S2.E1)\)–\([9](https://arxiv.org/html/2607.28989#S2.E9)\)\. The stimulus arrow issθ​\(x\)s\_\{\\theta\}\(x\), the local and global arrows areLθ​\(z~k\)L\_\{\\theta\}\(\\tilde\{z\}\_\{k\}\)andGθ​\(z~k\)G\_\{\\theta\}\(\\tilde\{z\}\_\{k\}\), and the gray arrow is the damped self\-persistence contribution\(1−α\)​zk\(1\-\\alpha\)z\_\{k\}\. The combination node represents the formation of the next solver iteratezk\+1z\_\{k\+1\}\. Figure[1](https://arxiv.org/html/2607.28989#S2.F1)b unfolds the same update across solver time\. Colored stimulus\-to\-state paths visualize the finite powers in Eq\. \([10](https://arxiv.org/html/2607.28989#S2.E10)\), and the final state is the approximate vector attractor read out by the task head\. The diagram links stimulus, interaction operators, damping, finite paths, and attractor readout using the same notation as the equations\.

![Refer to caption](https://arxiv.org/html/2607.28989v1/x1.png)Figure 1:Conceptual mechanism of a SILVA equilibrium layer\. \(a\) Stimulus, local interaction, global context, and self\-persistence in the damped update of Eqs\. \([1](https://arxiv.org/html/2607.28989#S2.E1)\)–\([9](https://arxiv.org/html/2607.28989#S2.E9)\)\. \(b\) Finite stimulus\-to\-state paths and vector\-attractor readout associated with Eq\. \([10](https://arxiv.org/html/2607.28989#S2.E10)\)\.#### Linear equilibrium limit\.

The simplest limit of the SILVA update makes the stimulus\-response structure explicit\. Letb∈ℝNb\\in\\mathbb\{R\}^\{N\}be a state vector, letm=sθ​\(x\)=Wstim​x\+cm=s\_\{\\theta\}\(x\)=W\_\{\\mathrm\{stim\}\}x\+cbe the stimulus, and letA∈ℝN×NA\\in\\mathbb\{R\}^\{N\\times N\}be a receiver\-first interaction matrix whose entryAi​jA\_\{ij\}carries influence from source coordinatejjto receiver coordinateii\. The linear equilibrium condition is

Rearranging gives\(I−A\)​b=m\(I\-A\)b=mand therefore

b=\(I−A\)−1​m,b=\(I\-A\)^\{\-1\}m,\(12\)wheneverI−AI\-Ais invertible\. Ifρ​\(A\)<1\\rho\(A\)<1, the inverse is the Neumann series\(I−A\)−1=∑t=0∞At\(I\-A\)^\{\-1\}=\\sum\_\{t=0\}^\{\\infty\}A^\{t\}, so the hidden state is the sum of the direct stimulus plus all repeated intra\-layer interaction paths\. SILVA extends this linear limit by replacing the single matrixAAwith the locally linearized Jacobian of a domain\-specific nonlinear update map whose terms are stimulus, local interaction, and global interaction, with the nonlinearities and normalization stated in Eqs\. \([3](https://arxiv.org/html/2607.28989#S2.E3)\)–\([5](https://arxiv.org/html/2607.28989#S2.E5)\)\. The state is reached by the damped solver of Eq\. \([1](https://arxiv.org/html/2607.28989#S2.E1)\) rather than by forming an explicit inverse\.

By the Banach fixed\-point theorem\[[10](https://arxiv.org/html/2607.28989#bib.bib22)\], a unique fixed point exists and damped iteration converges wheneverfθ​\(⋅,x\)f\_\{\\theta\}\(\\cdot,x\)is a contraction\. In the locally linearized analysis used for the diagnostic, the corresponding fixed\-point stability condition is

ρ​\(Jf​\(z∗\)\)\\displaystyle\\rho\\big\(J\_\{f\}\(z^\{\\ast\}\)\\big\)<1,\\displaystyle<1,\(13\)Jf​\(z∗\)\\displaystyle J\_\{f\}\(z^\{\\ast\}\)=∂fθ/∂z\|z∗,\\displaystyle=\\partial f\_\{\\theta\}/\\partial z\\,\\big\|\_\{z^\{\\ast\}\},\(14\)whereρ\\rhois the spectral radius of the Jacobian at the fixed point\. The JacobianJf​\(z∗\)∈ℝN​d×N​dJ\_\{f\}\(z^\{\\ast\}\)\\in\\mathbb\{R\}^\{Nd\\times Nd\}is taken with respect to the flattened state coordinates ofzz, withxxheld fixed\. The damped solver itself has local iteration Jacobian

J𝒯α​\(z∗\)=\(1−α\)​I\+α​Jf​\(z∗\)\.J\_\{\\mathcal\{T\}\_\{\\alpha\}\}\(z^\{\\ast\}\)=\(1\-\\alpha\)I\+\\alpha J\_\{f\}\(z^\{\\ast\}\)\.\(15\)The spectral\-radius diagnostic estimatesρ𝒯=ρ​\(J𝒯α\)\\rho\_\{\\mathcal\{T\}\}=\\rho\(J\_\{\\mathcal\{T\}\_\{\\alpha\}\}\), the spectral radius of the actual damped update used by the finite solver\. This diagnostic coincides withρ​\(Jf\)\\rho\(J\_\{f\}\)whenα=1\\alpha=1, but otherwise measures the local spectral stability of the executed iteration rather than the undamped update alone\. Stacked SILVA architectures use multiple equilibria\. For a stack ofLLSILVA layers, layerℓ\\ellhas its own parametersθℓ\\theta\_\{\\ell\}, statezℓz\_\{\\ell\}, damping rateαℓ\\alpha\_\{\\ell\}, solver budgetKℓK\_\{\\ell\}, update mapfθℓ\(ℓ\)f\_\{\\theta\_\{\\ell\}\}^\{\(\\ell\)\}, and connectorψℓ\\psi\_\{\\ell\}to the next layer\.

zℓ,0\\displaystyle z\_\{\\ell,0\}=0,\\displaystyle=0,\(16\)zℓ,k\+1\\displaystyle z\_\{\\ell,k\+1\}=\(1−αℓ\)​zℓ,k\+αℓ​fθℓ\(ℓ\)​\(zℓ,k,xℓ\),\\displaystyle=\(1\-\\alpha\_\{\\ell\}\)z\_\{\\ell,k\}\+\\alpha\_\{\\ell\}f\_\{\\theta\_\{\\ell\}\}^\{\(\\ell\)\}\(z\_\{\\ell,k\},x\_\{\\ell\}\),\(17\)xℓ\+1\\displaystyle x\_\{\\ell\+1\}=ψℓ​\(zℓ,Kℓ\)\.\\displaystyle=\\psi\_\{\\ell\}\(z\_\{\\ell,K\_\{\\ell\}\}\)\.\(18\)Hereℓ\\ellindexes equilibrium layers,zℓ,kz\_\{\\ell,k\}is thekkth solver iterate inside layerℓ\\ell,xℓx\_\{\\ell\}is the input supplied to that layer, andψℓ\\psi\_\{\\ell\}is a learned or fixed connector that maps the final state of layerℓ\\ellinto the input space of layerℓ\+1\\ell\+1\. Stacking is therefore a composition of separately solved vector attractors, not an increase in Picard steps inside one solver\. The first attractor receives the encoded data\. Later attractors receive transformed states from earlier attractors and solve their own local–global interaction fields\. This construction is related to multiscale DEQs\[[7](https://arxiv.org/html/2607.28989#bib.bib6)\], but differs in the equilibrium that is solved\. MDEQ solves one enlarged equilibriumZ∗=FΘ​\(Z∗,x\)Z^\{\\ast\}=F\_\{\\Theta\}\(Z^\{\\ast\},x\)in which several feature resolutions are coupled and driven to equilibrium simultaneously\. The sequential SILVA stack solves each stage as its own root problem, with distinct parameters and a distinct timescale, before its vector attractor is passed to the next stage\. In a local linearization, MDEQ has one coupled response operator over the joint multiscale state, whereas a SILVA stack composes stagewise responses,Rstack≈RL​D​ψL−1​⋯​R2​D​ψ1​R1R\_\{\\mathrm\{stack\}\}\\approx R\_\{L\}D\\psi\_\{L\-1\}\\cdots R\_\{2\}D\\psi\_\{1\}R\_\{1\}, withRℓR\_\{\\ell\}the finite\- or infinite\-iteration response of layerℓ\\ellandD​ψℓD\\psi\_\{\\ell\}the Jacobian of the connectorψℓ\\psi\_\{\\ell\}at the passed state\. Interaction histories are first accumulated within one stimulus\-driven attractor and then handed to the next, instead of all scales being solved as one simultaneous equilibrium\. In the main two\-layer stack used for CIFAR\-10, ZINC, citation networks, and CLUSTER, Layer 1 \(α1=0\.5\\alpha\_\{1\}=0\.5, fast\) feedstanh⁡\(z1∗\)\\tanh\(z\_\{1\}^\{\\ast\}\)to Layer 2 \(α2=0\.2\\alpha\_\{2\}=0\.2, slow\)\. CIFAR\-10 and ZINC useK1=K2=20K\_\{1\}=K\_\{2\}=20Picard iterations\. The node\-classification and graph\-benchmark experiments useKℓ=15K\_\{\\ell\}=15unless a specific ablation changes the stack\. The MNIST dynamic\-channel\-graph experiments use a single equilibrium layer withα=0\.25\\alpha=0\.25andK=20K=20\. More generally, any number of such equilibrium layers can be stacked, each at its ownαℓ\\alpha\_\{\\ell\}, giving a hierarchy of timescales analogous in spirit to clockwork and hierarchical recurrent architectures\[[35](https://arxiv.org/html/2607.28989#bib.bib14),[21](https://arxiv.org/html/2607.28989#bib.bib15)\]\. Here, every rung is a full equilibrium solve over the interaction structure below, not a fixed\-depth recurrent step\. The stacked\-equilibrium result below evaluates a three\-rung instance\.

#### Relation to classical attractor and interaction models\.

Several established models provide useful reference points for this formulation\. Hopfield networks introduced recurrent collective computation as convergence to associative memory states, usually under symmetry assumptions that support an energy interpretation\[[27](https://arxiv.org/html/2607.28989#bib.bib23)\]\. Wilson\-Cowan and Amari models describe population dynamics and continuous neural fields with biologically interpretable activity variables and connection kernels\[[56](https://arxiv.org/html/2607.28989#bib.bib24),[2](https://arxiv.org/html/2607.28989#bib.bib25)\]\. Predictive\-coding and feedback accounts of visual processing emphasize interactions between sensory drive and contextual or top\-down signals\[[44](https://arxiv.org/html/2607.28989#bib.bib28),[22](https://arxiv.org/html/2607.28989#bib.bib29)\]\. Interaction networks and message\-passing networks learn relations among objects or graph nodes\[[11](https://arxiv.org/html/2607.28989#bib.bib30),[23](https://arxiv.org/html/2607.28989#bib.bib34)\]\. These lines of work motivate the language of recurrent interaction, but the model studied here is a learned implicit neural layer rather than a biological circuit or a symmetric associative memory\. Its attractor is the fixed point of the learned update map reached by damped numerical iteration\. Its local and global operators are domain\-specific neural modules trained from data\. The additional structure is the explicit separation of stimulus, self\-persistence, local interactions, and global interactions inside the solver, which makes the resulting finite path\-response analysis possible\.

#### Evaluated model families\.

The same damped outer iteration in Eq\. \([1](https://arxiv.org/html/2607.28989#S2.E1)\) is used throughout, but the state space, local graph, global term, nonlinearity, and backward mode are domain\-specific\. In the vision experiments on MNIST and CIFAR\-10, state nodes are the hidden channels of one sample\. The local graph is a dynamic channelkk\-nearest\-neighbor graph\. The local term is a degree\-normalized neighbor average through a shared matrixWLW\_\{L\}, the global term is per\-sample channel self\-attention, and both consumetanh⁡\(z\)\\tanh\(z\)\. The update map is a raw sum\. The finite\-solver vision sweeps use Picard iteration with truncated backpropagation\. Anderson acceleration\[[3](https://arxiv.org/html/2607.28989#bib.bib20),[54](https://arxiv.org/html/2607.28989#bib.bib36)\]appears only in solver comparisons, and the separate adjoint\-trained MNIST diagnostic uses an undamped two\-layer interacting equilibrium model with a hidden\-channel nearest\-neighbor attention local branch, a learned static global interaction matrix, and a generalized minimal residual \(GMRES\) adjoint\[[47](https://arxiv.org/html/2607.28989#bib.bib37)\]\. This diagnostic is used for the hidden\-state geometry and stability analyses rather than as the principal MNIST accuracy model\.

For citation networks, state nodes are graph nodes and the local graph is the observed citation edge set\. The local term is edge\-based graph attention ontanh⁡\(z\)\\tanh\(z\)\. The global term is a scalar\-gated mean\-field broadcast ontanh⁡\(z\)\\tanh\(z\), with bounded top\-kkattention used only in the optional global variant\. For CLUSTER, state nodes are benchmark graph nodes, the local graph is the benchmark edge set, the local term is graph attention ontanh⁡\(z\)\\tanh\(z\), and the global term is a mean\-field broadcast ontanh⁡\(z\)\\tanh\(z\)unless it is explicitly ablated\. These graph families use the LayerNorm\-ReLU update block, Picard iteration, and truncated backpropagation through the unrolled solver\. For ZINC, state nodes are atoms, the local graph is the molecular bond graph, and the enhanced model includes bond attributes\. Both local and global terms consume rawzz\. The local term is graph attention or bond\-aware graph attention, while the global term is graph mean pooling followed by a learned broadcast\. ZINC also uses the LayerNorm\-ReLU update block, Picard iteration, and truncated backpropagation\.

Algorithm 1Forward pass for one damped local–global equilibrium layer1:Encode raw input as

xxand initialize

z0z\_\{0\}\(zeros or a learned/input\-derived initial state\)

2:for

k=0,…,K−1k=0,\\dots,K\-1do

3:Compute the stimulus

s←Wstim​xs\\leftarrow W\_\{\\mathrm\{stim\}\}x
4:Compute the recurrent signal

z~k←χ𝒟​\(zk\)\\tilde\{z\}\_\{k\}\\leftarrow\\chi\_\{\\mathcal\{D\}\}\(z\_\{k\}\)
5:Build or read the local graph

𝒩\(k\)\\mathcal\{N\}^\{\(k\)\}for this domain

6:Compute

hlocal←L​\(z~k,𝒩\(k\)\)h\_\{\\mathrm\{local\}\}\\leftarrow L\(\\tilde\{z\}\_\{k\},\\mathcal\{N\}^\{\(k\)\}\)
7:Compute

hglobal←G​\(z~k\)h\_\{\\mathrm\{global\}\}\\leftarrow G\(\\tilde\{z\}\_\{k\}\)
8:Form

fθ​\(zk,x\)f\_\{\\theta\}\(z\_\{k\},x\)using the domain\-specific composition above

9:Update

zk\+1←\(1−α\)​zk\+α​fθ​\(zk,x\)z\_\{k\+1\}\\leftarrow\(1\-\\alpha\)z\_\{k\}\+\\alpha f\_\{\\theta\}\(z\_\{k\},x\), or replace this step by the Anderson extrapolation of Eq\. \([52](https://arxiv.org/html/2607.28989#S2.E52)\)

10:endfor

11:Return

zKz\_\{K\}as the approximate fixed point and apply the task readout

Wout​zKW\_\{\\mathrm\{out\}\}z\_\{K\}

Algorithm[1](https://arxiv.org/html/2607.28989#alg1)uses the same notation as Eq\. \([1](https://arxiv.org/html/2607.28989#S2.E1)\)\. The neighborhood object𝒩\(k\)\\mathcal\{N\}^\{\(k\)\}is either a fixed graph, as in molecules and citation networks, or a state\-dependent graph, as in the vision channelkk\-nearest\-neighbor construction\. The final matrixWoutW\_\{\\mathrm\{out\}\}denotes the task readout introduced in the Training subsection\.

### 2\.2Interaction\-operator decomposition

Linearizing Eq\. \([1](https://arxiv.org/html/2607.28989#S2.E1)\) at fixed attention/gating weights gives a block\-structured operatorΩ\\Omegaon the flattened state, decomposed into three mechanistically distinct components,

Ω=Ωself\+Ωlocal\+Ωglobal,\\Omega=\\Omega\_\{\\mathrm\{self\}\}\+\\Omega\_\{\\mathrm\{local\}\}\+\\Omega\_\{\\mathrm\{global\}\},\(19\)HereΩ\\Omegais the Jacobian\-level interaction operator, not a separate stored weight matrix\. The three terms areΩself\\Omega\_\{\\mathrm\{self\}\}\(diagonal, node\-wise persistence\),Ωlocal\\Omega\_\{\\mathrm\{local\}\}\(sparse, neighborhood message passing\), andΩglobal\\Omega\_\{\\mathrm\{global\}\}\(dense low\-rank, mean\-field context\)\. Node persistence is provided by the\(1−α\)​zk\(1\-\\alpha\)z\_\{k\}term of the damped iteration \(Eq\. \([1](https://arxiv.org/html/2607.28989#S2.E1)\)\) rather than a separate learned self\-weight, so the model comprises two directly\-parameterized interaction components \(local, global\) plus this implicit self\-term\.

#### Solver\-level operator interpretation\.

These operators are more than diagnostic labels\. They correspond to the neural computations executed at each solver step\. In the locally linearized form used for the diagnostics, one damped step can be written

zk\+1\\displaystyle z\_\{k\+1\}=α​s\+\(Pself\+Plocal\+Pglobal\)​zk,\\displaystyle=\\alpha s\+\\big\(P\_\{\\mathrm\{self\}\}\+P\_\{\\mathrm\{local\}\}\+P\_\{\\mathrm\{global\}\}\\big\)z\_\{k\},\(20\)Pself\\displaystyle P\_\{\\mathrm\{self\}\}=\(1−α\)​I,\\displaystyle=\(1\-\\alpha\)I,\(21\)Plocal\\displaystyle P\_\{\\mathrm\{local\}\}=α​ML,\\displaystyle=\\alpha M\_\{L\},\(22\)Pglobal\\displaystyle P\_\{\\mathrm\{global\}\}=α​MG,\\displaystyle=\\alpha M\_\{G\},\(23\)wheres=Wstim​xs=W\_\{\\mathrm\{stim\}\}xis the injected stimulus\. The matricesMLM\_\{L\}andMGM\_\{G\}are the local linearizations of the evaluated local and global branches\. They include the mapz↦z~z\\mapsto\\tilde\{z\}and any active normalization/nonlinearity factors\. The operatorsPselfP\_\{\\mathrm\{self\}\},PlocalP\_\{\\mathrm\{local\}\}, andPglobalP\_\{\\mathrm\{global\}\}all act on the same flattened state space asMLM\_\{L\}andMGM\_\{G\}\. Starting fromz0=0z\_\{0\}=0,

zK\\displaystyle z\_\{K\}=α​∑t=0K−1\(Pself\+Plocal\+Pglobal\)t​s,\\displaystyle=\\alpha\\sum\_\{t=0\}^\{K\-1\}\\big\(P\_\{\\mathrm\{self\}\}\+P\_\{\\mathrm\{local\}\}\+P\_\{\\mathrm\{global\}\}\\big\)^\{t\}s,\(24\)zK\\displaystyle z\_\{K\}=α​∑t=0K−1∑w∈\{self,local,global\}tPwt​⋯​Pw1​s\.\\displaystyle=\\alpha\\sum\_\{t=0\}^\{K\-1\}\\sum\_\{w\\in\\\{\\mathrm\{self\},\\mathrm\{local\},\\mathrm\{global\}\\\}^\{t\}\}P\_\{w\_\{t\}\}\\cdots P\_\{w\_\{1\}\}s\.\(25\)Each wordw=\(w1,…,wt\)w=\(w\_\{1\},\\ldots,w\_\{t\}\)is a concrete stimulus\-to\-state interaction history through the evaluated layer\. A self symbol keeps information through damping, a local symbol applies the neighborhood update, and a global symbol applies the broadcast or attention\-mediated context update\. The coordinate\(Pwt​⋯​Pw1​s\)i\(P\_\{w\_\{t\}\}\\cdots P\_\{w\_\{1\}\}s\)\_\{i\}is the contribution of that particular interaction history to node or channelii\. The later Neumann\-series path\-sum form \(Section[2\.10](https://arxiv.org/html/2607.28989#S2.SS10)\) is the same expansion after collecting the self\-persistence factors and writing the remaining interaction operator asM=ML\+MGM=M\_\{L\}\+M\_\{G\}\. Fort=0t=0, the word is empty and the product is the identity, giving the direct stimulus contribution\.

Written per\-channel for the undamped mapfθf\_\{\\theta\}, the update for feature channelccof nodeiiis

fθ​\(z\(k\),x\)i,c=σ​\(ui,c⏟stimulus\+∑j∈𝒩​\(i\)αi​j\(k\)​∑fWc​flocal​z~j,f\(k\)⏟local,​Ωlocal\+1N​∑j∑fWc​fglobal​z~j,f\(k\)⏟global,​Ωglobal\)\.f\_\{\\theta\}\(z^\{\(k\)\},x\)\_\{i,c\}=\\sigma\\Bigg\(\\underbrace\{u\_\{i,c\}\}\_\{\\text\{stimulus\}\}\+\\underbrace\{\\sum\_\{j\\in\\mathcal\{N\}\(i\)\}\\alpha\_\{ij\}^\{\(k\)\}\\sum\_\{f\}W\_\{cf\}^\{\\mathrm\{local\}\}\\tilde\{z\}\_\{j,f\}^\{\(k\)\}\}\_\{\\text\{local, \}\\Omega\_\{\\mathrm\{local\}\}\}\+\\underbrace\{\\tfrac\{1\}\{N\}\\sum\_\{j\}\\sum\_\{f\}W\_\{cf\}^\{\\mathrm\{global\}\}\\tilde\{z\}\_\{j,f\}^\{\(k\)\}\}\_\{\\text\{global, \}\\Omega\_\{\\mathrm\{global\}\}\}\\Bigg\)\.\(26\)In Eq\. \([26](https://arxiv.org/html/2607.28989#S2.E26)\),u=sθ​\(x\)u=s\_\{\\theta\}\(x\)is the stimulus,iiis the receiver node,jjis a source node,ccis the receiver feature channel, andffis the source feature channel being summed through the channel\-mixing matrices\. The coefficientαi​j\(k\)\\alpha\_\{ij\}^\{\(k\)\}is either a learned attention weight or a normalized edge weight at solver iterationkk, depending on the domain\-specific local operator\. The recurrent signal isz~\(k\)=χ𝒟​\(z\(k\)\)\\tilde\{z\}^\{\(k\)\}=\\chi\_\{\\mathcal\{D\}\}\(z^\{\(k\)\}\), andσ\\sigmadenotes the domain’s outer normalization/nonlinearity block or the identity map when no such block is used\. The global sum uses allNNnodes in the same sample, molecule, or graph\. The matricesWlocal,Wglobal∈ℝd×dW\_\{\\mathrm\{local\}\},W\_\{\\mathrm\{global\}\}\\in\\mathbb\{R\}^\{d\\times d\}mix feature channels, whileσ\\sigmaacts on the resulting node\-feature vector\. The damped solver then combines this map with the previous state according to Eq\. \([1](https://arxiv.org/html/2607.28989#S2.E1)\)\.

#### Unified local–global interaction sum\.

The split does not introduce two unrelated design pieces\. It can be written as one per\-node interaction sum whose sparse local channel and dense global channel have different support\. Define a one\-step propagatorΠi←j\\Pi\_\{i\\leftarrow j\}describing how source nodejjcontributes to receiver nodeii,

Πi←j​\(z~j\)=𝟏\{j∈𝒩​\(i\)\}​αi​j​Wlocal​z~j\+1N​Wglobal​z~j,\\Pi\_\{i\\leftarrow j\}\(\\tilde\{z\}\_\{j\}\)=\\mathbf\{1\}\_\{\\\{j\\in\\mathcal\{N\}\(i\)\\\}\}\\alpha\_\{ij\}\\,W\_\{\\mathrm\{local\}\}\\,\\tilde\{z\}\_\{j\}\+\\tfrac\{1\}\{N\}\\,W\_\{\\mathrm\{global\}\}\\,\\tilde\{z\}\_\{j\},\(27\)whereΠi←j\\Pi\_\{i\\leftarrow j\}maps the source feature vectorz~j\\tilde\{z\}\_\{j\}to a contribution at receiverii\. The indicator selects the local neighborhood term and the second term is the mean\-field contribution from every node\. The matricesWlocalW\_\{\\mathrm\{local\}\}andWglobalW\_\{\\mathrm\{global\}\}mix feature channels\. The update can then be written as a single sum over*all*NNnodes,fθ​\(z\(k\),x\)i=σ​\(ui\+∑j=1NΠi←j​\(z~j\(k\)\)\)f\_\{\\theta\}\(z^\{\(k\)\},x\)\_\{i\}=\\sigma\\big\(u\_\{i\}\+\\sum\_\{j=1\}^\{N\}\\Pi\_\{i\\leftarrow j\}\(\\tilde\{z\}\_\{j\}^\{\(k\)\}\)\\big\)\. The local/global split is therefore a support and operator choice within one interaction sum\. Collecting the same contributions at the level of the full state vector gives the linear operator’s closed form directly, using the adjacency matrixAA\(Ai​j=αi​jA\_\{ij\}=\\alpha\_\{ij\}forj∈𝒩​\(i\)j\\in\\mathcal\{N\}\(i\), else0\) and the Kronecker product⊗\\otimesto carry the per\-node channel mixingW∙W\_\{\\bullet\}across allNNnodes simultaneously\.

Jf​\(z∗\)≈A⊗Wlocal⏟Ωlocal\+1N​11⊤⊗Wglobal⏟Ωglobal,J\_\{f\}\(z^\{\\ast\}\)\\,\\approx\\,\\underbrace\{A\\otimes W\_\{\\mathrm\{local\}\}\}\_\{\\Omega\_\{\\mathrm\{local\}\}\}\\,\+\\,\\underbrace\{\\tfrac\{1\}\{N\}\\,\\mathbf\{1\}\\mathbf\{1\}^\{\\top\}\\otimes W\_\{\\mathrm\{global\}\}\}\_\{\\Omega\_\{\\mathrm\{global\}\}\},\(28\)Equation \([28](https://arxiv.org/html/2607.28989#S2.E28)\) is Eq\. \([19](https://arxiv.org/html/2607.28989#S2.E19)\) written out as an explicitN​d×N​dNd\\times Ndmatrix\. The first term is sparse with the local graph’s sparsity pattern\. The second term is rank\-1 dense in node space\. Since𝟏𝟏⊤\\mathbf\{1\}\\mathbf\{1\}^\{\\top\}is the all\-ones matrix, this term repeats the sameWglobalW\_\{\\mathrm\{global\}\}block at every\(i,j\)\(i,j\)pair, which is exactly what mean\-field means at the matrix level\. The expression is a linearization at fixed attention/gating weightsαi​j\\alpha\_\{ij\}, frozen at their currently evaluated value\. It does not makeLayerNorm​\(ReLU​\(⋅\)\)\\mathrm\{LayerNorm\}\(\\mathrm\{ReLU\}\(\\cdot\)\)in Eq\. \([3](https://arxiv.org/html/2607.28989#S2.E3)\) itself linear\.

### 2\.3Local interaction \(Ωlocal\\Omega\_\{\\mathrm\{local\}\}\)

For chemistry and graph\-node experiments, letyydenote the signal supplied to the local map\. Inside the solver,y=z~=χ𝒟​\(z\)y=\\tilde\{z\}=\\chi\_\{\\mathcal\{D\}\}\(z\)\. The local term follows GAT\-style attention\[[53](https://arxiv.org/html/2607.28989#bib.bib10)\],

L​\(y\)i\\displaystyle L\(y\)\_\{i\}=∑j∈𝒩​\(i\)ai​j​W​yj,\\displaystyle=\\sum\_\{j\\in\\mathcal\{N\}\(i\)\}a\_\{ij\}Wy\_\{j\},\(29\)ai​j\\displaystyle a\_\{ij\}=softmaxj​\(LeakyReLU​\(u⊤​\[W​yi∥W​yj\]\)\),\\displaystyle=\\mathrm\{softmax\}\_\{j\}\\\!\\left\(\\mathrm\{LeakyReLU\}\(u^\{\\top\}\[Wy\_\{i\}\\\|Wy\_\{j\}\]\)\\right\),\(30\)The neighborhood𝒩\\mathcal\{N\}is the molecular bond graph for chemistry or the citation graph for node classification, an instance of the general neural message\-passing framework\[[23](https://arxiv.org/html/2607.28989#bib.bib34)\], of which graph convolutional networks\[[34](https://arxiv.org/html/2607.28989#bib.bib18)\]are the special case with fixed, non\-attentive edge weights\. In this local\-attention formula,WWmaps node features into the attention feature space,uuis the learned attention vector,∥\\\|denotes feature concatenation, and the softmax is normalized over source nodesj∈𝒩​\(i\)j\\in\\mathcal\{N\}\(i\)for each receiverii\. Vision’s local term,LvisL\_\{\\mathrm\{vis\}\}\(Eq\. \([5](https://arxiv.org/html/2607.28989#S2.E5)\)\), is architecturally different\. It changes both the neighborhood construction and the update formula\. It is a degree\-normalized average through one shared linear map, with no per\-edge attention weightai​ja\_\{ij\}at all, over a symmetrizedkk\-nearest\-neighbor graph built in*channel space*\. The neighbors of channeliiare other hidden channels of the same sample, not other samples or graph nodes\.

### 2\.4Domain\-specific input encoders

Raw input reaches the equilibrium layer through a domain\-specific encoder,x=ℛϕ​\(xraw\)x=\\mathcal\{R\}\_\{\\phi\}\(x\_\{\\mathrm\{raw\}\}\)\. For MNIST the encoder is a single learned linear projection,

x\\displaystyle x=Wproj​xraw,\\displaystyle=W\_\{\\mathrm\{proj\}\}\\,x\_\{\\mathrm\{raw\}\},\(31\)xraw\\displaystyle x\_\{\\mathrm\{raw\}\}∈ℝ784,\\displaystyle\\in\\mathbb\{R\}^\{784\},\(32\)Wproj\\displaystyle W\_\{\\mathrm\{proj\}\}∈ℝd×784\.\\displaystyle\\in\\mathbb\{R\}^\{d\\times 784\}\.\(33\)Herexrawx\_\{\\mathrm\{raw\}\}is the flattened28×2828\\times 28image,WprojW\_\{\\mathrm\{proj\}\}is the learned encoder matrix, andx∈ℝdx\\in\\mathbb\{R\}^\{d\}is the hidden vector passed to the first equilibrium layer\. For CIFAR\-10 the encoder is a two\-block convolutional stack with batch normalization, max\-pooling, and dropout\[[28](https://arxiv.org/html/2607.28989#bib.bib48),[50](https://arxiv.org/html/2607.28989#bib.bib49)\]\. Fori=1,2i=1,2,

hi\\displaystyle h\_\{i\}=MaxPool2​\(ReLU​\(BNi​\(Convi​\(hi−1\)\)\)\),\\displaystyle=\\mathrm\{MaxPool\}\_\{2\}\\big\(\\mathrm\{ReLU\}\(\\mathrm\{BN\}\_\{i\}\(\\mathrm\{Conv\}\_\{i\}\(h\_\{i\-1\}\)\)\)\\big\),\(34\)h0\\displaystyle h\_\{0\}=xraw\.\\displaystyle=x\_\{\\mathrm\{raw\}\}\.\(35\)x=Wproj​Flatten​\(Dropoutp=0\.3​\(h2\)\),x=W\_\{\\mathrm\{proj\}\}\\,\\mathrm\{Flatten\}\\big\(\\mathrm\{Dropout\}\_\{p=0\.3\}\(h\_\{2\}\)\\big\),\(36\)Herehih\_\{i\}is the feature map after convolution blockii,Convi\\mathrm\{Conv\}\_\{i\}andBNi\\mathrm\{BN\}\_\{i\}are the convolution and batch normalization in that block, andMaxPool2\\mathrm\{MaxPool\}\_\{2\}denotes2×22\\times 2max\-pooling\. The matrixWprojW\_\{\\mathrm\{proj\}\}maps the flattened post\-dropout feature map to the hidden vectorx∈ℝdx\\in\\mathbb\{R\}^\{d\}\. The stack uses3→32→643\\\!\\to\\\!32\\\!\\to\\\!64output channels,3×33\\times 3kernels with padding 1, and dropout applied only during training\. The encoder is a shallow feature extractor that hands its representation to the equilibrium dynamics\.

### 2\.5Domain\-specific interaction graphs

The architecture’s motivating premise is that vision and chemistry can share the same damped local\+global equilibrium template even though neither the node set nor the local mechanism is identical between them\. Chemistry’s nodes are atoms with GAT\-style local attention over the fixed bond graph \(Section[2\.3](https://arxiv.org/html/2607.28989#S2.SS3)\), while vision’s “nodes,” for the purpose of its own local termLvisL\_\{\\mathrm\{vis\}\}, are theddhidden channels of one sample, with a degree\-normalized average \(not attention\) over

𝒩\(k\)​\(i\)=arg​top​\-​kj≠i⁡\(−‖z~i\(k\)−z~j\(k\)‖2\),\\mathcal\{N\}^\{\(k\)\}\(i\)=\\operatorname\*\{arg\\,top\\text\{\-\}k\}\_\{j\\neq i\}\\big\(\-\\\|\\tilde\{z\}\_\{i\}^\{\(k\)\}\-\\tilde\{z\}\_\{j\}^\{\(k\)\}\\\|\_\{2\}\\big\),\(37\)wherez~\(k\)=tanh⁡\(z\(k\)\)\\tilde\{z\}^\{\(k\)\}=\\tanh\(z^\{\(k\)\}\)for the vision models\. The neighborhood is recomputed from the current recurrent signal, not fromz∗z^\{\\ast\}, at every solver iteration\. Equation \([37](https://arxiv.org/html/2607.28989#S2.E37)\) is dynamic and content\-dependent, uses Euclidean distance over the current hidden state’s channels, and is symmetrized \(Section[2\.3](https://arxiv.org/html/2607.28989#S2.SS3)above\)\. The operatorarg​top​\-​k\\operatorname\*\{arg\\,top\\text\{\-\}k\}returns the set ofkksource\-channel indices with largest score, equivalently the nearest neighbors because the score is the negative distance\. The construction resembles dynamic graph CNNs, which recompute feature\-space neighborhoods inside the network, but the entities here are hidden channels rather than points in a point cloud\[[55](https://arxiv.org/html/2607.28989#bib.bib50)\]\. For chemistry the neighborhood graph𝒩​\(i\)\\mathcal\{N\}\(i\)is the molecule’s fixed bond graph, identical at every iteration\. The shared object is the damped local\+global equilibrium template, while the graph construction is domain\-specific\. Vision infers a dynamic hidden\-channel graph from recurrent state similarity, and chemistry uses the fixed molecular bond graph\.

### 2\.6Global interaction \(Ωglobal\\Omega\_\{\\mathrm\{global\}\}\)

For chemistry and graph\-node experiments, the global term computes a mean\-field vector within one system \(one molecule’sNNatoms or one graph’sNNnodes\) and broadcasts it back to that system’s own nodes\. For a single system with interaction signalyy,

G​\(y\)=1N​𝟏N​𝟏N⊤​y​Wg⊤,G\(y\)=\\tfrac\{1\}\{N\}\\mathbf\{1\}\_\{N\}\\mathbf\{1\}\_\{N\}^\{\\top\}yW\_\{g\}^\{\\top\},\(38\)wherey=z~y=\\tilde\{z\}inside the solver andWgW\_\{g\}is the learned channel\-mixing matrix\. The vector𝟏N\\mathbf\{1\}\_\{N\}is the all\-ones vector inℝN\\mathbb\{R\}^\{N\}, soN−1​𝟏N​𝟏N⊤​yN^\{\-1\}\\mathbf\{1\}\_\{N\}\\mathbf\{1\}\_\{N\}^\{\\top\}yreplaces each row by the graph or molecule mean state\. Equation \([38](https://arxiv.org/html/2607.28989#S2.E38)\) is a fixed averaging operator over nodes followed by learned channel mixing\. Vision’s global term,GvisG\_\{\\mathrm\{vis\}\}\(Eq\. \([5](https://arxiv.org/html/2607.28989#S2.E5)\)\), is not an instance of this mean\-field broadcast\. It is genuine self\-attention, but over one sample’s ownddhidden channels rather than pooled across samples or nodes, so there is noNN\-indexed pooling step inGvisG\_\{\\mathrm\{vis\}\}at all\. For node classification on real graphs, this reduces to anO​\(N\)O\(N\)scalar\-gated mean\-field broadcast\. The graph\-wide mean stateg=1N​∑jyjg=\\tfrac\{1\}\{N\}\\sum\_\{j\}y\_\{j\}is computed once, gated by a learned, data\-dependent scalar and redistributed to every node,

β\\displaystyle\\beta=σ​\(\(Wq​g\)⊤​\(Wk​g\)d\),\\displaystyle=\\sigma\\\!\\left\(\\frac\{\(W\_\{q\}g\)^\{\\top\}\(W\_\{k\}g\)\}\{\\sqrt\{d\}\}\\right\),\(39\)G​\(y\)i\\displaystyle G\(y\)\_\{i\}=β​Wg​g\.\\displaystyle=\\beta\\,W\_\{g\}\\,g\.\(40\)This broadcast is applied to every nodeii\. Hereg∈ℝdg\\in\\mathbb\{R\}^\{d\}is the mean state,β\\betais a scalar gate,ddis the hidden dimension, andWq,Wk,WgW\_\{q\},W\_\{k\},W\_\{g\}are learned matrices distinct from the local term’s own parameters\. The sigmoidσ\\sigmamaps the scalar attention score toβ∈\(0,1\)\\beta\\in\(0,1\)\. The gate uses the same bilinear,d\\sqrt\{d\}\-scaled construction as standard dot\-product attention \(Eq\. \([3](https://arxiv.org/html/2607.28989#S2.E3)\)’s own attention mechanism\), here applied once to the single graph\-wide vectorgginstead of pairwise between nodes\. This gives anO​\(N\)O\(N\)global term at graph scale\[[52](https://arxiv.org/html/2607.28989#bib.bib11)\]\. Full pairwise attention over a real citation graph produced an exactly uniform attention pattern in the citation\-network diagnostic below\. A static \(non\-gated\) variant fixesβ=1\\beta=1in Eq\. \([39](https://arxiv.org/html/2607.28989#S2.E39)\)\. A bounded top\-kkvariant replaces the single\-vector gate of Eq\. \([39](https://arxiv.org/html/2607.28989#S2.E39)\) with genuine node\-to\-node dot\-product attention, restricted to each node’skknearest candidates \(withO​\(N2\)O\(N^\{2\}\)score computation, avoiding the uniform\-attention collapse noted above since the softmax is always taken over onlykkcandidates regardless ofNN\),

G​\(y\)i\\displaystyle G\(y\)\_\{i\}=∑j∈𝒯k​\(i\)softmaxj​\(\(Wq​yi\)⊤​\(Wk​yj\)d\)​Wv​yj,\\displaystyle=\\sum\_\{j\\in\\mathcal\{T\}\_\{k\}\(i\)\}\\mathrm\{softmax\}\_\{j\}\\\!\\left\(\\frac\{\(W\_\{q\}y\_\{i\}\)^\{\\top\}\(W\_\{k\}y\_\{j\}\)\}\{\\sqrt\{d\}\}\\right\)W\_\{v\}y\_\{j\},\(41\)𝒯k​\(i\)\\displaystyle\\mathcal\{T\}\_\{k\}\(i\)=arg​top​\-​kj\(Wqyi\)⊤\(Wkyj\),\\displaystyle=\\operatorname\*\{arg\\,top\\text\{\-\}k\}\_\{j\}\(W\_\{q\}y\_\{i\}\)^\{\\top\}\(W\_\{k\}y\_\{j\}\),\(42\)recomputed from the current interaction signal every solver iteration, with𝒯k​\(i\)\\mathcal\{T\}\_\{k\}\(i\)the selected source\-node set for receiveriiandWq,Wk,WvW\_\{q\},W\_\{k\},W\_\{v\}a third, independent set of learned matrices for this variant\. The softmax in Eq\. \([41](https://arxiv.org/html/2607.28989#S2.E41)\) is normalized only over sourcesj∈𝒯k​\(i\)j\\in\\mathcal\{T\}\_\{k\}\(i\)\. The indexiiis the receiving node,jjis a candidate source node, andWvW\_\{v\}maps the selected source state into the value vector that is summed\. These variants define the global\-term alternatives used in the citation\-network ablations below\.

### 2\.7Spectral\-radius diagnostic at the fixed point

The damped\-step spectral radiusρ𝒯=ρ​\(J𝒯α​\(z∗\)\)\\rho\_\{\\mathcal\{T\}\}=\\rho\(J\_\{\\mathcal\{T\}\_\{\\alpha\}\}\(z^\{\\ast\}\)\)is estimated by power iteration on the vector–Jacobian product of the executed solver update \(Algorithm[2](https://arxiv.org/html/2607.28989#alg2)\)\. For a normal linearization, values below one indicate local one\-step contraction of the damped iteration\. More generally, because the Jacobian can be non\-normal,ρ𝒯\\rho\_\{\\mathcal\{T\}\}is used here as a local asymptotic stability indicator rather than as a singular\-value contraction certificate\. It is related to the Jacobian quantity whose crossing of 1 explains vanishing or exploding gradients in ordinary recurrent networks\[[42](https://arxiv.org/html/2607.28989#bib.bib31)\], here tracked on the equilibrium itself instead of on an unrolled training trajectory\. In DEQ models, growth of the equilibrium Jacobian is also known to affect both forward and backward solver stability\[[8](https://arxiv.org/html/2607.28989#bib.bib7)\]\. A Lyapunov exponent would require trajectory averaging oflog⁡‖J‖\\log\\\|J\\\|, withlog⁡ρ𝒯\\log\\rho\_\{\\mathcal\{T\}\}serving only as a local linearized proxy under the estimator used here\.

Algorithm 2Spectral radius estimation at the fixed point1:

v0∼𝒩​\(0,I\)v\_\{0\}\\sim\\mathcal\{N\}\(0,I\),

‖v0‖=1\\\|v\_\{0\}\\\|=1
2:for

t=1,…,Tt=1,\\dots,Tdo

3:

vt←J𝒯α​\(z∗\)⊤​vt−1v\_\{t\}\\leftarrow J\_\{\\mathcal\{T\}\_\{\\alpha\}\}\(z^\{\\ast\}\)^\{\\top\}v\_\{t\-1\}via the vector–Jacobian product \(VJP\) of the damped update

4:

ρ←‖vt‖\\rho\\leftarrow\\\|v\_\{t\}\\\|
5:

vt←vt/ρv\_\{t\}\\leftarrow v\_\{t\}/\\rho
6:endfor

In Algorithm[2](https://arxiv.org/html/2607.28989#alg2),vtv\_\{t\}is a unit vector in the flattened state space,TTis the number of power\-iteration steps, and VJP denotes a vector–Jacobian product evaluated at the fixed pointz∗z^\{\\ast\}\. The scalarρ\\rhoreturned after normalization is a power\-iteration estimate of the dominant eigenvalue magnitude ofJ𝒯α​\(z∗\)J\_\{\\mathcal\{T\}\_\{\\alpha\}\}\(z^\{\\ast\}\)under the usual spectral assumptions for that estimator\.

### 2\.8Energy diagnostic and local descent surrogate

A companion diagnostic toρ\\rhois a scalar quantity tracking alignment between the activated state and the global interaction acting on it, in the tradition of Hopfield\-style energy functions for recurrent networks\[[27](https://arxiv.org/html/2607.28989#bib.bib23)\]and the broader energy\-based view of learning\[[38](https://arxiv.org/html/2607.28989#bib.bib42)\]\. The MNIST solver\-dynamics analysis uses the evaluated diagnostic

ab​\(z\)\\displaystyle a\_\{b\}\(z\)=tanh⁡\(zb\),\\displaystyle=\\tanh\(z\_\{b\}\),\(43\)Ediag​\(z\)\\displaystyle E\_\{\\mathrm\{diag\}\}\(z\)=1B​∑b=1B\[ab​\(z\)⋅ab​\(z\)−ab​\(z\)⋅\(ab​\(z\)​Aglobal\(b\)\)\]\.\\displaystyle=\\frac\{1\}\{B\}\\sum\_\{b=1\}^\{B\}\\left\[a\_\{b\}\(z\)\\cdot a\_\{b\}\(z\)\-a\_\{b\}\(z\)\\cdot\\big\(a\_\{b\}\(z\)A\_\{\\mathrm\{global\}\}^\{\(b\)\}\\big\)\\right\]\.\(44\)Herebbindexes samples in a batch of sizeBB,ab​\(z\)a\_\{b\}\(z\)is the activated hidden state for samplebb, andAglobal\(b\)A\_\{\\mathrm\{global\}\}^\{\(b\)\}is the row\-normalized global\-attention matrix computed for that sample at the current solver iterate\. This is the quantity plotted in Figure[2](https://arxiv.org/html/2607.28989#S3.F2)\. It measures global\-attention self\-consistency in the vision model and does not include the local term\.

Treating a DEQ’s fixed point as the equilibrium of a nonlinear dynamical system and analyzing its stability directly is itself an active direction\. Chu et al\.\[[17](https://arxiv.org/html/2607.28989#bib.bib43)\]enforce Lyapunov stability of DEQ fixed points as an explicit training objective and show this improves adversarial robustness\. Here, the energy diagnostic andρ\\rhoare evaluated on networks trained without a stability constraint\. The robustness\-versus\-spectral\-radius analysis evaluates whether observed variation inρ\\rhopredicts robustness under input perturbations, which is a different question from the stability\-enforcing objective of\[[17](https://arxiv.org/html/2607.28989#bib.bib43)\]\.

#### Local descent surrogate for the damped iteration\.

The scalar diagnostic plotted in Figure[2](https://arxiv.org/html/2607.28989#S3.F2)and the quadratic surrogate below serve distinct roles\. For the linearized interpretation of the damped Picard step, hold the update vectorh=fθ​\(zk,x\)h=f\_\{\\theta\}\(z\_\{k\},x\)fixed at the current iterate and define

ℰh​\(z\)=12​‖z‖2−z⊤​h\.\\mathcal\{E\}\_\{h\}\(z\)=\\tfrac\{1\}\{2\}\\\|z\\\|^\{2\}\-z^\{\\top\}h\.\(45\)The state is interpreted as a flattened vector when the dot product is formed\. The surrogate is the local quadratic model obtained by freezing the nonlinear update, the same local\-linearization move used to defineΩ\\Omegain Section[2\.2](https://arxiv.org/html/2607.28989#S2.SS2)\. Its gradient is

∇zℰh​\(z\)=z−h\.\\nabla\_\{z\}\\mathcal\{E\}\_\{h\}\(z\)=z\-h\.\(46\)Setting∇zℰh​\(z\)=0\\nabla\_\{z\}\\mathcal\{E\}\_\{h\}\(z\)=0recovers exactly the fixed\-point conditionz∗=hz^\{\\ast\}=h, the frozen\-update version ofz∗=fθ​\(z∗,x\)z^\{\\ast\}=f\_\{\\theta\}\(z^\{\\ast\},x\)\. A single gradient\-descent step onℰh\\mathcal\{E\}\_\{h\}with step sizeη\\eta,

zk\+1\\displaystyle z\_\{k\+1\}=zk−η​∇zℰh​\(zk\),\\displaystyle=z\_\{k\}\-\\eta\\,\\nabla\_\{z\}\\mathcal\{E\}\_\{h\}\(z\_\{k\}\),\(47\)zk\+1\\displaystyle z\_\{k\+1\}=zk−η​\(zk−h\),\\displaystyle=z\_\{k\}\-\\eta\(z\_\{k\}\-h\),\(48\)zk\+1\\displaystyle z\_\{k\+1\}=\(1−η\)​zk\+η​h,\\displaystyle=\(1\-\\eta\)z\_\{k\}\+\\eta\\,h,\(49\)is, term for term, the damped Picard update of Eq\. \([1](https://arxiv.org/html/2607.28989#S2.E1)\) withη=α\\eta=\\alphaandhhidentified withfθ​\(zk,x\)f\_\{\\theta\}\(z\_\{k\},x\)\. Under this local surrogate,ℰh\\mathcal\{E\}\_\{h\}is non\-increasing under the usual gradient\-descent step\-size condition for a locally quadratic objective\. The global energy function of the trained nonlinear network is not assumed\.

### 2\.9Fixed\-point acceleration

Under appropriate contractive assumptions, plain Picard iteration \(Eq\. \([1](https://arxiv.org/html/2607.28989#S2.E1)\)\) converges linearly\. Several experiments instead use Anderson acceleration\[[3](https://arxiv.org/html/2607.28989#bib.bib20),[54](https://arxiv.org/html/2607.28989#bib.bib36)\], the same fixed\-point extrapolation technique used in the original DEQ solver\[[6](https://arxiv.org/html/2607.28989#bib.bib5)\], which extrapolates from themAm\_\{\\mathrm\{A\}\}most recent iterates by solving a constrained least\-squares problem for mixing coefficientsγ\\gammathat minimize the residualrk=fθ​\(zk,x\)−zkr\_\{k\}=f\_\{\\theta\}\(z\_\{k\},x\)\-z\_\{k\}in their span\.

γ\\displaystyle\\gamma=arg​minγ~⁡‖∑j=0mA−1γ~j​rk−j‖2,\\displaystyle=\\operatorname\*\{arg\\,min\}\_\{\\tilde\{\\gamma\}\}\\Big\\\|\\textstyle\\sum\_\{j=0\}^\{m\_\{\\mathrm\{A\}\}\-1\}\\tilde\{\\gamma\}\_\{j\}r\_\{k\-j\}\\Big\\\|\_\{2\},\(50\)∑jγ~j\\displaystyle\\sum\_\{j\}\\tilde\{\\gamma\}\_\{j\}=1,\\displaystyle=1,\(51\)zk\+1\\displaystyle z\_\{k\+1\}=∑j=0mA−1γj​fθ​\(zk−mA\+j\+1,x\)\.\\displaystyle=\\textstyle\\sum\_\{j=0\}^\{m\_\{\\mathrm\{A\}\}\-1\}\\gamma\_\{j\}\\,f\_\{\\theta\}\(z\_\{k\-m\_\{\\mathrm\{A\}\}\+j\+1\},x\)\.\(52\)HeremAm\_\{\\mathrm\{A\}\}is the Anderson memory length, distinct from the stimulus symbolmmused in the scalar examples\. The vectorrkr\_\{k\}is the residual of the fixed\-point equation at iteratekk,γ~\\tilde\{\\gamma\}is the optimization variable in the least\-squares problem, andγj\\gamma\_\{j\}are the resulting mixing weights over the lastmAm\_\{\\mathrm\{A\}\}residuals\. The constraint in Eq\. \([51](https://arxiv.org/html/2607.28989#S2.E51)\) makes the extrapolation affine, so the next iterate is a weighted combination of recent update\-map evaluations rather than a newly parameterized layer\.

### 2\.10Path\-sum \(Neumann series\) interpretation

The architecture’s global\-context design is motivated by a path\-sum argument\. IfΩ\\Omegaacted only once, indirect influence between nodes would be lost\. Summing every interaction path givesΩ\+Ω2\+Ω3\+⋯=\(I−Ω\)−1​Ω\\Omega\+\\Omega^\{2\}\+\\Omega^\{3\}\+\\cdots=\(I\-\\Omega\)^\{\-1\}\\Omega, converging wheneverρ​\(Ω\)<1\\rho\(\\Omega\)<1as in Eq\. \([13](https://arxiv.org/html/2607.28989#S2.E13)\)\. The derivation then connects this infinite\-series picture to the finite\-iteration solver, since convergence alone does not say how much of the series a damped budget ofKKiterations resolves\.

#### Three\-node scalar path expansion\.

Consider three nodes\{1,2,3\}\\\{1,2,3\\\}, each receiving a stimulusmim\_\{i\}and interacting through an elementary3×33\\times 3matrixAAwith entriesai​ja\_\{ij\}\(the influence of nodejjon nodeii\)\. The update ruleb=m\+A​tanh⁡\(b\)b=m\+A\\tanh\(b\)says nodeiiintegrates its own stimulusmim\_\{i\}plus a weighted sum of every other node’s nonlinear signal,∑jai​j​tanh⁡\(bj\)\\sum\_\{j\}a\_\{ij\}\\tanh\(b\_\{j\}\), including a self\-loop atj=ij=i\. Linearizing \(tanh⁡\(z\)≈z\\tanh\(z\)\\approx zfor small activations\) makes the path structure explicit\. Starting fromb\(0\)=mb^\{\(0\)\}=mand iteratingb\(t\+1\)≈m\+A​b\(t\)b^\{\(t\+1\)\}\\approx m\+Ab^\{\(t\)\},

b\(1\)\\displaystyle b^\{\(1\)\}≈m\+A​m,\\displaystyle\\approx m\+Am,\(53\)b\(2\)\\displaystyle b^\{\(2\)\}≈m\+A​m\+A2​m,\\displaystyle\\approx m\+Am\+A^\{2\}m,\(54\)b\(3\)\\displaystyle b^\{\(3\)\}≈m\+A​m\+A2​m\+A3​m\.\\displaystyle\\approx m\+Am\+A^\{2\}m\+A^\{3\}m\.\(55\)and by inductionb\(T\)≈∑r=0TAr​mb^\{\(T\)\}\\approx\\sum\_\{r=0\}^\{T\}A^\{r\}m\. Therr\-step term\(Ar​m\)i\(A^\{r\}m\)\_\{i\}sums every length\-rrpath into nodeii\. Collecting all path lengths gives the effective interactionAeff=A\+A2\+A3\+⋯=\(I−A\)−1​AA\_\{\\mathrm\{eff\}\}=A\+A^\{2\}\+A^\{3\}\+\\cdots=\(I\-A\)^\{\-1\}Awheneverρ​\(A\)<1\\rho\(A\)<1\. In the linear regime,b∗≈\(I−A\)−1​mb^\{\\ast\}\\approx\(I\-A\)^\{\-1\}mexactly, so the path interpretation is the closed form rather than an analogy\. The same statement holds for the locally linearized damped finite\-KKsolver, withAAgeneralized toM=ML\+MGM=M\_\{L\}\+M\_\{G\}\.

For the evaluated damped solver, the local linear form in Eq\. \([10](https://arxiv.org/html/2607.28989#S2.E10)\) can be written

zt\+1\\displaystyle z\_\{t\+1\}=α​s\+Tα​zt,\\displaystyle=\\alpha s\+T\_\{\\alpha\}z\_\{t\},\(56\)Tα\\displaystyle T\_\{\\alpha\}=\(1−α\)​I\+α​M,\\displaystyle=\(1\-\\alpha\)I\+\\alpha M,\(57\)M\\displaystyle M=ML\+MG\.\\displaystyle=M\_\{L\}\+M\_\{G\}\.\(58\)This giveszK=α​∑t=0K−1Tαt​sz\_\{K\}=\\alpha\\sum\_\{t=0\}^\{K\-1\}T\_\{\\alpha\}^\{t\}s\. ExpandingTαtT\_\{\\alpha\}^\{t\}gives all stimulus\-to\-state paths made of self\-persistence steps and local/global interaction steps\. Because the self\-persistence term is a scalar multiple of the identity, collecting those self steps yields the weighted powers ofMMbelow\.

M=Jf​\(z∗\)M=J\_\{f\}\(z^\{\\ast\}\)decomposes at the local\-linear level asM=ML\+MGM=M\_\{L\}\+M\_\{G\}\. The factors inMLM\_\{L\}andMGM\_\{G\}include the derivative ofz↦z~z\\mapsto\\tilde\{z\}and any active outer normalization or nonlinearity, matchingΩlocal\+Ωglobal\\Omega\_\{\\mathrm\{local\}\}\+\\Omega\_\{\\mathrm\{global\}\}of Eq\. \([19](https://arxiv.org/html/2607.28989#S2.E19)\)\. Solving the linearized systemz=s\+M​zz=s\+Mz\(s=Wstim​xs=W\_\{\\mathrm\{stim\}\}x\) by damped Picard iteration fromz0=0z\_\{0\}=0gives, in closed form,

zK\\displaystyle z\_\{K\}=∑r=0K−1wr\(K,α\)​Mr​s,\\displaystyle=\\sum\_\{r=0\}^\{K\-1\}w\_\{r\}^\{\(K,\\alpha\)\}M^\{r\}s,\(59\)wr\(K,α\)\\displaystyle w\_\{r\}^\{\(K,\\alpha\)\}=α​∑q=rK−1\(qr\)​\(1−α\)q−r​αr\.\\displaystyle=\\alpha\\sum\_\{q=r\}^\{K\-1\}\\binom\{q\}\{r\}\(1\-\\alpha\)^\{q\-r\}\\alpha^\{r\}\.\(60\)Hererris the interaction length,qqis the solver\-step index being summed out, andwr\(K,α\)w\_\{r\}^\{\(K,\\alpha\)\}is the finite\-budget weight assigned to all length\-rrinteraction products\. SinceM=ML\+MGM=M\_\{L\}\+M\_\{G\}and matrix products do not commute,Mr=∑w∈\{L,G\}rMwr​⋯​Mw1M^\{r\}=\\sum\_\{w\\in\\\{L,G\\\}^\{r\}\}M\_\{w\_\{r\}\}\\cdots M\_\{w\_\{1\}\}\. It sums every ordered length\-rrpath through the two\-generator interaction algebra, where path means an ordered product of local and global propagators acting on the stimulus\. Equation \([59](https://arxiv.org/html/2607.28989#S2.E59)\) is an exact finite path\-sum expansion with two propagator types and external sourcess\. AsK→∞K\\to\\infty,wr\(K,α\)→1w\_\{r\}^\{\(K,\\alpha\)\}\\to 1and Eq\. \([59](https://arxiv.org/html/2607.28989#S2.E59)\) recovers the Neumann seriesz∗=\(I−M\)−1​s=∑r=0∞Mr​sz^\{\\ast\}=\(I\-M\)^\{\-1\}s=\\sum\_\{r=0\}^\{\\infty\}M^\{r\}s, independent ofα\\alpha, which only governs which partial sum is realized at finiteKK\.

Asρ​\(M\)→1−\\rho\(M\)\\to 1^\{\-\}, termsMr​sM^\{r\}sdecay more slowly withrr\. Longer interaction histories matter, and finiteKKincreasingly under\-resolves them\. This gives the edge\-of\-chaos diagnostic \(Figure[7](https://arxiv.org/html/2607.28989#S3.F7)\) a precise meaning, in the sense already established for recurrent and deep networks at the boundary between ordered and chaotic dynamics\[[13](https://arxiv.org/html/2607.28989#bib.bib33),[43](https://arxiv.org/html/2607.28989#bib.bib32)\]\. It tracks how close the network sits to the boundary at which long\-range interaction paths stop being negligible\.

### 2\.11Linear response and effective influence

Equation \([59](https://arxiv.org/html/2607.28989#S2.E59)\) also defines the finite\-iteration response operator\. In the linear\-response sense, the Green’s operator maps an injected source or stimulus to the state produced by the locally linearized system\.

RK,α\\displaystyle R\_\{K,\\alpha\}=∑r=0K−1wr\(K,α\)​Mr,\\displaystyle=\\sum\_\{r=0\}^\{K\-1\}w\_\{r\}^\{\(K,\\alpha\)\}M^\{r\},\(61\)zK\\displaystyle z\_\{K\}=RK,α​s\.\\displaystyle=R\_\{K,\\alpha\}s\.\(62\)The finite solver budget realizesRK,αR\_\{K,\\alpha\}as its Green’s operator\. The indicesiiandjjnow index flattened state coordinates, which combine node and feature\-channel indices whenz∈ℝN×dz\\in\\mathbb\{R\}^\{N\\times d\}\. Its entry\(RK,α\)i​j\(R\_\{K,\\alpha\}\)\_\{ij\}is the total linearized effect on state coordinateiiof a unit stimulus at coordinatejj, after all interaction paths resolved byKKdamped iterations have been summed\. In the contractive infinite\-iteration limit,

R\\displaystyle R=limK→∞RK,α,\\displaystyle=\\lim\_\{K\\to\\infty\}R\_\{K,\\alpha\},\(63\)R\\displaystyle R=\(I−M\)−1,\\displaystyle=\(I\-M\)^\{\-1\},\(64\)independent ofα\\alpha\. This operator is an interpretability object attached to the locally linearized trained layer, derived from the update map rather than learned as an additional parameter\. It is connected to the fitted layer through the same update map and Jacobian\-vector machinery used for Algorithm[2](https://arxiv.org/html/2607.28989#alg2)\. The computation can use the linear action ofMMthrough vector–Jacobian or Jacobian–vector products without materializingMMorRR\.

### 2\.12Training objectives

Models are trained with negative log\-likelihood loss for classification,

ℒcls\\displaystyle\\mathcal\{L\}\_\{\\mathrm\{cls\}\}=−log⁡pθ​\(y∣z∗\),\\displaystyle=\-\\log p\_\{\\theta\}\(y\\mid z^\{\\ast\}\),\(65\)pθ\(⋅∣z∗\)\\displaystyle p\_\{\\theta\}\(\\cdot\\mid z^\{\\ast\}\)=softmax​\(Wout​z∗\),\\displaystyle=\\mathrm\{softmax\}\(W\_\{\\mathrm\{out\}\}z^\{\\ast\}\),\(66\)or mean absolute error for the ZINC regression target,

ℒreg\\displaystyle\\mathcal\{L\}\_\{\\mathrm\{reg\}\}=\|y^​\(z∗\)−y\|,\\displaystyle=\\big\|\\,\\hat\{y\}\(z^\{\\ast\}\)\-y\\,\\big\|,\(67\)y^​\(z∗\)\\displaystyle\\hat\{y\}\(z^\{\\ast\}\)=Wout​z∗,\\displaystyle=W\_\{\\mathrm\{out\}\}z^\{\\ast\},\(68\)Hereℒcls\\mathcal\{L\}\_\{\\mathrm\{cls\}\}is the per\-example classification loss,pθ​\(y∣z∗\)p\_\{\\theta\}\(y\\mid z^\{\\ast\}\)is the predicted probability assigned to the true classyy, andpθ\(⋅∣z∗\)p\_\{\\theta\}\(\\cdot\\mid z^\{\\ast\}\)is the full softmax distribution over classes\. The matrixWoutW\_\{\\mathrm\{out\}\}is the learned final readout\. It maps a node state, pooled graph state, or image hidden state to logits for classification, and maps the pooled molecular state to a scalar predictiony^\\hat\{y\}for regression\. In Eq\. \([67](https://arxiv.org/html/2607.28989#S2.E67)\),yyis the ground\-truth scalar property andℒreg\\mathcal\{L\}\_\{\\mathrm\{reg\}\}is the absolute error for one molecule\. For node classification, the readout is applied row\-wise to node states\. For graph\-level classification or regression,z∗z^\{\\ast\}is pooled before the final readout\. For image classification,z∗z^\{\\ast\}is the final hidden representation returned by the vision equilibrium stack\. Optimization uses Adam\[[33](https://arxiv.org/html/2607.28989#bib.bib40)\], gradient clipping\[[42](https://arxiv.org/html/2607.28989#bib.bib31)\], and step\-decay or plateau scheduling\. The losses contain no explicit physics term or fixed\-point residual term\. The dynamical quantities above are diagnostics evaluated on trained networks\.

### 2\.13Backward differentiation

SILVA models are differentiated in two ways in these experiments\. Most vision, chemistry, node\-classification, and graph\-benchmark models differentiate through the finiteKKsolver iterations\. The separate adjoint\-trained MNIST diagnostic uses a GMRES linear solve over the adjoint operator described in the Supplementary Methods\. The two backward modes share the same fixed\-point motivation, but they correspond to different numerical training procedures and are reported with different scope\.

#### Finite unrolled backward\.

Let

𝒯α,θ​\(z,x\)=\(1−α\)​z\+α​fθ​\(z,x\)\\mathcal\{T\}\_\{\\alpha,\\theta\}\(z,x\)=\(1\-\\alpha\)z\+\\alpha f\_\{\\theta\}\(z,x\)\(69\)be the damped solver map and letzk\+1=𝒯α,θ​\(zk,x\)z\_\{k\+1\}=\\mathcal\{T\}\_\{\\alpha,\\theta\}\(z\_\{k\},x\)fork=0,…,K−1k=0,\\ldots,K\-1\. The loss used in a finite\-solver experiment isℒK=ℒ​\(rψ​\(zK\),y\)\\mathcal\{L\}\_\{K\}=\\mathcal\{L\}\(r\_\{\\psi\}\(z\_\{K\}\),y\), whererψr\_\{\\psi\}is the task readout andyyis the target\. Reverse differentiation through the unrolled solver obeys

z¯K\\displaystyle\\bar\{z\}\_\{K\}=∇zKℒK,\\displaystyle=\\nabla\_\{z\_\{K\}\}\\mathcal\{L\}\_\{K\},\(70\)z¯k\\displaystyle\\bar\{z\}\_\{k\}=J𝒯α,θ​\(zk,x\)⊤​z¯k\+1,k=K−1,…,0\.\\displaystyle=J\_\{\\mathcal\{T\}\_\{\\alpha,\\theta\}\}\(z\_\{k\},x\)^\{\\top\}\\bar\{z\}\_\{k\+1\},\\qquad k=K\-1,\\ldots,0\.\(71\)Herez¯k\\bar\{z\}\_\{k\}is the adjoint gradient with respect to thekkth solver state, andJ𝒯α,θ​\(zk,x\)J\_\{\\mathcal\{T\}\_\{\\alpha,\\theta\}\}\(z\_\{k\},x\)is the Jacobian of the damped update map with respect tozzat that iterate\. The parameter and input gradients accumulate over all solver steps,

∇θℒK\\displaystyle\\nabla\_\{\\theta\}\\mathcal\{L\}\_\{K\}=∑k=0K−1\(∂θ𝒯α,θ​\(zk,x\)\)⊤​z¯k\+1,\\displaystyle=\\sum\_\{k=0\}^\{K\-1\}\\left\(\\partial\_\{\\theta\}\\mathcal\{T\}\_\{\\alpha,\\theta\}\(z\_\{k\},x\)\\right\)^\{\\top\}\\bar\{z\}\_\{k\+1\},\(72\)∇xℒK\\displaystyle\\nabla\_\{x\}\\mathcal\{L\}\_\{K\}=∑k=0K−1\(∂x𝒯α,θ​\(zk,x\)\)⊤​z¯k\+1\.\\displaystyle=\\sum\_\{k=0\}^\{K\-1\}\\left\(\\partial\_\{x\}\\mathcal\{T\}\_\{\\alpha,\\theta\}\(z\_\{k\},x\)\\right\)^\{\\top\}\\bar\{z\}\_\{k\+1\}\.\(73\)In the evaluated SILVA models, parameters enter the solver through stimulus, local\-interaction, global\-interaction, and normalization maps, while the task loss also differentiates the final readoutrψr\_\{\\psi\}\. Since∂θ𝒯α,θ=α​∂θfθ\\partial\_\{\\theta\}\\mathcal\{T\}\_\{\\alpha,\\theta\}=\\alpha\\,\\partial\_\{\\theta\}f\_\{\\theta\}for parameters entering the update map, Eqs\. \([71](https://arxiv.org/html/2607.28989#S2.E71)\) –\([73](https://arxiv.org/html/2607.28989#S2.E73)\) are the exact chain rule for the finite Picard loop used in those experiments\.

#### Implicit adjoint backward\.

For the damped equilibrium condition𝒯α​\(z,x\)=\(1−α\)​z\+α​fθ​\(z,x\)=z\\mathcal\{T\}\_\{\\alpha\}\(z,x\)=\(1\-\\alpha\)z\+\\alpha f\_\{\\theta\}\(z,x\)=z, the implicit function theorem givesJ𝒯=\(1−α\)​I\+α​JfθJ\_\{\\mathcal\{T\}\}=\(1\-\\alpha\)I\+\\alpha J\_\{f\_\{\\theta\}\}\. Substituting intoI−J𝒯I\-J\_\{\\mathcal\{T\}\}term by term,

I−J𝒯\\displaystyle I\-J\_\{\\mathcal\{T\}\}=I−\[\(1−α\)​I\+α​Jfθ\],\\displaystyle=I\-\\big\[\(1\-\\alpha\)I\+\\alpha J\_\{f\_\{\\theta\}\}\\big\],\(74\)I−J𝒯\\displaystyle I\-J\_\{\\mathcal\{T\}\}=α​\(I−Jfθ\),\\displaystyle=\\alpha\\big\(I\-J\_\{f\_\{\\theta\}\}\\big\),\(75\)HereJ𝒯J\_\{\\mathcal\{T\}\}is the Jacobian of the damped solver map𝒯α\\mathcal\{T\}\_\{\\alpha\}with respect tozz,JfθJ\_\{f\_\{\\theta\}\}is the corresponding Jacobian offθf\_\{\\theta\}, andIIis the identity on the flattened state space\. All Jacobians in this subsection are evaluated at the fixed point unless otherwise stated\. Consequently,\(I−J𝒯\)−1=1α​\(I−Jfθ\)−1\(I\-J\_\{\\mathcal\{T\}\}\)^\{\-1\}=\\tfrac\{1\}\{\\alpha\}\(I\-J\_\{f\_\{\\theta\}\}\)^\{\-1\}wheneverI−Jfθ​\(z∗\)I\-J\_\{f\_\{\\theta\}\}\(z^\{\\ast\}\)is invertible, which Eq\. \([13](https://arxiv.org/html/2607.28989#S2.E13)\) guarantees\. Applying the implicit function theorem to𝒯α​\(z∗,x\)=z∗\\mathcal\{T\}\_\{\\alpha\}\(z^\{\\ast\},x\)=z^\{\\ast\}with respect to any input or parameter quantityξ\\xientering throughfθf\_\{\\theta\}, differentiating both sides gives\(I−J𝒯\)​∂z∗/∂ξ=α​∂fθ/∂ξ\(I\-J\_\{\\mathcal\{T\}\}\)\\,\\partial z^\{\\ast\}/\\partial\\xi=\\alpha\\,\\partial f\_\{\\theta\}/\\partial\\xi, hence

∂z∗∂ξ=\(I−Jfθ\)−1​∂fθ∂ξ,\\frac\{\\partial z^\{\\ast\}\}\{\\partial\\xi\}=\\big\(I\-J\_\{f\_\{\\theta\}\}\\big\)^\{\-1\}\\frac\{\\partial f\_\{\\theta\}\}\{\\partial\\xi\},\(76\)Hereξ\\xidenotes any input, weight, or other differentiable quantity entering the update map throughfθf\_\{\\theta\}\. The derivatives are understood in vectorized state coordinates, so∂z∗/∂ξ\\partial z^\{\\ast\}/\\partial\\xiis the sensitivity of the fixed point to that quantity\. The damping coefficientα\\alphacancels exactly\. The fixed pointz∗z^\{\\ast\}itself does not depend onα\\alpha\(the fixed\-point condition holds for anyα≠0\\alpha\\neq 0\), so its derivative cannot depend onα\\alphaeither\. Damping changes only the solver trajectory used to reachz∗z^\{\\ast\}, not the destination\. This specializes the general equilibrium\-gradient result of\[[6](https://arxiv.org/html/2607.28989#bib.bib5)\], stated there for the undamped root conditiongθ​\(z\)=fθ​\(z,x\)−z=0g\_\{\\theta\}\(z\)=f\_\{\\theta\}\(z,x\)\-z=0directly\. The derivation above shows the same conclusion survives introducing the damped iteration of Eq\. \([1](https://arxiv.org/html/2607.28989#S2.E1)\) as an explicit solver choice\. The damping coefficient is absent from the resulting gradient\. The backward pass then solves theα\\alpha\-free adjoint system

\(I−Jfθ​\(z∗\)⊤\)​Δ=∇z∗ℒ\\left\(I\-J\_\{f\_\{\\theta\}\}\(z^\{\\ast\}\)^\{\\top\}\\right\)\\Delta=\\nabla\_\{z^\{\\ast\}\}\\mathcal\{L\}\(77\)HereΔ\\Deltais the adjoint vector and∇z∗ℒ\\nabla\_\{z^\{\\ast\}\}\\mathcal\{L\}is the gradient of the loss with respect to the fixed\-point state\. The corresponding fixed\-point gradients are

∇θℒ\\displaystyle\\nabla\_\{\\theta\}\\mathcal\{L\}=\(∂θfθ​\(z∗,x\)\)⊤​Δ,\\displaystyle=\\left\(\\partial\_\{\\theta\}f\_\{\\theta\}\(z^\{\\ast\},x\)\\right\)^\{\\top\}\\Delta,\(78\)∇xℒ\\displaystyle\\nabla\_\{x\}\\mathcal\{L\}=\(∂xfθ​\(z∗,x\)\)⊤​Δ\.\\displaystyle=\\left\(\\partial\_\{x\}f\_\{\\theta\}\(z^\{\\ast\},x\)\\right\)^\{\\top\}\\Delta\.\(79\)Here∂θfθ\\partial\_\{\\theta\}f\_\{\\theta\}collects the derivatives of the stimulus, local\-interaction, global\-interaction, and normalization maps with respect to their parameters\. Equivalently, one may solve\(I−J𝒯α​\(z∗\)⊤\)​μ=∇z∗ℒ\(I\-J\_\{\\mathcal\{T\}\_\{\\alpha\}\}\(z^\{\\ast\}\)^\{\\top\}\)\\mu=\\nabla\_\{z^\{\\ast\}\}\\mathcal\{L\}and apply\(∂θ𝒯α,θ\)⊤​μ\(\\partial\_\{\\theta\}\\mathcal\{T\}\_\{\\alpha,\\theta\}\)^\{\\top\}\\mu\. BecauseI−J𝒯α=α​\(I−Jfθ\)I\-J\_\{\\mathcal\{T\}\_\{\\alpha\}\}=\\alpha\(I\-J\_\{f\_\{\\theta\}\}\)and∂θ𝒯α,θ=α​∂θfθ\\partial\_\{\\theta\}\\mathcal\{T\}\_\{\\alpha,\\theta\}=\\alpha\\partial\_\{\\theta\}f\_\{\\theta\}, the damping factors cancel in the exact fixed\-point gradient\. This is the memory\-efficient equilibrium\-gradient identity used by DEQs\[[6](https://arxiv.org/html/2607.28989#bib.bib5)\]\. It requires the fixed point and Jacobian\-vector products, rather than the stored forward trajectory\. The exact identity above is the reference DEQ gradient\. The adjoint\-trained MNIST diagnostic uses the GMRES adjoint described in the Supplementary Methods, with hatted gradients used there to distinguish the executed adjoint path from the formal identity in Eqs\. \([77](https://arxiv.org/html/2607.28989#S2.E77)\)–\([79](https://arxiv.org/html/2607.28989#S2.E79)\)\. The remaining experiments use truncated backpropagation through the unrolledKK\-step Picard loop, a weight\-tied\-network training regime that does not carry theO​\(1\)O\(1\)\-memory guarantee\.

#### Static global\-matrix adjoint\.

For a single learned, non\-gated global matrixAA\(Section[2\.6](https://arxiv.org/html/2607.28989#S2.SS6)\), Eq\. \([77](https://arxiv.org/html/2607.28989#S2.E77)\) can be solved for∇Aℒ\\nabla\_\{A\}\\mathcal\{L\}in closed form\. Using the receiver\-first column convention of Eq\. \([27](https://arxiv.org/html/2607.28989#S2.E27)\), write the fixed point as a root condition onb∗b^\{\\ast\}\(pre\-activation state,z∗=tanh⁡\(b∗\)z^\{\\ast\}=\\tanh\(b^\{\\ast\}\)\),

g​\(b∗,A\)\\displaystyle g\(b^\{\\ast\},A\)=b∗−m−A​tanh⁡\(b∗\)=0,\\displaystyle=b^\{\\ast\}\-m\-A\\tanh\(b^\{\\ast\}\)=0,\(80\)m\\displaystyle m=Wstim​x\.\\displaystyle=W\_\{\\mathrm\{stim\}\}x\.\(81\)Hereb∗,m∈ℝNb^\{\\ast\},m\\in\\mathbb\{R\}^\{N\},A∈ℝN×NA\\in\\mathbb\{R\}^\{N\\times N\}, andtanh\\tanhis applied elementwise\. The root functiong​\(b∗,A\)g\(b^\{\\ast\},A\)is zero exactly when the pre\-activation equilibrium condition is satisfied\. Contraction follows fromtanh\\tanh’s 1\-Lipschitz property,\|tanh⁡\(u\)−tanh⁡\(v\)\|≤\|u−v\|\|\\tanh\(u\)\-\\tanh\(v\)\|\\leq\|u\-v\|\. For anyb1,b2b\_\{1\},b\_\{2\},

‖fθ​\(b1,x\)−fθ​\(b2,x\)‖\\displaystyle\\\|f\_\{\\theta\}\(b\_\{1\},x\)\-f\_\{\\theta\}\(b\_\{2\},x\)\\\|=‖A​\(tanh⁡\(b1\)−tanh⁡\(b2\)\)‖,\\displaystyle=\\\|A\(\\tanh\(b\_\{1\}\)\-\\tanh\(b\_\{2\}\)\)\\\|,\(82\)‖fθ​\(b1,x\)−fθ​\(b2,x\)‖\\displaystyle\\\|f\_\{\\theta\}\(b\_\{1\},x\)\-f\_\{\\theta\}\(b\_\{2\},x\)\\\|≤‖A‖2​‖b1−b2‖,\\displaystyle\\leq\\\|A\\\|\_\{2\}\\,\\\|b\_\{1\}\-b\_\{2\}\\\|,\(83\)Hereb1b\_\{1\}andb2b\_\{2\}are two candidate pre\-activation states,∥⋅∥\\\|\\cdot\\\|is the Euclidean vector norm, and‖A‖2\\\|A\\\|\_\{2\}is the spectral norm of the interaction matrix\. Thus‖A‖2<1\\\|A\\\|\_\{2\}<1is a checkable\-at\-initialization sufficient condition for this case to be a contraction, stronger than Eq\. \([13](https://arxiv.org/html/2607.28989#S2.E13)\)’s exact criterion \(the operator norm dominates the spectral radius for non\-normalAA, so this bound can be conservative\. The same pattern appears in Figure[9](https://arxiv.org/html/2607.28989#S3.F9)showing several configurations withρ\\rhocomfortably below 1 at initialization despite no explicit norm constraint onAA\)\. Totally differentiatingg=0g=0at fixedmm, withJ=A​diag​\(1−tanh2⁡\(b∗\)\)J=A\\,\\mathrm\{diag\}\(1\-\\tanh^\{2\}\(b^\{\\ast\}\)\),

d​g\\displaystyle dg=\(I−J\)​d​b∗−d​A​tanh⁡\(b∗\)=0,\\displaystyle=\(I\-J\)\\,db^\{\\ast\}\-dA\\,\\tanh\(b^\{\\ast\}\)=0,\(84\)d​b∗\\displaystyle db^\{\\ast\}=\(I−J\)−1​d​A​tanh⁡\(b∗\)\.\\displaystyle=\(I\-J\)^\{\-1\}dA\\,\\tanh\(b^\{\\ast\}\)\.\(85\)Hered​gdg,d​b∗db^\{\\ast\}, andd​AdAare differentials of the root function, fixed point, and interaction matrix\. The matrixJ=A​diag​\(1−tanh2⁡\(b∗\)\)J=A\\,\\mathrm\{diag\}\(1\-\\tanh^\{2\}\(b^\{\\ast\}\)\)is the Jacobian ofA​tanh⁡\(b∗\)A\\tanh\(b^\{\\ast\}\)with respect tob∗b^\{\\ast\}\. Given the loss gradientv=∇b∗ℒv=\\nabla\_\{b^\{\\ast\}\}\\mathcal\{L\}, the adjoint vectoruusolves\(I−J\)⊤​u=v\(I\-J\)^\{\\top\}u=v, the same linear solve as Eq\. \([77](https://arxiv.org/html/2607.28989#S2.E77)\) specialized to this static matrix case\. The parameter gradient follows by substituting Eq\. \([84](https://arxiv.org/html/2607.28989#S2.E84)\) intod​ℒ=v⊤​d​b∗d\\mathcal\{L\}=v^\{\\top\}db^\{\\ast\}\.

∇Aℒ=utanh\(b∗\)⊤\.\\nabla\_\{A\}\\mathcal\{L\}=u\\,\\tanh\(b^\{\\ast\}\)^\{\\top\}\.\(86\)The right\-hand side is an outer product, so the resulting gradient has the sameN×NN\\times Nshape asAA\. Its\(i,j\)\(i,j\)entry is the sensitivity of the loss to the influence from source nodejjto receiver nodeii\. The equivalent row\-oriented tensor convention writes the same interaction asϕ​\(b∗\)​Arow\\phi\(b^\{\\ast\}\)A\_\{\\mathrm\{row\}\}, which stores the transpose of the receiver\-first convention above\. In that notation, the per\-example gradient isϕ​\(b∗\)⊤​u\\phi\(b^\{\\ast\}\)^\{\\top\}ubefore summing over a minibatch or differentiating through an attention\-generatedArowA\_\{\\mathrm\{row\}\}\.

When the global matrix is itself generated from the stimulus, the outer\-product adjoint first gives the sensitivity with respect toAglobalA\_\{\\mathrm\{global\}\}\. Differentiating the attention mapAglobal​\(m\)A\_\{\\mathrm\{global\}\}\(m\)then propagates this sensitivity to the stimulus and to the parameters that define the global interaction\. The resulting adjoint is the sum of the direct stimulus contribution and the contribution mediated by the global\-attention map\.

#### Head\-averaged attention adjoint\.

This case applies to local\-attention layers in which theHHhead outputs are averaged to a single per\-node value, including the chemistry configuration with averaged rather than concatenated heads and the adjoint\-differentiated layer\. Let𝒜φ​\(y\)=\{ηi,h\}i,h\\mathcal\{A\}\_\{\\varphi\}\(y\)=\\\{\\eta\_\{i,h\}\\\}\_\{i,h\}denote the multi\-head attention map before head aggregation, whereφ\\varphiare the attention parameters andηi,h\\eta\_\{i,h\}is the output of headhhat nodeii\. The averaged local output is

ℓi\\displaystyle\\ell\_\{i\}=1H​∑h=1Hηi,h,\\displaystyle=\\frac\{1\}\{H\}\\sum\_\{h=1\}^\{H\}\\eta\_\{i,h\},\(87\)η¯i,h\\displaystyle\\bar\{\\eta\}\_\{i,h\}=1H​ℓ¯i,i=1,…,N,h=1,…,H,\\displaystyle=\\frac\{1\}\{H\}\\bar\{\\ell\}\_\{i\},\\qquad i=1,\\ldots,N,\\quad h=1,\\ldots,H,\(88\)y¯local\\displaystyle\\bar\{y\}\_\{\\mathrm\{local\}\}=\(Dy​𝒜φ​\(y\)\)⊤​η¯,\\displaystyle=\\left\(D\_\{y\}\\mathcal\{A\}\_\{\\varphi\}\(y\)\\right\)^\{\\top\}\\bar\{\\eta\},\(89\)∇φℒ\\displaystyle\\nabla\_\{\\varphi\}\\mathcal\{L\}=\(Dφ​𝒜φ​\(y\)\)⊤​η¯\.\\displaystyle=\\left\(D\_\{\\varphi\}\\mathcal\{A\}\_\{\\varphi\}\(y\)\\right\)^\{\\top\}\\bar\{\\eta\}\.\(90\)Hereℓ¯i=∂ℒ/∂ℓi\\bar\{\\ell\}\_\{i\}=\\partial\\mathcal\{L\}/\\partial\\ell\_\{i\}is the incoming adjoint for the averaged local output,η¯i,h=∂ℒ/∂ηi,h\\bar\{\\eta\}\_\{i,h\}=\\partial\\mathcal\{L\}/\\partial\\eta\_\{i,h\}is the adjoint for an individual head output, andy¯local\\bar\{y\}\_\{\\mathrm\{local\}\}is the local\-attention contribution to the adjoint of the signalyy\. The matricesDy​𝒜φD\_\{y\}\\mathcal\{A\}\_\{\\varphi\}andDφ​𝒜φD\_\{\\varphi\}\\mathcal\{A\}\_\{\\varphi\}are the Jacobians of the multi\-head attention map with respect to its input signal and parameters\. Thus head averaging contributes only the factor1/H1/Hin Eq\. \([88](https://arxiv.org/html/2607.28989#S2.E88)\); the remaining local adjoint is the standard vector–Jacobian product of the attention map\. The local and global contributions then combine as in Eq\. \([77](https://arxiv.org/html/2607.28989#S2.E77)\)\. The global term is unaffected by the head\-averaging reduction\. Citation\-network and CLUSTER layers use concatenated attention heads in the finite unrolled training path, so their gradients are governed by Eqs\. \([71](https://arxiv.org/html/2607.28989#S2.E71)\)–\([72](https://arxiv.org/html/2607.28989#S2.E72)\)\.

### 2\.14Spectral normalization

Spectral normalizationWSN=W/σmax​\(W\)W\_\{\\mathrm\{SN\}\}=W/\\sigma\_\{\\max\}\(W\)\[[41](https://arxiv.org/html/2607.28989#bib.bib51)\]on a weight matrix appearing linearly infθf\_\{\\theta\}makes that linear map non\-expansive \(‖WSN​z‖≤‖z‖\\\|W\_\{\\mathrm\{SN\}\}z\\\|\\leq\\\|z\\\|, by submultiplicativity of the operator norm\), so that each component force is individually bounded\. Non\-expansiveness of individual linear maps does not by itself make their sum strictly contractive\. This mirrors the weight\-norm\-constrained view of stability used to design provably stable deep and ordinary\-differential\-equation\-inspired architectures\[[25](https://arxiv.org/html/2607.28989#bib.bib19)\], and the same argument as Eq\. \([83](https://arxiv.org/html/2607.28989#S2.E83)\) applied per weight matrix instead of to a single static matrix\. Hereσmax​\(W\)\\sigma\_\{\\max\}\(W\)is the largest singular value ofWW, andWSNW\_\{\\mathrm\{SN\}\}is the normalized version used in the forward map\. In these experiments, spectral normalization is applied to the molecular\-regression model’s global and stimulus weight matrices\. Only the state\-dependent parts of the update contribute to the Jacobian with respect tozz, so the resulting stability is monitored through the damped\-step diagnosticρ𝒯\\rho\_\{\\mathcal\{T\}\}\(Section[2\.7](https://arxiv.org/html/2607.28989#S2.SS7)\) rather than assumed from weight normalization alone\.

## 3Results

#### MNIST solver dynamics\.

Two MNIST training regimes are used for different purposes\. The principal finite\-unrolled accuracy experiment uses a 90/10 train/validation split with a held\-out test set never used for model selection\. A five\-seed sweep of the single\-layer MNIST local\+global architecture\[[37](https://arxiv.org/html/2607.28989#bib.bib52)\]\(α=0\.25\\alpha=0\.25,K=20K=20\) reaches97\.23±0\.13%97\.23\\pm 0\.13\\%test accuracy\. Figure[2](https://arxiv.org/html/2607.28989#S3.F2)compares Picard and Anderson solvers on the same five trained models\. In Figure[2](https://arxiv.org/html/2607.28989#S3.F2)a, the Anderson trajectory has already reduced the residual by its first plotted accelerated step and stays below Picard through the early function evaluations\. By 20–30 evaluations, both methods reach the same low\-residual regime, with Picard marginally lower at the end of this sweep\. This supports the interpretation that Anderson changes the route to the fixed point rather than the learned equilibrium being evaluated\. Figure[2](https://arxiv.org/html/2607.28989#S3.F2)b gives the paired solver\-state diagnostic\. Picard rises from the initialized state, whereas Anderson relaxes from a higher initial diagnostic value, and both trajectories approach the same plateau\. The diagnostic therefore tracks convergence toward a common trained state\. It is not used as a proof of global Lyapunov descent\.

![Refer to caption](https://arxiv.org/html/2607.28989v1/x2.png)Figure 2:Convergence diagnostics for trained SILVA image classifiers\. \(a\) Damped fixed\-point residual‖𝒯θ​\(zk;x\)−zk‖2\\\|\\mathcal\{T\}\_\{\\theta\}\(z\_\{k\};x\)\-z\_\{k\}\\\|\_\{2\}over solver function evaluations\. \(b\) Energy diagnostic over the same iterations\. Lines show means and bands show standard deviations across five independently trained MNIST models\. Axes report function\-evaluation counts and arbitrary units for the residual and energy\-diagnostic magnitudes\.
#### Implicit\-adjoint MNIST equilibrium\.

A separate MNIST experiment trained with a GMRES adjoint, rather than truncated backpropagation through the unrolled solver, learned a non\-trivial classifier but remained below the finite\-unrolled MNIST configuration\. It is used here as an equilibrium\-state separability diagnostic\. Figure[3](https://arxiv.org/html/2607.28989#S3.F3)shows principal component analysis \(PCA\) and t\-SNE projections\[[31](https://arxiv.org/html/2607.28989#bib.bib53),[51](https://arxiv.org/html/2607.28989#bib.bib54)\]of the converged hidden state across the MNIST test set at 89\.9% validation accuracy\. For this diagnostic, the fixed\-point solver is run on each test image and the terminal statez∗z^\{\\ast\}is stored before the classification readout\. PCA displays the two largest variance directions of these full hidden states, while t\-SNE displays their local neighborhood organization in a nonlinear embedding\. The plotted points are therefore projections of vector attractors rather than logits or input pixels\. In Figure[3](https://arxiv.org/html/2607.28989#S3.F3)a, the first two linear principal components separate several digit regions but also retain substantial overlap\. The equilibrium state therefore contains class structure without collapsing into ten linearly disjoint regions in this projection\. In Figure[3](https://arxiv.org/html/2607.28989#S3.F3)b, the nonlinear neighborhood embedding exposes more coherent class neighborhoods, with residual mixing between visually related or ambiguous classes\. The silhouette score0\.3820\.382and Davies–Bouldin score1\.0981\.098\[[45](https://arxiv.org/html/2607.28989#bib.bib55),[19](https://arxiv.org/html/2607.28989#bib.bib56)\]are computed on the fullz∗z^\{\\ast\}vectors before two\-dimensional projection\. They quantify this moderate, rather than complete, attractor separation\.

![Refer to caption](https://arxiv.org/html/2607.28989v1/figures/fig_attractor_landscape_real.png)Figure 3:Input\-conditioned equilibrium\-state geometry from the converged hidden statez∗z^\{\\ast\}on MNIST test data from the adjoint\-trained diagnostic model\. \(a\) PCA projection\. \(b\) t\-SNE projection\. Each point is one test image represented by its terminal hidden state before readout\. Points are colored by digit class, and projection coordinates are reported in arbitrary units\.
#### Stacked equilibrium dynamics on CIFAR\-10\.

The same damped fixed\-point template, with the CIFAR\-specific convolutional input encoder \(Section[2](https://arxiv.org/html/2607.28989#S2)\), reaches74\.77±1\.59%74\.77\\pm 1\.59\\%test accuracy on CIFAR\-10\[[36](https://arxiv.org/html/2607.28989#bib.bib57)\]over five independent seeds at hidden width 64\. Figure[4](https://arxiv.org/html/2607.28989#S3.F4)traces the Picard\-iteration trajectory, not only the converged endpoint, for both equilibrium layers across all five trained seeds, one held\-out test image per class\. Each curve follows the hidden state sequencez0,z1,…,zKz\_\{0\},z\_\{1\},\\ldots,z\_\{K\}through solver time after projection into the panel’s PCA plane\. Color indicates CIFAR\-10 class\. The large terminal markers highlight representative approximate equilibria, and the length of a curve reports change in the projected hidden state rather than distance in the image space\. In Figure[4](https://arxiv.org/html/2607.28989#S3.F4)a, the first layer withα=0\.5\\alpha=0\.5moves through a broad state\-space range before settling, and the endpoints remain class\-dependent across seeds\. In Figure[4](https://arxiv.org/html/2607.28989#S3.F4)b, the second layer withα=0\.2\\alpha=0\.2follows shorter trajectories in a more compact projected region while preserving class\-conditioned endpoints\. The stacked model therefore does not merely add depth\. It applies a second equilibrium transform with a different relaxation scale to the representation produced by the first\.

![Refer to caption](https://arxiv.org/html/2607.28989v1/figures/fig_hierarchy_real.png)Figure 4:Stacked equilibrium hierarchy on CIFAR\-10\. PCA\-projected Picard\-iteration trajectories are shown for five independently\-trained seeds and one held\-out image per class\. Curves show solver\-time hidden\-state trajectories, and large markers highlight representative terminal equilibria\. \(a\) Fast first layer withα=0\.5\\alpha=0\.5\. \(b\) Slow second layer withα=0\.2\\alpha=0\.2\. PCA coordinates are reported in arbitrary units\.
#### Molecular property regression on ZINC\.

Property\-regression targets are drawn from the ZINC molecular\-property benchmark\[[29](https://arxiv.org/html/2607.28989#bib.bib21)\]\(12k subset\)\. Table[1](https://arxiv.org/html/2607.28989#S3.T1)gives mean absolute error \(MAE\) for graph convolutional network \(GCN\), graph attention network \(GAT\), graph isomorphism network \(GIN\), principal neighbourhood aggregation \(PNA\), and gated graph convolutional network \(GatedGCN\) baselines\[[34](https://arxiv.org/html/2607.28989#bib.bib18),[53](https://arxiv.org/html/2607.28989#bib.bib10),[57](https://arxiv.org/html/2607.28989#bib.bib12),[18](https://arxiv.org/html/2607.28989#bib.bib16),[14](https://arxiv.org/html/2607.28989#bib.bib58),[20](https://arxiv.org/html/2607.28989#bib.bib17)\], each trained under the same ZINC split, hidden width, and evaluation protocol\. The SILVA molecular\-regression model reaches0\.3893±0\.01320\.3893\\pm 0\.0132validation MAE across five seeds \(Table[1](https://arxiv.org/html/2607.28989#S3.T1)\), improving on GCN, GAT, GIN, and GatedGCN within the same training and evaluation protocol while trailing PNA \(0\.3475±0\.01300\.3475\\pm 0\.0130\)\. The SILVA model uses bond\-type input through the molecular edge encoder, whereas the baselines in this matched protocol use only the shared node\-feature and graph\-topology setting\. The ZINC table is therefore a same\-split and same\-budget comparison\. It does not isolate parameter count or bond\-type edge features from the SILVA interaction field\.

Table 1:ZINC\-12k validation mean absolute error \(MAE\)\. Lower is better\. Values are mean±\\pmstandard deviation across independent seeds\. GCN denotes a graph convolutional network, GAT a graph attention network, GIN a graph isomorphism network, PNA principal neighbourhood aggregation, and GatedGCN a gated graph convolutional network\. All rows share the same ZINC split, training protocol, and hidden dimensiond=64d=64\. The SILVA row additionally uses bond\-type edge input, and the rows are not parameter\-matched\.
#### Citation\-network node classification\.

SILVA casts images, molecules, and citation graphs as interaction systems solved to a fixed point under one structured template \(Section[2\.2](https://arxiv.org/html/2607.28989#S2.SS2)\)\. This representation predicts that the global interaction term should be useful when a task needs information mixing beyond immediate local structure and neutral when local structure already contains the relevant signal\. Node classification on citation networks tests this regime directly, since these tasks are transductive, single\-graph, and trained on a very small labeled budget \(140/120/60 nodes for Cora/Citeseer/Pubmed respectively\[[49](https://arxiv.org/html/2607.28989#bib.bib41),[58](https://arxiv.org/html/2607.28989#bib.bib13)\]\)\. The standard citation\-network graph convolutional network \(GCN\) and graph attention network \(GAT\) baselines were introduced in this same semi\-supervised, citation\-edge aggregation setting\[[34](https://arxiv.org/html/2607.28989#bib.bib18),[53](https://arxiv.org/html/2607.28989#bib.bib10)\]\. These datasets therefore probe a local\-structure\-dominated regime for the additional global mechanism\.

Table[2](https://arxiv.org/html/2607.28989#S3.T2)compares graph convolutional network \(GCN\) and graph attention network \(GAT\) baselines \(no equilibrium solve, no global term, five seeds each\) against SILVA \(local term instantiated as graph attention over the real citation graph, global term as a scalar\-gated mean\-field broadcast\)\. On Cora, SILVA remains below both local graph baselines and shows larger seed\-to\-seed dispersion\. SILVA reaches70\.80±4\.86%70\.80\\pm 4\.86\\%against77\.04±1\.21%77\.04\\pm 1\.21\\%\(graph convolution\) and78\.48±0\.95%78\.48\\pm 0\.95\\%\(graph attention\)\. Citeseer shows the same ordering, with SILVA at55\.96±1\.98%55\.96\\pm 1\.98\\%against67\.32±1\.22%67\.32\\pm 1\.22\\%and67\.60±1\.36%67\.60\\pm 1\.36\\%\. Together, Cora and Citeseer indicate that, in small local\-structure\-dominated transductive regimes, the additional equilibrium and global machinery does not improve over a plain graph\-attention classifier with no fixed\-point solve\. The corresponding longer\-range test is CLUSTER\[[20](https://arxiv.org/html/2607.28989#bib.bib17)\]\. Published 4\-to\-16\-layer graph\-convolution baselines without positional encoding span approximately 47\.8–69\.0%, whereas gated graph convolution spans approximately 60\.4–73\.8% under the same protocol, as discussed in the CLUSTER result below\.

A direct channel\-attention global term, obtained by applying the vision attention structure to the node dimension of a full graph, collapses at graph scale to an exactly uniform attention distribution \(measured row\-entropy ratio 1\.0000\) and accuracy near chance\. The scalar\-gated mean\-field broadcast of Section[2\.6](https://arxiv.org/html/2607.28989#S2.SS6)removes this collapse by construction\. A bounded top\-kkattention variant uses genuine node\-to\-node attention restricted to each node’skknearest candidates\. Across five seeds it reaches69\.12±4\.95%69\.12\\pm 4\.95\\%on Cora and55\.80±3\.53%55\.80\\pm 3\.53\\%on Citeseer, within the same range as the mean\-field variant\.

Table 2:Citation\-network node\-classification validation accuracy\. Higher is better\. Values are mean±\\pmstandard deviation acrossn=5n=5independent seeds on the public train/validation/test splits\[[58](https://arxiv.org/html/2607.28989#bib.bib13)\]matching Table 1 of\[[53](https://arxiv.org/html/2607.28989#bib.bib10)\]\. GCN denotes graph convolutional network and GAT denotes graph attention network\. Pubmed SILVA is shown per seed, reflecting bimodal outcomes\.*Note\. Pubmed SILVA lists three well\-converged seeds on the first line and two collapsed seeds on the second\.*

Term\-level attribution\.A seven\-variant ablation study on Cora and Citeseer \(Figure[5](https://arxiv.org/html/2607.28989#S3.F5), five seeds each\) sharpens this finding\. Figure[5](https://arxiv.org/html/2607.28989#S3.F5)a shows that on Cora, removing the local term is the only ablation that causes a large accuracy drop \(70\.8%→53\.0%70\.8\\%\\to 53\.0\\%\), whereas removing the global term slightly exceeds the full model \(72\.7%72\.7\\%vs\.70\.8%70\.8\\%\)\. Figure[5](https://arxiv.org/html/2607.28989#S3.F5)b repeats the same pattern on Citeseer\. Removing the local term again collapses accuracy \(56\.0%→48\.5%56\.0\\%\\to 48\.5\\%\), whereas the no\-global model again slightly exceeds the full model \(57\.9%57\.9\\%vs\.56\.0%56\.0\\%\)\. The static global matrix, equal fast/slow timescales, three\-layer stack, and top\-kkmechanism all remain near the full model on these two datasets\. The two panels therefore attribute the citation\-network limitation to the global term not adding useful information, not to the absence of local graph processing\.

On Pubmed, the evaluated variants do not form a single summary pattern and are therefore shown per seed\. Figure[5](https://arxiv.org/html/2607.28989#S3.F5)c separates the well\-converged seeds from collapsed seeds visually rather than compressing them into a single mean\. Four variants, the full model, no\-global, static, and same\-alpha, split across their five seeds into a well\-converged band \(70–78% accuracy\) and a collapsed band \(19–56% accuracy\)\. The remaining evaluated variants show different behavior\. The three\-layer stack collapses uniformly, all five seeds landing in the 39–42% range with none reaching the well\-converged band\. The no\-local \(global\-only\) variant instead converges tightly and unimodally in an intermediate 62–68% range and falls in neither band\. The Pubmed panel therefore reveals a seed\-dependent convergence bifurcation, not ordinary statistical noise around one operating point\. The bimodal cases indicate a convergence instability specific to Pubmed’s substantially larger graph \(N=19,717N=19\{,\}717vs\. Cora’sN=2,708N=2\{,\}708\), plausibly connected to the spectral\-radius framing of Section[2\.10](https://arxiv.org/html/2607.28989#S2.SS10)\. The uniformly\-collapsed and uniformly\-converged exceptions indicate this instability is not a fixed property of the Pubmed graph alone, but depends on which interaction terms are active\.

![Refer to caption](https://arxiv.org/html/2607.28989v1/figures/fig_silva_ablation_matrix_real.png)Figure 5:Citation\-network ablations\. \(a\) Cora\. \(b\) Citeseer\. Filled markers and whiskers show mean±\\pmstandard deviation across five seeds, with individual seeds overlaid as open markers\. \(c\) Individual Pubmed seed values for evaluated variants\.
#### Long\-range node classification on CLUSTER\.

On CLUSTER\[[20](https://arxiv.org/html/2607.28989#bib.bib17)\], the controlled comparison is between two baselines trained under the same data split and evaluation protocol, but not matched parameter counts\. A graph convolutional network \(GCN\) reaches53\.22±0\.50%53\.22\\pm 0\.50\\%across five seeds, and a gated graph convolutional network \(GatedGCN\) reaches58\.75±0\.23%58\.75\\pm 0\.23\\%across five seeds\. These baselines establish a consistent comparison gap in the direction the benchmark is designed to produce\. Dwivedi et al\.\[[20](https://arxiv.org/html/2607.28989#bib.bib17)\]report two distinct graph\-convolution baselines under this same protocol, differing in normalization convention\. They report a47\.847\.8–69\.0%69\.0\\%range \(4\-to\-16 layers, their “GCN” row\) and a53\.453\.4–68\.5%68\.5\\%range \(their “vanilla GCN” row\)\. The graph\-convolution baseline sits almost exactly on the latter’s 4\-layer point \(53\.22%53\.22\\%vs\.53\.45%53\.45\\%\), indicating that the configuration follows that convention\. For gated graph convolution, their plain \(no positional encoding\) range is60\.460\.4–73\.8%73\.8\\%\. A positional\-encoding\-augmented variant reaches76\.1%76\.1\\%at 16 layers but is architecturally distinct from the plain 4\-layer gated\-convolution baseline used here\. Across those published protocol and hyperparameter settings, the graph\-convolution baseline sits within its matching published range, and the gated graph\-convolution result sits just below its published range while preserving the same qualitative ordering\. Because CLUSTER is designed to reward a longer\-reach mechanism, it directly tests whether the full local\+global mechanism can help under longer\-range demands, complementing the local\-structure\-dominated regime identified in the citation\-network analysis above\. The SILVA local\+global mechanism, trained under the identical 50\-epoch budget as the two baselines above, reaches73\.04±0\.60%73\.04\\pm 0\.60\\%validation accuracy across five independent seeds\. It exceeds the same\-budget graph\-convolution rerun by19\.8219\.82percentage points and the same\-budget gated graph\-convolution rerun by14\.2914\.29percentage points\. A four\-seed ablation that removes the global interaction term while keeping the local interaction term reaches67\.55±0\.84%67\.55\\pm 0\.84\\%, isolating a5\.495\.49percentage\-point aggregate difference associated with the global term\. The CLUSTER benchmark is the evaluated setting in which the global interaction term is load\-bearing\.

Table 3:CLUSTER node\-classification validation accuracy over synthetic graph instances\. Values are mean±\\pmstandard deviation across independent seeds\. GCN denotes graph convolutional network and GatedGCN denotes gated graph convolutional network\.
#### Stacked equilibrium hierarchy\.

The two\-layer fast/slow stack \(Eq\. \([1](https://arxiv.org/html/2607.28989#S2.E1)\)\) used in the CIFAR\-10, ZINC, citation\-network, and CLUSTER experiments is the shallowest instance of a more general architectural idea \(Section[2](https://arxiv.org/html/2607.28989#S2)\)\. It is a hierarchy of equilibria at progressively slower timescales, each still solving the full interaction structure\. This composition differs from previously evaluated recurrent hierarchies by making every rung an implicit local\+global equilibrium\. It is a multi\-timescale hierarchy\[[35](https://arxiv.org/html/2607.28989#bib.bib14),[21](https://arxiv.org/html/2607.28989#bib.bib15)\]in which every rung is an implicit fixed\-point solve with its own local\+global interaction structure, not a single fixed\-depth recurrent step\.

A first three\-layer instance \(α=0\.5,0\.35,0\.2\\alpha=0\.5,0\.35,0\.2\) was evaluated over five independent seeds on Cora and Citeseer under the same data split and training budget as Table[2](https://arxiv.org/html/2607.28989#S3.T2)\. It reaches69\.76±2\.43%69\.76\\pm 2\.43\\%on Cora, within the two\-layer variant’s70\.80±4\.86%70\.80\\pm 4\.86\\%spread, and53\.96±3\.88%53\.96\\pm 3\.88\\%on Citeseer, below the two\-layer variant’s55\.96±1\.98%55\.96\\pm 1\.98\\%\. The seed\-0 traces show training loss continuing to fall while validation loss climbs past approximately epoch 30\. This overfitting signature is consistent with the citation\-network analysis above, which identifies citation networks as a regime in which additional equilibrium capacity is not expected to help\.

Depth can also be added*within*a single equilibrium layer’s own solver iteration by repeating the local term’s message\-passing steplocal​\_​depth\\mathrm\{local\\\_depth\}times, weight\-tied, before combining with the global term\. This differs from adding separate stacked equilibria as above\. A Cora sweep with one seed per local\-depth\-2 variant \(six variants atlocal​\_​depth=2\\mathrm\{local\\\_depth\}=2, plus the full model atlocal​\_​depth=3\\mathrm\{local\\\_depth\}=3\) found accuracy decreasing at every variant tested\. The deltas compare the five\-seed means atlocal​\_​depth=1\\mathrm\{local\\\_depth\}=1with the single completed seed atlocal​\_​depth=2\\mathrm\{local\\\_depth\}=2\. The largest decreases occur for the top\-kkglobal\-term variant \(69\.1%→55\.4%69\.1\\%\\to 55\.4\\%\) and the combination with three\-layer stacking \(69\.8%→63\.4%69\.8\\%\\to 63\.4\\%\), consistent with the over\-smoothing failure mode identified for deeper graph convolutional propagation\[[39](https://arxiv.org/html/2607.28989#bib.bib46)\]\. Added local depth is therefore unfavorable in this citation\-network setting\.

#### Comparison with multi\-domain latent\-attention architectures\.

Several architectures target multiple, heterogeneous data domains with one mechanism\. Perceiver IO\[[30](https://arxiv.org/html/2607.28989#bib.bib38)\]is a single cross\-attention architecture demonstrated across images \(ImageNet, 84\.5% top\-1\), language, optical flow, multimodal audio\-video tasks, and symbolic game state, by mapping arbitrary inputs into a fixed\-size latent array via cross\-attention and querying that latent flexibly for arbitrary outputs\. Graph Perceiver IO\[[5](https://arxiv.org/html/2607.28989#bib.bib39)\]extends the same architecture to graph\-structured data and to joint image\-graph few\-shot classification\.

Perceiver IO shares no benchmark with this evaluation suite, so the comparison is architectural\. Graph Perceiver IO overlaps with SILVA on the same three citation\-network datasets evaluated in the citation\-network analysis above under the public split and 100\-repeat protocol\. It reports Cora83\.9±0\.6%83\.9\\pm 0\.6\\%, Citeseer70\.1±1\.0%70\.1\\pm 1\.0\\%, Pubmed79\.9±0\.4%79\.9\\pm 0\.4\\%\[[5](https://arxiv.org/html/2607.28989#bib.bib39)\], exceeding both the graph\-convolution/graph\-attention baselines and the SILVA values in Table[2](https://arxiv.org/html/2607.28989#S3.T2)\. Its reported baselines on these splits are drawn from literature\-tuned values rather than a matched rerun under the present protocol\. The accuracy gap therefore reflects both architectural and protocol differences\.

On citation networks, Graph Perceiver IO is ahead under its reported protocol\. The citation\-network analysis above identifies citation networks as a local\-structure\-dominated regime for which the additional global machinery is weakly expressed\. SILVA contributes a falsifiable dynamical account of the mechanism itself\. This account includes a monitored spectral radius that tracks instability \(Section[2\.7](https://arxiv.org/html/2607.28989#S2.SS7)\), a solver\-state energy diagnostic and local descent surrogate \(Section[2\.8](https://arxiv.org/html/2607.28989#S2.SS8)\), a finite path\-sum decomposition for the local linearization \(Section[2\.10](https://arxiv.org/html/2607.28989#S2.SS10)\), and an explicit, graph\-structured local/global interaction split in place of an opaque learned latent bottleneck\.

#### Local–global interaction ablations\.

A four\-arm ablation on MNIST \(full local\+global model, local only, global only, and a stimulus\-only variant with neither interaction term, five seeds each, matched hidden width and training budget\) reaches96\.85±0\.13%96\.85\\pm 0\.13\\%,96\.83±0\.16%96\.83\\pm 0\.16\\%,96\.95±0\.17%96\.95\\pm 0\.17\\%, and96\.95±0\.22%96\.95\\pm 0\.22\\%validation accuracy respectively\. The four arms overlap within seed variation, with every pairwise difference within 1–2 standard deviations\. A direct diagnostic confirms that the local and global terms alter the computation when enabled\. The maximum output difference between the full and stimulus\-only variants is approximately 1\.01 at identical initialization and input\. The near\-identical accuracy therefore reflects this task and parameter budget rather than absence of a computational difference between variants\. Figure[6](https://arxiv.org/html/2607.28989#S3.F6)e shows the full comparison\. Together with the citation\-network result above, this extends the same finding to a second domain\. At this capacity budget, MNIST classification is easy enough that the additional local/global interaction machinery does not measurably improve on the stimulus\-only equilibrium variant, consistent with the representation’s own premise that these terms are expected to help specifically when a task needs richer structure than a single injected stimulus already provides\. The remaining panels of Figure[6](https://arxiv.org/html/2607.28989#S3.F6)provide an operator\-level MNIST diagnostic from the adjoint\-trained equilibrium model with local attention and a learned static global matrix\. Figure[6](https://arxiv.org/html/2607.28989#S3.F6)a shows the learned hidden\-neuronkk\-nearest\-neighbor graph used by the local term\. Its sparse binary structure fixes which state coordinates may exchange local messages\. Figure[6](https://arxiv.org/html/2607.28989#S3.F6)b shows the corresponding head\-averaged attention weights for one sample\. The weights are nonuniform over the allowed edges, showing that the local term is an input\-dependent interaction operator rather than a fixed averaging rule\. Figure[6](https://arxiv.org/html/2607.28989#S3.F6)c shows the static global matrix\. Its dense signed entries complement the sparse local graph by allowing every coordinate to contribute to every other coordinate through the global channel\. Figure[6](https://arxiv.org/html/2607.28989#S3.F6)d compares the diagonal and off\-diagonal entry distributions of this matrix\. The larger positive diagonal mean indicates an emergent self\-feedback component inside the learned global interaction\. Figure[6](https://arxiv.org/html/2607.28989#S3.F6)f closes the diagnostic by showing that the Picard solver reaches tolerance in a small finite number of iterations on the tested samples\. Together, these panels show the concrete objects represented by the SILVA decomposition: sparse local support, learned local weights, dense global coupling, self\-interaction in the global operator, and finite solver effort\.

![Refer to caption](https://arxiv.org/html/2607.28989v1/figures/fig_interaction_structure_real.png)Figure 6:MNIST local\-global diagnostics\. \(a–d,f\) Operator and solver diagnostics from an adjoint\-trained equilibrium model with local attention and a learned static global matrix\. \(a\) Learnedkk\-nearest\-neighbor graph in hidden\-neuron space for one test sample\. \(b\) Head\-averaged graph\-attention weights for the same sample and layer\. \(c\) Learned static global interaction matrix\. \(d\) Distributions of diagonal and off\-diagonal entries of the static global matrix, with filled markers denoting means\. \(e\) Four\-arm MNIST ablation across full, local\-only, global\-only, and stimulus\-only variants, shown as mean±\\pmstandard deviation with individual seeds overlaid\. \(f\) Picard iterations required to reach residual tolerance10−310^\{\-3\}over 64 test samples\.Parameter\-matched computational cost\.The ablation above holds the equilibrium solve fixed and varies only which interaction terms feed it\. A single\-hidden\-layer feedforward network with no fixed\-point iteration, no local or global interaction term, and a parameter count matched to the full equilibrium model to within 0\.4% \(63,610 vs\. 63,370 parameters\) reaches97\.29±0\.23%97\.29\\pm 0\.23\\%validation accuracy across five seeds, matching or marginally exceeding the full model’s96\.85±0\.13%96\.85\\pm 0\.13\\%, at 2\.5\- to 6\.5\-fold lower wall\-clock cost per epoch \(12–29s vs\. 72–80s, measured under identical conditions\)\. Taken together with the four\-arm ablation above, this indicates that at the MNIST task and parameter budget tested here, neither the choice of interaction term nor the equilibrium mechanism itself is measurably earning its additional computational cost over the feedforward baseline\. The result identifies a task regime in which the SILVA mechanism is unnecessary at the tested capacity\. Together with the citation\-network result above, it supports the interpretation that the representation’s value lies in tasks that need richer structure than a single\-shot stimulus already provides, and that MNIST classification at this capacity is not such a task\.

Linear decodability of frozen representations\.Classification accuracy requires the minimal statistic sufficient to separate ten classes, and has no reason to reward any additional structure a representation happens to preserve\. Two linear probes test this on the trained frozen representations, with no reference to digit identity\. The probes are a four\-way rotation probe \(0/90/180/270 degrees, exact array rotation, no interpolation\) and a strictly linear pixel\-space reconstruction decoder, both trained only on top of the frozen representation with the backbone untouched\. On both probes the plain feedforward baseline’s representation is more linearly decodable than the full equilibrium model’s\. Rotation probe accuracy80\.3%80\.3\\%vs\.74\.6%74\.6\\%\(chance25%25\\%\), and reconstruction mean squared error0\.02840\.0284vs\.0\.03180\.0318\(lower is better\)\. The comparison differs in representation width \(80 vs\. 64 dimensions, matched on parameter count rather than dimensionality, since the feedforward network concentrates its budget into one wide layer where the equilibrium model spreads it across several weight matrices\), so some of this gap may reflect probe capacity rather than information content\. These probes do not support the hypothesis that the interaction terms incidentally preserve more usable structure than a plain baseline, on the two probes tested here\.

Optimization diagnostics\.Mean validation accuracy across seeds can hide whether a flat ablation reflects task structure or optimization failure\. An optimization explanation would predict a landscape with variable\-quality local optima, the same signature already observed for the SILVA node\-classification variants on Pubmed in the citation\-network analysis above, where seeds split bimodally between well\-converged and collapsed outcomes\. Two checks distinguish these on the MNIST ablation matrix specifically\. First, final training losses across the four arms overlap tightly and remain close to zero \(final training loss0\.01470\.0147–0\.02050\.0205, no arm systematically higher\), so no arm is failing to fit the training data\. The flat result is not an underfitting symptom\. Second, comparing best\-of\-five rather than mean\-of\-five seeds does not favor the full model\. Its best seed \(96\.98%96\.98\\%\) is in fact the lowest best\-of\-five across all four arms, while the stimulus\-only variant’s best seed \(97\.20%97\.20\\%\) is the highest\. If some initializations of the full model reliably found a materially better configuration, its best\-of\-five would be expected to pull ahead of the simpler arms\. The observed best\-of\-five results therefore favor a genuine performance limit on this benchmark over an optimization\-landscape explanation, in contrast to the bimodal instability observed on Pubmed\. Truncated backpropagation through the solver, compounded by stacking, may still produce a harder\-to\-navigate landscape differently across tasks\.

#### Stability dynamics\.

Figure[7](https://arxiv.org/html/2607.28989#S3.F7)traces a training instability in the adjoint\-trained MNIST diagnostic described above\. In Figure[7](https://arxiv.org/html/2607.28989#S3.F7)a, validation accuracy is 89\.9% at the first saved checkpoint, then collapses to 14\.2% at the following epoch, recovering only partially to 53\.1% over six subsequent epochs\. The paired loss curve remains elevated after the collapse, indicating that the model did not simply move to a different high\-accuracy solution with changed calibration\. The spectral radius \(Algorithm[2](https://arxiv.org/html/2607.28989#alg2)\) wasρL​1=0\.67\\rho\_\{L1\}=0\.67,ρL​2=0\.63\\rho\_\{L2\}=0\.63\(below the localρ=1\\rho=1boundary\) immediately before the instability andρL​1=1\.15\\rho\_\{L1\}=1\.15,ρL​2=1\.89\\rho\_\{L2\}=1\.89\(above that boundary\) at the final epoch\. Figure[7](https://arxiv.org/html/2607.28989#S3.F7)b compares those measured checkpoints with layer\-wise random\-initialization values from five seeds for the same configuration\. The initialization values lie well below theρ=1\\rho=1boundary, whereas the trained checkpoint moves one layer moderately and the other layer strongly beyond it\. The instability is therefore associated with the training trajectory moving the local Jacobian beyond the measuredρ=1\\rho=1stability boundary, not with this configuration starting beyond that boundary at initialization\.

![Refer to caption](https://arxiv.org/html/2607.28989v1/figures/fig_edge_of_chaos_real.png)Figure 7:Loss of local stability in one adjoint\-trained MNIST diagnostic trajectory\. \(a\) Validation loss and accuracy over training epochs\. Dotted vertical lines mark the checkpoints at which the layer\-wise spectral radii were measured\. \(b\) Spectral radii at random initialization and at the two measured checkpoints\. Circles denote the first equilibrium layer and squares denote the second\. Open gray markers show the five random\-initialization seeds for each layer, and filled gray markers with whiskers show mean±\\pmstandard deviation\. Green checkpoint markers indicateρ<1\\rho<1, red checkpoint markers indicateρ\>1\\rho\>1, and the dotted vertical line marksρ=1\\rho=1\.
#### Epoch\-resolved near\-boundary trajectory\.

The configuration\-space analysis below identifies its largest\-ρ\\rhocell at random initialization \(multi\-head attention local term, node\-to\-node global\-attention term, meanρ=1\.27\\rho=1\.27across five seeds\)\. One instance of that configuration was trained for twenty epochs withρ\\rhoestimated at every epoch for both layers \(Figure[8](https://arxiv.org/html/2607.28989#S3.F8)\)\. This trajectory began with both layers aboveρ=1\\rho=1\(ρL​1=1\.19\\rho\_\{L1\}=1\.19,ρL​2=1\.05\\rho\_\{L2\}=1\.05at initialization\), the opposite starting condition from the instability traced above\. Training nonetheless converged normally\. Validation accuracy rose from 55\.3% \(epoch 1\) to a peak of 92\.3% \(epoch 18\), ending at 91\.0% \(epoch 20\), while validation loss fell from 1\.41 to 0\.30 over the same trajectory\.ρ\\rhofor both layers fell below 1 within the first three epochs and then oscillated in a narrow band \(0\.78–1\.10\) for the remainder of training\. It crossed back above 1 intermittently, in three of twenty epochs for each layer, including one epoch \(epoch 9\) where both layers exceeded 1 simultaneously \(ρL​1=1\.10\\rho\_\{L1\}=1\.10,ρL​2=1\.07\\rho\_\{L2\}=1\.07\), without any accompanying drop in validation accuracy\.

The epoch\-resolved trace shows that crossingρ=1\\rho=1is not sufficient, by itself, to trigger the kind of collapse traced above\. In Figure[8](https://arxiv.org/html/2607.28989#S3.F8)a, both layers start at or above the local stability boundary, fall below it within the first few epochs, and make only brief excursions above 1 later in training\. The marked epoch 9 is the strongest joint excursion, but both layers remain close to the boundary \(ρL​1=1\.10\\rho\_\{L1\}=1\.10,ρL​2=1\.07\\rho\_\{L2\}=1\.07\)\. In Figure[8](https://arxiv.org/html/2607.28989#S3.F8)b, validation accuracy continues rising through the same epoch, and validation loss continues its downward trend over the surrounding epochs\. Together, the two trajectories distinguish transient threshold crossings from sustained large excursions above the measured boundary\. The collapsed trajectory was substantially and persistently above the boundary in the second layer \(ρL​2=1\.89\\rho\_\{L2\}=1\.89\), whereas the probe’s excursions above 1 were small and self\-corrected within one to two epochs\. The evidence is consistent with brief, small excursions nearρ=1\\rho=1being tolerated, while a large, sustained departure from the stability boundary marks the observed instability\. This refines the spectral\-radius diagnostic within the probed configuration\.

![Refer to caption](https://arxiv.org/html/2607.28989v1/figures/fig_instability_probe_real.png)Figure 8:Epoch\-resolved spectral\-radius probe for the configuration\-space grid’s largest\-ρ\\rhocell at initialization\. \(a\) Layer\-wise spectral radius during training\. Epoch 0 denotes random initialization, the horizontal dotted line marksρ=1\\rho=1, and the red vertical line marks epoch 9, where both layers briefly exceedρ=1\\rho=1\. \(b\) Validation accuracy and validation loss for the same trajectory\.
#### Configuration\-space stability\.

Figure[9](https://arxiv.org/html/2607.28989#S3.F9)mapsρ𝒯\\rho\_\{\\mathcal\{T\}\}at random initialization across the full grid of local\- and global\-interaction variants, with five seeds per cell\. In Figure[9](https://arxiv.org/html/2607.28989#S3.F9)a, color encodes the across\-seed mean and each cell reports mean±\\pmstandard deviation\. The no\-global row remains belowρ𝒯=1\\rho\_\{\\mathcal\{T\}\}=1throughout the local\-interaction choices, including the configuration underlying the instability traced above \(ρ𝒯=0\.41±0\.08\\rho\_\{\\mathcal\{T\}\}=0\.41\\pm 0\.08\)\. Adding simple or multi\-head global attention raises the gain substantially\. The graph\-neural local column crosses the local stability boundary for both global\-attention variants \(ρ𝒯=1\.26±0\.06\\rho\_\{\\mathcal\{T\}\}=1\.26\\pm 0\.06and1\.27±0\.081\.27\\pm 0\.08\), while the GAT local column moves close to the boundary with larger seed variation \(0\.96±0\.210\.96\\pm 0\.21\)\. Figure[9](https://arxiv.org/html/2607.28989#S3.F9)b isolates that GAT local column and shows the five individual seed estimates for each global\-interaction mode\. Changing only the global interaction shifts the values from a safely below\-boundary regime to a near\-boundary regime in which some seeds exceedρ𝒯=1\\rho\_\{\\mathcal\{T\}\}=1\. Local stability is therefore not uniform across this architecture’s own configuration space and cannot be assumed by analogy from one working configuration to a nearby one\.

![Refer to caption](https://arxiv.org/html/2607.28989v1/figures/fig_configspace_rho.png)Figure 9:Configuration\-space map ofρ𝒯\\rho\_\{\\mathcal\{T\}\}at random initialization\. \(a\) Meanρ𝒯\\rho\_\{\\mathcal\{T\}\}for each pairing of global and local interaction modes\. Cells report mean±\\pmstandard deviation overn=5n=5seeds\. \(b\) Individual seed estimates in the GAT local\-interaction column, with diamonds and whiskers showing mean±\\pmstandard deviation\. Dotted lines markρ𝒯=1\\rho\_\{\\mathcal\{T\}\}=1\.
#### Robustness versus spectral radius\.

Section[2\.10](https://arxiv.org/html/2607.28989#S2.SS10)and Figure[7](https://arxiv.org/html/2607.28989#S3.F7)connectρ\\rhoapproaching 1 to loss of local linear stability and to one observed training instability\. A broader test evaluates whether a converged equilibrium closer toρ=1\\rho=1degrades faster under input noise, with no additional training\.ρ\\rhowas estimated \(Algorithm[2](https://arxiv.org/html/2607.28989#alg2)\) at each of twenty independently\-trained MNIST ablation\-matrix models’ own converged fixed point, and each model was evaluated with increasing fractions of input pixels replaced by noise \(Figure[10](https://arxiv.org/html/2607.28989#S3.F10)\)\.

There is no measurable correlation \(Pearsonr=−0\.026r=\-0\.026,n=20n=20\)\. In Figure[10](https://arxiv.org/html/2607.28989#S3.F10)a, the mean accuracy curves for the four ablation families degrade smoothly as the fraction of corrupted pixels increases, and their seed\-to\-seed bands overlap over most noise levels\. No ablation family separates into a distinctly more robust or more fragile group\. In Figure[10](https://arxiv.org/html/2607.28989#S3.F10)b, the measured fixed\-point spectral radii do not order the degradation slopes, and individual checkpoints from the four interaction variants remain intermingled\. The estimator itself is calibrated by the stimulus\-only variant, where both interaction terms are disabled\. This variant measuresρ=0\.75\\rho=0\.75for all five seeds\. It matches the analytically exact value\(1−α\)\(1\-\\alpha\)forα=0\.25\\alpha=0\.25that this configuration’s update provably reduces to, since with both interaction terms off the stimulus does not depend on the state at all and the true Jacobian is exactly\(1−α\)​I\(1\-\\alpha\)I\. This null result covers one perturbation type, one dataset, and aρ\\rhorange \(0\.75–0\.97\) below the instability threshold observed in training\. The stability\-dynamics traces provide the complementary epoch\-resolved view\. A system deliberately tracked through the transition acrossρ=1\\rho=1during training, rather than evaluated only at a safely below\-boundary operating point, shows that transient, small crossings ofρ=1\\rho=1are tolerated without collapse, while the unstable trajectory analyzed above involved a much larger and sustained departure from the local stability boundary\. Both results support a stability\-diagnostic role forρ\\rhoin SILVA\. Chu et al\.\[[17](https://arxiv.org/html/2607.28989#bib.bib43)\]report that*imposing*Lyapunov stability as an explicit constraint, rather than observing where a network’sρ\\rhohappens to fall, yields a robustness benefit for deep equilibrium models\.

![Refer to caption](https://arxiv.org/html/2607.28989v1/figures/fig_robustness_vs_rho_real.png)Figure 10:Robustness and spectral radius for twenty independently\-trained MNIST ablation\-matrix models\. \(a\) Accuracy as increasing fractions of input pixels are replaced with noise\. Lines show means and bands show standard deviations across five seeds per ablation family\. \(b\) Measuredρ𝒯\\rho\_\{\\mathcal\{T\}\}at each model’s converged fixed point versus the fitted accuracy\-degradation slope\. Slopes are reported in percentage points per unit noise fraction, and the gray line shows the least\-squares fit\.

## 4Discussion

The experiments support SILVA as a structured implicit\-layer framework in which input stimulus, local interaction, global context, and self\-persistence are separated as trainable components of the same equilibrium computation\. The model reads out the settled response of an interacting system driven by the input and relaxed by a damped solver\. This separation is the main architectural contribution\. It makes the stimulus entering the system, the routes through which that stimulus propagates, the damping that controls relaxation, and the vector attractor that is ultimately read out available for ablation, visualization, and stability analysis\.

The experiments establish three consequences of this formulation\. First, the same outer equilibrium template transfers across image classification, molecular property regression, citation\-network node classification, and long\-range node classification, while the local and global operators remain adapted to each domain\. Across these settings, the shared object is the decomposition into stimulus, local interaction, global interaction, and solver, with each domain specifying the state space and the corresponding interaction operators\. Second, the representation can be inspected through its own operators\. The MNIST results expose learned neighborhoods, attention weights, global interaction structure, solver residuals, and equilibrium\-state geometry\. The spectral\-radius and path\-sum derivations connect those measurements to the damped iteration rather than treating them as post\-hoc visualizations\. Third, stacked SILVA layers compose separately solved approximate vector attractor states, giving depth through a sequence of equilibrium responses rather than through a larger number of ordinary feedforward transformations\.

The empirical pattern is consistent with the intended role of the global term\. On MNIST and citation\-network node classification, the additional global interaction does not improve over a local\-only or no\-interaction baseline\. These are regimes in which local structure, task simplicity, or limited label budgets can already determine much of the solution\. On CLUSTER, where long\-range aggregation is central\[[20](https://arxiv.org/html/2607.28989#bib.bib17)\], the same decomposition identifies the opposite case\. The full local\+global mechanism reaches73\.04±0\.60%73\.04\\pm 0\.60\\%across five independent seeds, exceeding the same\-budget gated graph\-convolution rerun by14\.2914\.29percentage points\. A four\-seed no\-global ablation reaches67\.55±0\.84%67\.55\\pm 0\.84\\%, giving a5\.495\.49percentage\-point aggregate difference associated with the global term in the long\-range setting\. This contrast identifies the regime in which the extra mechanism contributes measurable signal\.

SILVA is also distinct from the related families that motivate it\. Deep equilibrium models and multiscale DEQs provide the fixed\-point view and the implicit\-gradient framework\[[6](https://arxiv.org/html/2607.28989#bib.bib5),[7](https://arxiv.org/html/2607.28989#bib.bib6)\]\. SILVA uses that view to expose the internal composition of the equilibrium update\. Graph neural networks and attention layers provide learned interaction operators\[[34](https://arxiv.org/html/2607.28989#bib.bib18),[53](https://arxiv.org/html/2607.28989#bib.bib10),[52](https://arxiv.org/html/2607.28989#bib.bib11)\]\. SILVA places those operators inside a repeated, damped fixed\-point field and separates the external stimulus from local and global interactions\. Hopfield and related attractor models provide useful language for settled states and collective dynamics\[[27](https://arxiv.org/html/2607.28989#bib.bib23)\], but SILVA does not require symmetric weights, a Hopfield energy, or biological interpretation\. Graph Perceiver IO remains a strong reference point for the same citation\-network benchmark family\[[5](https://arxiv.org/html/2607.28989#bib.bib39)\]\. SILVA contributes a complementary formulation, a general implicit interaction field whose operators and dynamics are explicitly defined and empirically interrogated\.

Several points delimit this interpretation\. Most reported experiments train the finite unrolled solver with truncated backpropagation through Picard iterations\. Only the adjoint\-trained MNIST diagnostic uses the GMRES adjoint described in the Supplementary Methods\. Accordingly, the principal reported results characterize weight\-tied iterative networks trained through finite solver trajectories, without an impliedO​\(1\)O\(1\)\-memory guarantee\. Even when an adjoint solve is used, it is computed numerically by GMRES to a fixed tolerance, so the resulting gradient is an approximation to an implicit\-function\-theorem gradient\. This is a known source of instability in DEQs and motivates recent reversible formulations that compute exact gradients directly\[[40](https://arxiv.org/html/2607.28989#bib.bib44)\]\.

The discrete neighborhoods used by the vision hidden\-channel graph and the node\-classification top\-kkglobal\-attention variant introduce an additional qualification\. The selected graph is recomputed from the current hidden state, and top\-kkselection is not everywhere differentiable\. Implicit differentiation for those cases treats the selected graph as frozen atz∗z^\{\\ast\}\. This qualification does not apply to ZINC, where the local graph is the fixed molecular bond graph\. The Pubmed instability and the three\-layer stacking result also have wider uncertainty than the five\-seed aggregate comparisons\. Finally, the Banach fixed\-point guarantee underlying Eq\. \([13](https://arxiv.org/html/2607.28989#S2.E13)\) and the adjoint derivation of Eqs\. \([75](https://arxiv.org/html/2607.28989#S2.E75)\)–\([77](https://arxiv.org/html/2607.28989#S2.E77)\) apply to the nonlinear update as stated\. The path\-sum decomposition in Eq\. \([59](https://arxiv.org/html/2607.28989#S2.E59)\) and the local descent surrogate in Eqs\. \([45](https://arxiv.org/html/2607.28989#S2.E45)\)– \([49](https://arxiv.org/html/2607.28989#S2.E49)\) are exact for the specified local linearization, withtanh⁡\(z\)≈z\\tanh\(z\)\\approx zorhlocal\+hglobalh\_\{\\mathrm\{local\}\}\+h\_\{\\mathrm\{global\}\}held fixed for the gradient\. They explain the observed solver dynamics through a controlled approximation rather than proving an exact identity for the full nonlinear system\.

Taken together, the results establish SILVA as a structured way to build and study implicit neural layers\. The method turns within\-layer interaction into an explicit modeling object, makes local and global influence separable inside the solver, and provides diagnostics that connect performance to the dynamics of the equilibrium computation\. The contribution is therefore methodological as well as empirical\. SILVA defines a learnable interaction field whose settled vector state can be trained, inspected, stacked, and compared across domains\.

## Methods summary

Full architectural, training, and evaluation detail is given in Section[2](https://arxiv.org/html/2607.28989#S2)\. The baseline families are graph convolutional networks \(GCNs\), graph attention networks \(GATs\), graph isomorphism networks \(GINs\), principal neighbourhood aggregation \(PNA\), and gated graph convolutional networks \(GatedGCNs\)\. They were trained under the same dataset splits, training and evaluation protocol, and hidden\-width setting as the corresponding SILVA configuration in each experiment\. They are not parameter\-matched unless explicitly stated\. Means and standard deviations are computed over five independent seeds unless stated otherwise\. Spectral radius \(Algorithm[2](https://arxiv.org/html/2607.28989#alg2)\) is estimated by power iteration on the vector–Jacobian product of the relevant equilibrium layer’s update map, evaluated either at random initialization in the configuration\-space analysis or at the converged fixed point of a trained model in the robustness analysis, as stated in each case\.

Table[4](https://arxiv.org/html/2607.28989#Sx1.T4)lists the training hyperparameters for SILVA in each experiment\. Baseline hyperparameters follow the corresponding baseline family under the same training scale, rather than being forced to identical parameter counts or identical optimizer settings\.

Table 4:Training hyperparameters for SILVA by experiment\. Hereddis hidden dimension,η\\etais learning rate,BBis batch size, andkkis localkk\-nearest\-neighbor size for vision or the number of attention heads for graph experiments\. Patience and gradient\-clip norm are as defined in Section[2](https://arxiv.org/html/2607.28989#S2)’s Training paragraph\.aCora/Citeseer/Pubmed are single\-graph transductive tasks\. Batch size is not meaningful and is fixed at 1 \(the whole graph\) under the evaluation protocol\. Damping is architecture\-specific\. MNIST uses one equilibrium layer withα=0\.25\\alpha=0\.25and 20 Picard iterations\. CIFAR\-10 and ZINC use two equilibrium layers withα1=0\.5\\alpha\_\{1\}=0\.5,α2=0\.2\\alpha\_\{2\}=0\.2, and 20 Picard iterations per layer\. Node\-classification and graph\-benchmark experiments use two equilibrium layers withα1=0\.5\\alpha\_\{1\}=0\.5,α2=0\.2\\alpha\_\{2\}=0\.2, and 15 Picard iterations per layer unless an ablation changes the stack, damping, or local\-interaction repetition depth\.

## Supplementary Methods for scalar and stacked equilibrium derivations

### Three scalar nodes

This scalar example is the minimal computational version of the interacting equilibrium layer\. Let the input bex∈ℝPx\\in\\mathbb\{R\}^\{P\}, and let one hidden equilibrium layer contain three nodes\. The same receiver\-first convention is used here\. The entryai​ja\_\{ij\}is the influence from source nodejjto receiver nodeii\. The feed\-forward stimulus for nodei∈\{1,2,3\}i\\in\\\{1,2,3\\\}is

mi=∑r=1PWi​r​xr\+ci,m\_\{i\}=\\sum\_\{r=1\}^\{P\}W\_\{ir\}x\_\{r\}\+c\_\{i\},\(91\)Hererrindexes input features,Wi​rW\_\{ir\}is the weight from input coordinaterrto hidden nodeii, andcic\_\{i\}is the bias for that node\. Equivalently, in matrix form,

m\\displaystyle m=\[W11⋯W1​PW21⋯W2​PW31⋯W3​P\]​x\+\[c1c2c3\],\\displaystyle=\\begin\{bmatrix\}W\_\{11\}&\\cdots&W\_\{1P\}\\\\ W\_\{21\}&\\cdots&W\_\{2P\}\\\\ W\_\{31\}&\\cdots&W\_\{3P\}\\end\{bmatrix\}x\+\\begin\{bmatrix\}c\_\{1\}\\\\ c\_\{2\}\\\\ c\_\{3\}\\end\{bmatrix\},\(92\)W\\displaystyle W∈ℝ3×P,\\displaystyle\\in\\mathbb\{R\}^\{3\\times P\},\(93\)c,m\\displaystyle c,m∈ℝ3\.\\displaystyle\\in\\mathbb\{R\}^\{3\}\.\(94\)HereWWstacks the three stimulus rows,ccstacks the three biases, andm=\(m1,m2,m3\)⊤m=\(m\_\{1\},m\_\{2\},m\_\{3\}\)^\{\\top\}is the stimulus vector injected into the equilibrium layer\. The intra\-layer interactions are

b1\\displaystyle b\_\{1\}=m1\+a11​ϕ​\(b1\)\+a12​ϕ​\(b2\)\+a13​ϕ​\(b3\),\\displaystyle=m\_\{1\}\+a\_\{11\}\\phi\(b\_\{1\}\)\+a\_\{12\}\\phi\(b\_\{2\}\)\+a\_\{13\}\\phi\(b\_\{3\}\),\(95\)b2\\displaystyle b\_\{2\}=m2\+a21​ϕ​\(b1\)\+a22​ϕ​\(b2\)\+a23​ϕ​\(b3\),\\displaystyle=m\_\{2\}\+a\_\{21\}\\phi\(b\_\{1\}\)\+a\_\{22\}\\phi\(b\_\{2\}\)\+a\_\{23\}\\phi\(b\_\{3\}\),\(96\)b3\\displaystyle b\_\{3\}=m3\+a31​ϕ​\(b1\)\+a32​ϕ​\(b2\)\+a33​ϕ​\(b3\),\\displaystyle=m\_\{3\}\+a\_\{31\}\\phi\(b\_\{1\}\)\+a\_\{32\}\\phi\(b\_\{2\}\)\+a\_\{33\}\\phi\(b\_\{3\}\),\(97\)Herebib\_\{i\}is the pre\-activation state of nodeii,ϕ\\phiis the scalar activation applied elementwise, and eachai​ja\_\{ij\}multiplies the nonlinear signal sent from source nodejjto receiver nodeii\. The compact vector form isb=m\+A​ϕ​\(b\)b=m\+A\\phi\(b\)with

A\\displaystyle A=\[a11a12a13a21a22a23a31a32a33\],\\displaystyle=\\begin\{bmatrix\}a\_\{11\}&a\_\{12\}&a\_\{13\}\\\\ a\_\{21\}&a\_\{22\}&a\_\{23\}\\\\ a\_\{31\}&a\_\{32\}&a\_\{33\}\\end\{bmatrix\},\(98\)A\\displaystyle A∈ℝ3×3\.\\displaystyle\\in\\mathbb\{R\}^\{3\\times 3\}\.\(99\)HereAAcollects all pairwise interaction coefficients, including the diagonal self\-interaction coefficientsa11,a22,a33a\_\{11\},a\_\{22\},a\_\{33\}\. The nonlinear SILVA model solves this system by iteration, either as the undamped map

bk\+1=m\+A​ϕ​\(bk\),b\_\{k\+1\}=m\+A\\phi\(b\_\{k\}\),\(100\)or as the damped map

bk\+1=\(1−α\)​bk\+α​\[m\+A​ϕ​\(bk\)\]\.b\_\{k\+1\}=\(1\-\\alpha\)b\_\{k\}\+\\alpha\[m\+A\\phi\(b\_\{k\}\)\]\.\(101\)Herebkb\_\{k\}is the solver iterate afterkkupdates andα\\alphais the same damping coefficient used in Eq\. \([1](https://arxiv.org/html/2607.28989#S2.E1)\)\. The bracketed term is the undamped nonlinear update, and the prefactor\(1−α\)​bk\(1\-\\alpha\)b\_\{k\}carries forward part of the previous state\. In the linear special caseϕ​\(b\)=b\\phi\(b\)=b, the exact solution is

b∗=\(I−A\)−1​mb^\{\\ast\}=\(I\-A\)^\{\-1\}m\(102\)wheneverρ​\(A\)<1\\rho\(A\)<1\. HereIIis the3×33\\times 3identity matrix,ρ​\(A\)\\rho\(A\)is the spectral radius ofAA, andb∗b^\{\\ast\}is the converged equilibrium state\.

### From three nodes toNNscalar nodes

The three\-node construction extends directly to an arbitrary number of scalar state nodes\. Let one equilibrium layer containNNnodes and let the encoded input bex∈ℝPx\\in\\mathbb\{R\}^\{P\}\. The stimulus at receiver nodeiiis

mi=∑r=1PWi​r\(s\)​xr\+ci\.m\_\{i\}=\\sum\_\{r=1\}^\{P\}W^\{\(s\)\}\_\{ir\}x\_\{r\}\+c\_\{i\}\.\(103\)Herei∈\{1,…,N\}i\\in\\\{1,\\ldots,N\\\}indexes state nodes,rrindexes input coordinates,W\(s\)∈ℝN×PW^\{\(s\)\}\\in\\mathbb\{R\}^\{N\\times P\}is the stimulus matrix, andc∈ℝNc\\in\\mathbb\{R\}^\{N\}is the stimulus bias\.

The full scalar interacting layer is

bi=mi\+∑j=1NAi​j​ϕ​\(bj\)\.b\_\{i\}=m\_\{i\}\+\\sum\_\{j=1\}^\{N\}A\_\{ij\}\\phi\(b\_\{j\}\)\.\(104\)Herejjindexes source nodes,Ai​jA\_\{ij\}is the receiver\-first interaction coefficient from sourcejjto receiverii, andϕ\\phiis the scalar nonlinearity\. In vector form,

b\\displaystyle b=m\+A​ϕ​\(b\),\\displaystyle=m\+A\\phi\(b\),\(105\)A\\displaystyle A=Alocal\+Aglobal\.\\displaystyle=A\_\{\\mathrm\{local\}\}\+A\_\{\\mathrm\{global\}\}\.\(106\)The matrixAlocalA\_\{\\mathrm\{local\}\}contains the sparse local interactions\. Its entry\(Alocal\)i​j\(A\_\{\\mathrm\{local\}\}\)\_\{ij\}is zero whenjjis not in the local neighborhood ofii\. The matrixAglobalA\_\{\\mathrm\{global\}\}contains dense or low\-rank context interactions\. For a scalar mean\-field broadcast, one example is

\(Aglobal\)i​j=βN,\(A\_\{\\mathrm\{global\}\}\)\_\{ij\}=\\frac\{\\beta\}\{N\},\(107\)whereβ\\betais the scalar global gate\. This is the scalar analogue of the mean\-field broadcast in Eq\. \([40](https://arxiv.org/html/2607.28989#S2.E40)\)\.

The damped solver for theNN\-node scalar system is

bk\+1=\(1−α\)​bk\+α​\[m\+\(Alocal\+Aglobal\)​ϕ​\(bk\)\]\.b\_\{k\+1\}=\(1\-\\alpha\)b\_\{k\}\+\\alpha\\big\[m\+\(A\_\{\\mathrm\{local\}\}\+A\_\{\\mathrm\{global\}\}\)\\phi\(b\_\{k\}\)\\big\]\.\(108\)Herebk∈ℝNb\_\{k\}\\in\\mathbb\{R\}^\{N\}is the solver state afterkkiterations andα\\alphais the damping coefficient\. Linearizing around a fixed point replacesϕ​\(bk\)\\phi\(b\_\{k\}\)by a local linear action, giving the finite\-solver transition

bk\+1\\displaystyle b\_\{k\+1\}=α​m\+Tα​bk,\\displaystyle=\\alpha m\+T\_\{\\alpha\}b\_\{k\},\(109\)Tα\\displaystyle T\_\{\\alpha\}=\(1−α\)​I\+α​A,\\displaystyle=\(1\-\\alpha\)I\+\\alpha A,\(110\)A\\displaystyle A=Alocal\+Aglobal\.\\displaystyle=A\_\{\\mathrm\{local\}\}\+A\_\{\\mathrm\{global\}\}\.\(111\)The identityIIhas sizeN×NN\\times N\. Starting fromb0=0b\_\{0\}=0, repeated substitution gives

b1\\displaystyle b\_\{1\}=α​m,\\displaystyle=\\alpha m,\(112\)b2\\displaystyle b\_\{2\}=α​m\+α​Tα​m,\\displaystyle=\\alpha m\+\\alpha T\_\{\\alpha\}m,\(113\)b3\\displaystyle b\_\{3\}=α​m\+α​Tα​m\+α​Tα2​m\.\\displaystyle=\\alpha m\+\\alpha T\_\{\\alpha\}m\+\\alpha T\_\{\\alpha\}^\{2\}m\.\(114\)Therefore, afterKKiterations,

bK=α​∑t=0K−1Tαt​m\.b\_\{K\}=\\alpha\\sum\_\{t=0\}^\{K\-1\}T\_\{\\alpha\}^\{t\}m\.\(115\)The termTαt​mT\_\{\\alpha\}^\{t\}mcontains every interaction history of lengthttin which the stimulus has been repeatedly carried by self\-persistence, local interaction, or global interaction\.

If the self\-persistence factors are collected separately, the same finite response can be expressed in powers of the learned interaction matrixAA,

bK\\displaystyle b\_\{K\}=∑r=0K−1wr\(K,α\)​Ar​m,\\displaystyle=\\sum\_\{r=0\}^\{K\-1\}w\_\{r\}^\{\(K,\\alpha\)\}A^\{r\}m,\(116\)wr\(K,α\)\\displaystyle w\_\{r\}^\{\(K,\\alpha\)\}=α​∑q=rK−1\(qr\)​\(1−α\)q−r​αr\.\\displaystyle=\\alpha\\sum\_\{q=r\}^\{K\-1\}\\binom\{q\}\{r\}\(1\-\\alpha\)^\{q\-r\}\\alpha^\{r\}\.\(117\)Hererris the number of local/global interaction applications andwr\(K,α\)w\_\{r\}^\{\(K,\\alpha\)\}is the weight assigned by a finite damped solver to all length\-rrinteraction histories\. SinceA=Alocal\+AglobalA=A\_\{\\mathrm\{local\}\}\+A\_\{\\mathrm\{global\}\},

Ar=∑w∈\{local,global\}rAwr​⋯​Aw1\.A^\{r\}=\\sum\_\{w\\in\\\{\\mathrm\{local\},\\mathrm\{global\}\\\}^\{r\}\}A\_\{w\_\{r\}\}\\cdots A\_\{w\_\{1\}\}\.\(118\)The wordw=\(w1,…,wr\)w=\(w\_\{1\},\\ldots,w\_\{r\}\)records which branch is applied at each interaction step\. This is the scalar multi\-node version of the path\-sum interpretation in Section[2\.10](https://arxiv.org/html/2607.28989#S2.SS10)\.

### From scalar nodes to vector\-feature nodes

The evaluated graph and molecule layers use vector features at each node\. LetB∈ℝN×dB\\in\\mathbb\{R\}^\{N\\times d\}be the state matrix, with rowBi∈ℝdB\_\{i\}\\in\\mathbb\{R\}^\{d\}holding theddfeatures of nodeii\. LetX∈ℝN×pX\\in\\mathbb\{R\}^\{N\\times p\}be the encoded node input\. A vector\-feature stimulus can be written

Si=Xi​Wstim⊤\+cstim\.S\_\{i\}=X\_\{i\}W\_\{\\mathrm\{stim\}\}^\{\\top\}\+c\_\{\\mathrm\{stim\}\}\.\(119\)HereSi∈ℝdS\_\{i\}\\in\\mathbb\{R\}^\{d\},Xi∈ℝpX\_\{i\}\\in\\mathbb\{R\}^\{p\},Wstim∈ℝd×pW\_\{\\mathrm\{stim\}\}\\in\\mathbb\{R\}^\{d\\times p\}, andcstim∈ℝdc\_\{\\mathrm\{stim\}\}\\in\\mathbb\{R\}^\{d\}\. LetY=χ𝒟​\(B\)Y=\\chi\_\{\\mathcal\{D\}\}\(B\)be the recurrent signal supplied to the interaction branches\. The local and global linearized branches can be written by blocks,

L​\(Y\)i\\displaystyle L\(Y\)\_\{i\}=∑j=1NAi​jL​Yj​WL⊤,\\displaystyle=\\sum\_\{j=1\}^\{N\}A^\{L\}\_\{ij\}Y\_\{j\}W\_\{L\}^\{\\top\},\(120\)G​\(Y\)i\\displaystyle G\(Y\)\_\{i\}=∑j=1NAi​jG​Yj​WG⊤\.\\displaystyle=\\sum\_\{j=1\}^\{N\}A^\{G\}\_\{ij\}Y\_\{j\}W\_\{G\}^\{\\top\}\.\(121\)HereALA^\{L\}is sparse over the local graph,AGA^\{G\}is dense or low\-rank over the global context, andWL,WG∈ℝd×dW\_\{L\},W\_\{G\}\\in\\mathbb\{R\}^\{d\\times d\}mix feature channels\. For ZINC,Y=BY=B\. For the evaluated vision and graph\-node models,Y=tanh⁡\(B\)Y=\\tanh\(B\)\.

Flatten the state by stacking all node\-feature vectors intoz=vec​\(B\)∈ℝN​dz=\\mathrm\{vec\}\(B\)\\in\\mathbb\{R\}^\{Nd\}and flatten the stimulus ass=vec​\(S\)s=\\mathrm\{vec\}\(S\)\. Before differentiating through the mapB↦YB\\mapsto Y, the branch matrices associated with Eqs\. \([120](https://arxiv.org/html/2607.28989#Sx2.E120)\)–\([121](https://arxiv.org/html/2607.28989#Sx2.E121)\) have receiver\-source blocks

M¯L,i​j\\displaystyle\\bar\{M\}\_\{L,ij\}=Ai​jL​WL,\\displaystyle=A^\{L\}\_\{ij\}W\_\{L\},\(122\)M¯G,i​j\\displaystyle\\bar\{M\}\_\{G,ij\}=Ai​jG​WG\.\\displaystyle=A^\{G\}\_\{ij\}W\_\{G\}\.\(123\)Each block has sized×dd\\times d\. Let

Dχ=∂vec​\(χ𝒟​\(B\)\)∂vec​\(B\)\.D\_\{\\chi\}=\\frac\{\\partial\\,\\mathrm\{vec\}\(\\chi\_\{\\mathcal\{D\}\}\(B\)\)\}\{\\partial\\,\\mathrm\{vec\}\(B\)\}\.\(124\)HereDχD\_\{\\chi\}is the identity for ZINC and the diagonal Jacobian oftanh\\tanhfor the evaluated vision and graph\-node models\. The full local linearizations with respect to the solver state are

ML\\displaystyle M\_\{L\}=M¯L​Dχ,\\displaystyle=\\bar\{M\}\_\{L\}D\_\{\\chi\},\(125\)MG\\displaystyle M\_\{G\}=M¯G​Dχ\.\\displaystyle=\\bar\{M\}\_\{G\}D\_\{\\chi\}\.\(126\)The matricesML,MG∈ℝN​d×N​dM\_\{L\},M\_\{G\}\\in\\mathbb\{R\}^\{Nd\\times Nd\}are therefore the local linearizations of the two interaction branches\. The linearized damped vector\-feature solver is

zk\+1≈α​s\+\[\(1−α\)​I\+α​ML\+α​MG\]​zk\.z\_\{k\+1\}\\approx\\alpha s\+\\big\[\(1\-\\alpha\)I\+\\alpha M\_\{L\}\+\\alpha M\_\{G\}\\big\]z\_\{k\}\.\(127\)This is the multi\-node, multi\-feature version of Eq\. \([10](https://arxiv.org/html/2607.28989#S2.E10)\)\. The same finite response follows,

zK=RK,α​s,z\_\{K\}=R\_\{K,\\alpha\}s,\(128\)with

RK,α=∑r=0K−1wr\(K,α\)​\(ML\+MG\)r\.R\_\{K,\\alpha\}=\\sum\_\{r=0\}^\{K\-1\}w\_\{r\}^\{\(K,\\alpha\)\}\(M\_\{L\}\+M\_\{G\}\)^\{r\}\.\(129\)HereRK,α∈ℝN​d×N​dR\_\{K,\\alpha\}\\in\\mathbb\{R\}^\{Nd\\times Nd\}is the finite response operator realized implicitly by the solver\. The matrix is not stored during the experiments\. Its action is represented by repeated calls to the learned update map\.

### Stacked and joint multi\-equilibrium form

The stacked SILVA architecture solves multiple equilibrium layers in sequence\. For layerℓ∈\{1,…,L\}\\ell\\in\\\{1,\\ldots,L\\\}, letNℓN\_\{\\ell\}be the number of state nodes,dℓd\_\{\\ell\}be the hidden width, andzℓ,k∈ℝNℓ​dℓz\_\{\\ell,k\}\\in\\mathbb\{R\}^\{N\_\{\\ell\}d\_\{\\ell\}\}be the flattened state at solver iterationkk\. The layer has its own stimulus mapsℓs\_\{\\ell\}, local operatorMℓ,LM\_\{\\ell,L\}, global operatorMℓ,GM\_\{\\ell,G\}, damping coefficientαℓ\\alpha\_\{\\ell\}, and solver budgetKℓK\_\{\\ell\}\. Its linearized recurrence is

zℓ,k\+1\\displaystyle z\_\{\\ell,k\+1\}=αℓ​sℓ​\(xℓ\)\+Tℓ,α​zℓ,k,\\displaystyle=\\alpha\_\{\\ell\}s\_\{\\ell\}\(x\_\{\\ell\}\)\+T\_\{\\ell,\\alpha\}z\_\{\\ell,k\},\(130\)Tℓ,α\\displaystyle T\_\{\\ell,\\alpha\}=\(1−αℓ\)​Iℓ\+αℓ​Mℓ,\\displaystyle=\(1\-\\alpha\_\{\\ell\}\)I\_\{\\ell\}\+\\alpha\_\{\\ell\}M\_\{\\ell\},\(131\)Mℓ\\displaystyle M\_\{\\ell\}=Mℓ,L\+Mℓ,G\.\\displaystyle=M\_\{\\ell,L\}\+M\_\{\\ell,G\}\.\(132\)HereIℓI\_\{\\ell\}is the identity on the flattened state space of layerℓ\\ell\. With zero initialization in each layer, the finite response of layerℓ\\ellis

zℓ,Kℓ\\displaystyle z\_\{\\ell,K\_\{\\ell\}\}=Rℓ,Kℓ,αℓ​sℓ​\(xℓ\),\\displaystyle=R\_\{\\ell,K\_\{\\ell\},\\alpha\_\{\\ell\}\}s\_\{\\ell\}\(x\_\{\\ell\}\),\(133\)Rℓ,Kℓ,αℓ\\displaystyle R\_\{\\ell,K\_\{\\ell\},\\alpha\_\{\\ell\}\}=αℓ​∑t=0Kℓ−1Tℓ,αt\.\\displaystyle=\\alpha\_\{\\ell\}\\sum\_\{t=0\}^\{K\_\{\\ell\}\-1\}T\_\{\\ell,\\alpha\}^\{t\}\.\(134\)The connector to the next equilibrium is

xℓ\+1=ψℓ​\(zℓ,Kℓ\)\.x\_\{\\ell\+1\}=\\psi\_\{\\ell\}\(z\_\{\\ell,K\_\{\\ell\}\}\)\.\(135\)Hereψℓ\\psi\_\{\\ell\}is the inter\-layer transformation\. In the evaluated stacks it includes the fixedtanh\\tanhconnector used between solved equilibria\.

For a local sensitivity derivation, linearize the stimulus and connector maps,

Sℓ\\displaystyle S\_\{\\ell\}=∂sℓ∂xℓ,\\displaystyle=\\frac\{\\partial s\_\{\\ell\}\}\{\\partial x\_\{\\ell\}\},\(136\)Cℓ\\displaystyle C\_\{\\ell\}=∂ψℓ∂zℓ,Kℓ\.\\displaystyle=\\frac\{\\partial\\psi\_\{\\ell\}\}\{\\partial z\_\{\\ell,K\_\{\\ell\}\}\}\.\(137\)HereSℓS\_\{\\ell\}maps perturbations in the layer input to perturbations in the layer stimulus, whileCℓC\_\{\\ell\}maps perturbations in one layer’s attractor to the next layer’s input\. Perturbations propagate through one solved equilibrium as

δ​zℓ,Kℓ=Rℓ,Kℓ,αℓ​Sℓ​δ​xℓ\.\\delta z\_\{\\ell,K\_\{\\ell\}\}=R\_\{\\ell,K\_\{\\ell\},\\alpha\_\{\\ell\}\}S\_\{\\ell\}\\delta x\_\{\\ell\}\.\(138\)For a two\-equilibrium stack,

δ​z2,K2=R2,K2,α2​S2​C1​R1,K1,α1​S1​δ​x1\.\\delta z\_\{2,K\_\{2\}\}=R\_\{2,K\_\{2\},\\alpha\_\{2\}\}S\_\{2\}C\_\{1\}R\_\{1,K\_\{1\},\\alpha\_\{1\}\}S\_\{1\}\\delta x\_\{1\}\.\(139\)For anLL\-equilibrium stack,

δ​zL,KL=RL,KL,αL​SL​CL−1​⋯​R2,K2,α2​S2​C1​R1,K1,α1​S1​δ​x1\.\\delta z\_\{L,K\_\{L\}\}=R\_\{L,K\_\{L\},\\alpha\_\{L\}\}S\_\{L\}C\_\{L\-1\}\\cdots R\_\{2,K\_\{2\},\\alpha\_\{2\}\}S\_\{2\}C\_\{1\}R\_\{1,K\_\{1\},\\alpha\_\{1\}\}S\_\{1\}\\delta x\_\{1\}\.\(140\)This expression is the multi\-node, multi\-layer version of the finite response operator\. EachRℓ,Kℓ,αℓR\_\{\\ell,K\_\{\\ell\},\\alpha\_\{\\ell\}\}sums paths inside one equilibrium\. Each connector JacobianCℓC\_\{\\ell\}transfers the settled vector attractor of one layer into the stimulus of the next\.

The finite unrolled backward pass through a sequential stack applies the same chain rule layer by layer\. Let𝒯ℓ,αℓ​\(zℓ,xℓ\)\\mathcal\{T\}\_\{\\ell,\\alpha\_\{\\ell\}\}\(z\_\{\\ell\},x\_\{\\ell\}\)be the damped update map of layerℓ\\ell\. The reverse recursion inside layerℓ\\ellis

z¯ℓ,k\\displaystyle\\bar\{z\}\_\{\\ell,k\}=J𝒯ℓ,αℓ​\(zℓ,k,xℓ\)⊤​z¯ℓ,k\+1,k=Kℓ−1,…,0,\\displaystyle=J\_\{\\mathcal\{T\}\_\{\\ell,\\alpha\_\{\\ell\}\}\}\(z\_\{\\ell,k\},x\_\{\\ell\}\)^\{\\top\}\\bar\{z\}\_\{\\ell,k\+1\},\\qquad k=K\_\{\\ell\}\-1,\\ldots,0,\(141\)∇θℓℒ\\displaystyle\\nabla\_\{\\theta\_\{\\ell\}\}\\mathcal\{L\}=∑k=0Kℓ−1\(∂θℓ𝒯ℓ,αℓ​\(zℓ,k,xℓ\)\)⊤​z¯ℓ,k\+1\.\\displaystyle=\\sum\_\{k=0\}^\{K\_\{\\ell\}\-1\}\\left\(\\partial\_\{\\theta\_\{\\ell\}\}\\mathcal\{T\}\_\{\\ell,\\alpha\_\{\\ell\}\}\(z\_\{\\ell,k\},x\_\{\\ell\}\)\\right\)^\{\\top\}\\bar\{z\}\_\{\\ell,k\+1\}\.\(142\)Herez¯ℓ,k\\bar\{z\}\_\{\\ell,k\}is the loss adjoint for thekkth iterate of layerℓ\\ell, andθℓ\\theta\_\{\\ell\}denotes only the parameters of that layer\. The terminal adjointz¯ℓ,Kℓ\\bar\{z\}\_\{\\ell,K\_\{\\ell\}\}contains the readout gradient for the last layer and the connector contributionCℓ⊤​x¯ℓ\+1C\_\{\\ell\}^\{\\top\}\\bar\{x\}\_\{\\ell\+1\}for earlier layers\.

If each layer is differentiated as an exact implicit equilibrium instead, the stagewise adjoint solves

\(I−Jℓ⊤\)​λℓ\\displaystyle\\left\(I\-J\_\{\\ell\}^\{\\top\}\\right\)\\lambda\_\{\\ell\}=z¯ℓ∗,\\displaystyle=\\bar\{z\}\_\{\\ell\}^\{\\ast\},\(143\)∇θℓℒ\\displaystyle\\nabla\_\{\\theta\_\{\\ell\}\}\\mathcal\{L\}=\(∂θℓfθℓ\(ℓ\)​\(zℓ∗,xℓ\)\)⊤​λℓ,\\displaystyle=\\left\(\\partial\_\{\\theta\_\{\\ell\}\}f\_\{\\theta\_\{\\ell\}\}^\{\(\\ell\)\}\(z\_\{\\ell\}^\{\\ast\},x\_\{\\ell\}\)\\right\)^\{\\top\}\\lambda\_\{\\ell\},\(144\)x¯ℓ\\displaystyle\\bar\{x\}\_\{\\ell\}=Hℓ⊤​λℓ\.\\displaystyle=H\_\{\\ell\}^\{\\top\}\\lambda\_\{\\ell\}\.\(145\)HereJℓ=∂fθℓ\(ℓ\)/∂zℓJ\_\{\\ell\}=\\partial f\_\{\\theta\_\{\\ell\}\}^\{\(\\ell\)\}/\\partial z\_\{\\ell\}andHℓ=∂fθℓ\(ℓ\)/∂xℓH\_\{\\ell\}=\\partial f\_\{\\theta\_\{\\ell\}\}^\{\(\\ell\)\}/\\partial x\_\{\\ell\}are evaluated at the layer fixed point\. The upstream gradient for the previous attractor isCℓ−1⊤​x¯ℓC\_\{\\ell\-1\}^\{\\top\}\\bar\{x\}\_\{\\ell\}\. These equations describe the implicit analogue of the sequential stack\. The experiments with stacked SILVA layers use the finite unrolled rule above unless explicitly identified as adjoint\-trained\.

The same sequential stack can be written as one block fixed\-point system for analysis\. Define the joint state

𝐳=\(z1,z2,…,zL\)\.\\mathbf\{z\}=\(z\_\{1\},z\_\{2\},\\ldots,z\_\{L\}\)\.\(146\)The feed\-forward\-coupled fixed\-point equations are

z1\\displaystyle z\_\{1\}=fθ1\(1\)​\(z1,x1\),\\displaystyle=f\_\{\\theta\_\{1\}\}^\{\(1\)\}\(z\_\{1\},x\_\{1\}\),\(147\)z2\\displaystyle z\_\{2\}=fθ2\(2\)​\(z2,ψ1​\(z1\)\),\\displaystyle=f\_\{\\theta\_\{2\}\}^\{\(2\)\}\(z\_\{2\},\\psi\_\{1\}\(z\_\{1\}\)\),\(148\)zL\\displaystyle z\_\{L\}=fθL\(L\)​\(zL,ψL−1​\(zL−1\)\)\.\\displaystyle=f\_\{\\theta\_\{L\}\}^\{\(L\)\}\(z\_\{L\},\\psi\_\{L\-1\}\(z\_\{L\-1\}\)\)\.\(149\)The omitted middle equations follow the same pattern forℓ=3,…,L−1\\ell=3,\\ldots,L\-1\. Linearizing this joint system gives a block lower bidiagonal Jacobian,

J𝐅=\[J100⋯0H2​C1J20⋯00H3​C2J3⋯0⋮⋮⋮⋱⋮000HL​CL−1JL\]\.J\_\{\\mathbf\{F\}\}=\\begin\{bmatrix\}J\_\{1\}&0&0&\\cdots&0\\\\ H\_\{2\}C\_\{1\}&J\_\{2\}&0&\\cdots&0\\\\ 0&H\_\{3\}C\_\{2\}&J\_\{3\}&\\cdots&0\\\\ \\vdots&\\vdots&\\vdots&\\ddots&\\vdots\\\\ 0&0&0&H\_\{L\}C\_\{L\-1\}&J\_\{L\}\\end\{bmatrix\}\.\(150\)HereJℓ=∂fθℓ\(ℓ\)/∂zℓJ\_\{\\ell\}=\\partial f\_\{\\theta\_\{\\ell\}\}^\{\(\\ell\)\}/\\partial z\_\{\\ell\}is the within\-equilibrium Jacobian,Hℓ=∂fθℓ\(ℓ\)/∂xℓH\_\{\\ell\}=\\partial f\_\{\\theta\_\{\\ell\}\}^\{\(\\ell\)\}/\\partial x\_\{\\ell\}is the input\-to\-update Jacobian, andCℓC\_\{\\ell\}is the connector Jacobian from Eq\. \([137](https://arxiv.org/html/2607.28989#Sx2.E137)\)\. The off\-diagonal blockHℓ\+1​CℓH\_\{\\ell\+1\}C\_\{\\ell\}is nonzero because layerℓ\+1\\ell\+1receives its input from the attractor of layerℓ\\ell\.

The joint implicit response is

δ​𝐳∗=\(I−J𝐅\)−1​𝐒​δ​x1,\\delta\\mathbf\{z\}^\{\\ast\}=\(I\-J\_\{\\mathbf\{F\}\}\)^\{\-1\}\\mathbf\{S\}\\,\\delta x\_\{1\},\(151\)where𝐒\\mathbf\{S\}injects the external input perturbation into the first layer’s stimulus block\. For the sequential SILVA stack, the lower\-bidiagonal structure of Eq\. \([150](https://arxiv.org/html/2607.28989#Sx2.E150)\) makes this inverse equivalent to the stagewise response in Eq\. \([140](https://arxiv.org/html/2607.28989#Sx2.E140)\)\. This is different from a simultaneous multiscale DEQ, where the off\-diagonal coupling blocks are solved together as one enlarged equilibrium rather than as separately settled vector attractors\.

### Adjoint path for the MNIST equilibrium diagnostic

The experiments use two backward paths\. Finite\-solver experiments differentiate the executed Picard iterations as in Eqs\. \([71](https://arxiv.org/html/2607.28989#S2.E71)\)–\([73](https://arxiv.org/html/2607.28989#S2.E73)\)\. The separate adjoint\-trained MNIST diagnostic uses a linear adjoint solve at the terminal equilibrium state of the adjoint\-trained model\. For one such layer, letm=sθ​\(x\)m=s\_\{\\theta\}\(x\)be the stimulus and let the pre\-activation fixed point be

b∗=m\+Lθ​\(ϕ​\(b∗\)\)\+Gθ​\(ϕ​\(b∗\)\),b^\{\\ast\}=m\+L\_\{\\theta\}\(\\phi\(b^\{\\ast\}\)\)\+G\_\{\\theta\}\(\\phi\(b^\{\\ast\}\)\),\(152\)withϕ​\(b\)=tanh⁡\(b\)\\phi\(b\)=\\tanh\(b\)\. The forward Picard loop used for this diagnostic is undamped\. No factorα\\alphaappears in the layer’s forward recurrence or in its GMRES backward operator\. The layer output isB∗=ϕ​\(b∗\)B^\{\\ast\}=\\phi\(b^\{\\ast\}\)\. IfB¯=∇B∗ℒ\\bar\{B\}=\\nabla\_\{B^\{\\ast\}\}\\mathcal\{L\}, then the adjoint entering the pre\-activation equilibrium state is

b¯=B¯⊙\(1−tanh2⁡\(b∗\)\)\.\\bar\{b\}=\\bar\{B\}\\odot\\big\(1\-\\tanh^\{2\}\(b^\{\\ast\}\)\\big\)\.\(153\)Here⊙\\odotis elementwise multiplication, andb¯\\bar\{b\}has the same shape asb∗b^\{\\ast\}\.

The GMRES solve is applied to the local and global linearized actions,

𝒜​\(v\)\\displaystyle\\mathcal\{A\}\(v\)=v−𝒥L​\(v\)−𝒥G​\(v\),\\displaystyle=v\-\\mathcal\{J\}\_\{L\}\(v\)\-\\mathcal\{J\}\_\{G\}\(v\),\(154\)𝒜​\(u\)\\displaystyle\\mathcal\{A\}\(u\)=b¯\.\\displaystyle=\\bar\{b\}\.\(155\)Herevvanduuare flattened state\-shaped adjoint vectors\. The maps𝒥L\\mathcal\{J\}\_\{L\}and𝒥G\\mathcal\{J\}\_\{G\}are the local\- and global\-interaction linear actions evaluated at the terminal equilibrium state\. They are not materialized as matrices\. For the static global matrix in the row\-oriented convention, the global action has coordinates

𝒥G​\(v\)i=∑j=1NAj​i​vj​ϕ′​\(bj∗\)\.\\mathcal\{J\}\_\{G\}\(v\)\_\{i\}=\\sum\_\{j=1\}^\{N\}A\_\{ji\}\\,v\_\{j\}\\,\\phi^\{\\prime\}\(b^\{\\ast\}\_\{j\}\)\.\(156\)HereAAis the learned global interaction matrix,iiindexes the output coordinate of the linear action, andjjindexes the source coordinate in the summation\. For the local attention branch,𝒥L​\(v\)\\mathcal\{J\}\_\{L\}\(v\)is obtained by distributing the incoming per\-node adjoint evenly over theHHattention heads, as in Eq\. \([88](https://arxiv.org/html/2607.28989#S2.E88)\), and then applying the local attention layer’s backward map\.

The solutionuuis the adjoint used to accumulate the diagnostic gradients\. For a static global matrix, the per\-example matrix sensitivity is

∇Aℒ^i​j=ϕ​\(bi∗\)​uj\.\\widehat\{\\nabla\_\{A\}\\mathcal\{L\}\}\_\{ij\}=\\phi\(b^\{\\ast\}\_\{i\}\)u\_\{j\}\.\(157\)For a stimulus matrixWWand biascc, the direct gradients are

∇Wℒ^\\displaystyle\\widehat\{\\nabla\_\{W\}\\mathcal\{L\}\}=X⊤​\(u\+dmattn\),\\displaystyle=X^\{\\top\}\(u\+d\_\{m\}^\{\\mathrm\{attn\}\}\),\(158\)∇cℒ^\\displaystyle\\widehat\{\\nabla\_\{c\}\\mathcal\{L\}\}=∑n\(un\+dm,nattn\)\.\\displaystyle=\\sum\_\{n\}\(u\_\{n\}\+d\_\{m,n\}^\{\\mathrm\{attn\}\}\)\.\(159\)HereXXis the minibatch input to the layer,nnindexes minibatch elements, anddmattnd\_\{m\}^\{\\mathrm\{attn\}\}is the additional stimulus adjoint obtained when the global interaction matrix is generated by attention frommm\. The hat on these gradients distinguishes this diagnostic adjoint path from the formal fixed\-point gradient in Eqs\. \([77](https://arxiv.org/html/2607.28989#S2.E77)\)–\([79](https://arxiv.org/html/2607.28989#S2.E79)\)\. The conditioning of such forward and backward linearized solves is controlled by the equilibrium Jacobian, motivating both the spectral\-radius diagnostics used here and Jacobian\-regularized DEQ training in related work\[[8](https://arxiv.org/html/2607.28989#bib.bib7)\]\.

### Readout after the equilibrium state

For binary classification in the three\-node scalar guide, the readout can be written

y^\\displaystyle\\hat\{y\}=σsigmoid​\(v⊤​ϕ​\(b∗\)\+d\),\\displaystyle=\\sigma\_\{\\mathrm\{sigmoid\}\}\(v^\{\\top\}\\phi\(b^\{\\ast\}\)\+d\),\(160\)v\\displaystyle v∈ℝ3,\\displaystyle\\in\\mathbb\{R\}^\{3\},\(161\)d\\displaystyle d∈ℝ,\\displaystyle\\in\\mathbb\{R\},\(162\)Herey^\\hat\{y\}is the predicted probability for the positive class,vvis the readout vector,ddis the scalar readout bias, andσsigmoid\\sigma\_\{\\mathrm\{sigmoid\}\}is the logistic sigmoid\. The model is trained with binary cross\-entropy\. The multi\-class experiments in the main text use the same principle with a linear map followed by softmax, applied to the approximate fixed point returned by the finite solver\.

## Acknowledgements

The present work introduces the SILVA Networks methodology and its structured implicit interaction field\. The research program from which this study emerged was initiated within the Swedish National Infrastructure for Computing \(SNIC\) Small Compute project*Artificial Intelligence for Physics and Engineering, Modeling and Simulation*, Project No\. SNIC 2022/22\-843\. The project was conducted at Linköping University under the author’s principal investigatorship\. That project established a broader research direction connecting artificial intelligence, graph neural networks, deep reinforcement learning, neuromorphic computing, and the modeling and simulation of physics and engineering systems\.

This work was supported by the Brazilian National Council for Scientific and Technological Development \(CNPq\) under grant No\. 445344/2024\-5\. The author also acknowledges financial support provided through the Program Talentos Brasil, as Project Investigator\.

## Declarations

- •Funding\.This work was supported by the Brazilian National Council for Scientific and Technological Development \(CNPq\) under grant No\. 445344/2024\-5\. The author acknowledges financial support provided through the Program Talentos Brasil, as Project Investigator\.
- •Competing interests\.The author declares no competing interests\.
- •Ethics approval\.Not applicable\.
- •Data and code availability\.MNIST\[[37](https://arxiv.org/html/2607.28989#bib.bib52)\], CIFAR\-10\[[36](https://arxiv.org/html/2607.28989#bib.bib57)\], ZINC\[[29](https://arxiv.org/html/2607.28989#bib.bib21)\], Cora, Citeseer, Pubmed\[[49](https://arxiv.org/html/2607.28989#bib.bib41),[58](https://arxiv.org/html/2607.28989#bib.bib13)\], and CLUSTER\[[20](https://arxiv.org/html/2607.28989#bib.bib17)\]are public benchmarks\. Exact split and preprocessing details are given in Section[2](https://arxiv.org/html/2607.28989#S2)and the Methods summary\. The computational implementation used Python with standard scientific\-computing and deep\-learning software\. SILVA Networks research materials, processed experiment outputs, and additional supporting materials will be available from the corresponding author upon reasonable request\.
- •Author contributions\.J\.L\.L\.J\.S\. conceived the research direction, developed the SILVA Networks model and benchmark suite, designed and executed the experiments, interpreted the results, and prepared the manuscript\.

## References

- \[1\]\(2021\)On the bottleneck of graph neural networks and its practical implications\.InInternational Conference on Learning Representations \(ICLR\),External Links:[Link](https://openreview.net/forum?id=i80OPhOCVH2)Cited by:[§1](https://arxiv.org/html/2607.28989#S1.p4.1)\.
- \[2\]S\. Amari\(1977\)Dynamics of pattern formation in lateral\-inhibition type neural fields\.Biological Cybernetics27\(2\),pp\. 77–87\.External Links:[Document](https://dx.doi.org/10.1007/BF00337259)Cited by:[§1](https://arxiv.org/html/2607.28989#S1.p5.1),[§2\.1](https://arxiv.org/html/2607.28989#S2.SS1.SSS0.Px2.p1.1)\.
- \[3\]D\. G\. Anderson\(1965\)Iterative procedures for nonlinear integral equations\.Journal of the ACM12\(4\),pp\. 547–560\.External Links:[Document](https://dx.doi.org/10.1145/321296.321305)Cited by:[§2\.1](https://arxiv.org/html/2607.28989#S2.SS1.SSS0.Px3.p1.3),[§2\.9](https://arxiv.org/html/2607.28989#S2.SS9.p1.3)\.
- \[4\]L\. J\. Ba, J\. R\. Kiros, and G\. E\. Hinton\(2016\)Layer normalization\.arXiv preprint arXiv:1607\.06450\.External Links:[Document](https://dx.doi.org/10.48550/arXiv.1607.06450),[Link](https://arxiv.org/abs/1607.06450)Cited by:[§2\.1](https://arxiv.org/html/2607.28989#S2.SS1.p2.7)\.
- \[5\]S\. Bae, H\. Byun, C\. Oh, Y\. Cho, and K\. Song\(2022\)Graph perceiver IO: a general architecture for graph structured data\.arXiv preprint arXiv:2209\.06418\.External Links:[Document](https://dx.doi.org/10.48550/arXiv.2209.06418),[Link](https://arxiv.org/abs/2209.06418)Cited by:[§1](https://arxiv.org/html/2607.28989#S1.p4.1),[§3](https://arxiv.org/html/2607.28989#S3.SS0.SSS0.Px8.p1.1),[§3](https://arxiv.org/html/2607.28989#S3.SS0.SSS0.Px8.p2.3),[§4](https://arxiv.org/html/2607.28989#S4.p4.1)\.
- \[6\]S\. Bai, J\. Z\. Kolter, and V\. Koltun\(2019\)Deep equilibrium models\.InAdvances in Neural Information Processing Systems \(NeurIPS\),Vol\.32\.External Links:[Link](https://papers.nips.cc/paper/8358-deep-equilibrium-models)Cited by:[§1](https://arxiv.org/html/2607.28989#S1.p1.1),[§1](https://arxiv.org/html/2607.28989#S1.p3.1),[§2\.13](https://arxiv.org/html/2607.28989#S2.SS13.SSS0.Px2.p1.26),[§2\.13](https://arxiv.org/html/2607.28989#S2.SS13.SSS0.Px2.p1.35),[§2\.9](https://arxiv.org/html/2607.28989#S2.SS9.p1.3),[§4](https://arxiv.org/html/2607.28989#S4.p4.1)\.
- \[7\]S\. Bai, V\. Koltun, and J\. Z\. Kolter\(2020\)Multiscale deep equilibrium models\.InAdvances in Neural Information Processing Systems \(NeurIPS\),Vol\.33\.External Links:[Link](https://papers.nips.cc/paper/2020/hash/3812f9a59b634c2a9c574610eaba5bed-Abstract.html)Cited by:[§1](https://arxiv.org/html/2607.28989#S1.p3.1),[§2\.1](https://arxiv.org/html/2607.28989#S2.SS1.SSS0.Px1.p2.38),[§4](https://arxiv.org/html/2607.28989#S4.p4.1)\.
- \[8\]S\. Bai, V\. Koltun, and J\. Z\. Kolter\(2021\)Stabilizing equilibrium models by jacobian regularization\.arXiv preprint arXiv:2106\.14342\.External Links:[Document](https://dx.doi.org/10.48550/arXiv.2106.14342),[Link](https://arxiv.org/abs/2106.14342)Cited by:[§1](https://arxiv.org/html/2607.28989#S1.p3.1),[§2\.7](https://arxiv.org/html/2607.28989#S2.SS7.p1.4),[Adjoint path for the MNIST equilibrium diagnostic](https://arxiv.org/html/2607.28989#Sx2.SSx5.p3.7)\.
- \[9\]J\. Baker, Q\. Wang, C\. D\. Hauck, and B\. Wang\(2023\)Implicit graph neural networks: a monotone operator viewpoint\.InProceedings of the 40th International Conference on Machine Learning,Proceedings of Machine Learning Research, Vol\.202,pp\. 1521–1548\.External Links:[Link](https://proceedings.mlr.press/v202/baker23a.html)Cited by:[§1](https://arxiv.org/html/2607.28989#S1.p3.1)\.
- \[10\]S\. Banach\(1922\)Sur les opérations dans les ensembles abstraits et leur application aux équations intégrales\.Fundamenta Mathematicae3,pp\. 133–181\.External Links:[Document](https://dx.doi.org/10.4064/fm-3-1-133-181),[Link](https://www.impan.pl/en/publishing-house/journals-and-series/fundamenta-mathematicae/all/3/0/92453/sur-les-operations-dans-les-ensembles-abstraits-et-leur-application-aux-equations-integrales)Cited by:[§2\.1](https://arxiv.org/html/2607.28989#S2.SS1.SSS0.Px1.p2.1),[§2\.1](https://arxiv.org/html/2607.28989#S2.SS1.p1.7)\.
- \[11\]P\. W\. Battaglia, R\. Pascanu, M\. Lai, D\. J\. Rezende, and K\. Kavukcuoglu\(2016\)Interaction networks for learning about objects, relations and physics\.InAdvances in Neural Information Processing Systems \(NeurIPS\),Vol\.29,pp\. 4502–4510\.External Links:[Link](https://papers.nips.cc/paper/6418-interaction-networks-for-learning-about-objects-relations-and-physics)Cited by:[§1](https://arxiv.org/html/2607.28989#S1.p1.1),[§1](https://arxiv.org/html/2607.28989#S1.p4.1),[§2\.1](https://arxiv.org/html/2607.28989#S2.SS1.SSS0.Px2.p1.1)\.
- \[12\]R\. Ben\-Yishai, R\. L\. Bar\-Or, and H\. Sompolinsky\(1995\)Theory of orientation tuning in visual cortex\.Proceedings of the National Academy of Sciences92\(9\),pp\. 3844–3848\.External Links:[Document](https://dx.doi.org/10.1073/pnas.92.9.3844)Cited by:[§1](https://arxiv.org/html/2607.28989#S1.p5.1)\.
- \[13\]N\. Bertschinger and T\. Natschläger\(2004\)Real\-time computation at the edge of chaos in recurrent neural networks\.Neural Computation16\(7\),pp\. 1413–1436\.External Links:[Document](https://dx.doi.org/10.1162/089976604323057443)Cited by:[§2\.10](https://arxiv.org/html/2607.28989#S2.SS10.SSS0.Px1.p4.4)\.
- \[14\]X\. Bresson and T\. Laurent\(2017\)Residual gated graph convnets\.arXiv preprint arXiv:1711\.07553\.External Links:[Document](https://dx.doi.org/10.48550/arXiv.1711.07553),[Link](https://arxiv.org/abs/1711.07553)Cited by:[§3](https://arxiv.org/html/2607.28989#S3.SS0.SSS0.Px4.p1.2)\.
- \[15\]M\. M\. Bronstein, J\. Bruna, Y\. LeCun, A\. Szlam, and P\. Vandergheynst\(2017\)Geometric deep learning: going beyond Euclidean data\.IEEE Signal Processing Magazine34\(4\),pp\. 18–42\.External Links:[Document](https://dx.doi.org/10.1109/MSP.2017.2693418)Cited by:[§1](https://arxiv.org/html/2607.28989#S1.p4.1)\.
- \[16\]R\. T\. Q\. Chen, Y\. Rubanova, J\. Bettencourt, and D\. K\. Duvenaud\(2018\)Neural ordinary differential equations\.InAdvances in Neural Information Processing Systems \(NeurIPS\),Vol\.31\.External Links:[Link](https://papers.nips.cc/paper/7892-neural-ordinary-differential-equations)Cited by:[§1](https://arxiv.org/html/2607.28989#S1.p2.1)\.
- \[17\]H\. Chu, S\. Wei, T\. Liu, Y\. Zhao, and Y\. Miyatake\(2024\)Lyapunov\-stable deep equilibrium models\.InProceedings of the AAAI Conference on Artificial Intelligence,Vol\.38,pp\. 11615–11623\.External Links:[Document](https://dx.doi.org/10.1609/aaai.v38i10.29044),[Link](https://ojs.aaai.org/index.php/AAAI/article/view/29044)Cited by:[§2\.8](https://arxiv.org/html/2607.28989#S2.SS8.p2.2),[§3](https://arxiv.org/html/2607.28989#S3.SS0.SSS0.Px13.p2.11)\.
- \[18\]G\. Corso, L\. Cavalleri, D\. Beaini, P\. Liò, and P\. Veličković\(2020\)Principal neighbourhood aggregation for graph nets\.InAdvances in Neural Information Processing Systems \(NeurIPS\),Vol\.33\.External Links:[Link](https://papers.nips.cc/paper/2020/hash/99cad265a1768cc2dd013f0e740300ae-Abstract.html)Cited by:[§3](https://arxiv.org/html/2607.28989#S3.SS0.SSS0.Px4.p1.2),[Table 1](https://arxiv.org/html/2607.28989#S3.T1.10.6.2)\.
- \[19\]D\. L\. Davies and D\. W\. Bouldin\(1979\)A cluster separation measure\.IEEE Transactions on Pattern Analysis and Machine IntelligencePAMI\-1\(2\),pp\. 224–227\.External Links:[Document](https://dx.doi.org/10.1109/TPAMI.1979.4766909)Cited by:[§3](https://arxiv.org/html/2607.28989#S3.SS0.SSS0.Px2.p1.4)\.
- \[20\]V\. P\. Dwivedi, C\. K\. Joshi, T\. Laurent, Y\. Bengio, X\. Bresson, and A\. T\. Luu\(2023\)Benchmarking graph neural networks\.Journal of Machine Learning Research24\(43\),pp\. 1–48\.External Links:[Link](https://www.jmlr.org/papers/v24/22-0567.html)Cited by:[§1](https://arxiv.org/html/2607.28989#S1.p6.1),[§2\.1](https://arxiv.org/html/2607.28989#S2.SS1.p1.18),[§2\.1](https://arxiv.org/html/2607.28989#S2.SS1.p2.27),[§3](https://arxiv.org/html/2607.28989#S3.SS0.SSS0.Px4.p1.2),[§3](https://arxiv.org/html/2607.28989#S3.SS0.SSS0.Px5.p2.6),[§3](https://arxiv.org/html/2607.28989#S3.SS0.SSS0.Px6.p1.16),[Table 1](https://arxiv.org/html/2607.28989#S3.T1.11.7.2),[§4](https://arxiv.org/html/2607.28989#S4.p3.4),[4th item](https://arxiv.org/html/2607.28989#Sx4.I1.i4.p1.1)\.
- \[21\]S\. El Hihi and Y\. Bengio\(1995\)Hierarchical recurrent neural networks for long\-term dependencies\.InAdvances in Neural Information Processing Systems \(NeurIPS\),Vol\.8\.External Links:[Link](https://papers.nips.cc/paper/1102-hierarchical-recurrent-neural-networks-for-long-term-dependencies)Cited by:[§2\.1](https://arxiv.org/html/2607.28989#S2.SS1.SSS0.Px1.p2.38),[§3](https://arxiv.org/html/2607.28989#S3.SS0.SSS0.Px7.p1.1)\.
- \[22\]C\. D\. Gilbert and W\. Li\(2013\)Top\-down influences on visual processing\.Nature Reviews Neuroscience14,pp\. 350–363\.External Links:[Document](https://dx.doi.org/10.1038/nrn3476)Cited by:[§1](https://arxiv.org/html/2607.28989#S1.p5.1),[§2\.1](https://arxiv.org/html/2607.28989#S2.SS1.SSS0.Px2.p1.1)\.
- \[23\]J\. Gilmer, S\. S\. Schoenholz, P\. F\. Riley, O\. Vinyals, and G\. E\. Dahl\(2017\)Neural message passing for quantum chemistry\.InProceedings of the 34th International Conference on Machine Learning,Proceedings of Machine Learning Research, Vol\.70,pp\. 1263–1272\.External Links:[Link](https://proceedings.mlr.press/v70/gilmer17a.html)Cited by:[§1](https://arxiv.org/html/2607.28989#S1.p1.1),[§1](https://arxiv.org/html/2607.28989#S1.p4.1),[§2\.1](https://arxiv.org/html/2607.28989#S2.SS1.SSS0.Px2.p1.1),[§2\.3](https://arxiv.org/html/2607.28989#S2.SS3.p1.12)\.
- \[24\]F\. Gu, H\. Chang, W\. Zhu, S\. Sojoudi, and L\. El Ghaoui\(2020\)Implicit graph neural networks\.InAdvances in Neural Information Processing Systems \(NeurIPS\),Vol\.33,pp\. 11984–11995\.External Links:[Link](https://papers.nips.cc/paper/2020/hash/8b5c8441a8ff8e151b191c53c1842a38-Abstract.html)Cited by:[§1](https://arxiv.org/html/2607.28989#S1.p1.1),[§1](https://arxiv.org/html/2607.28989#S1.p3.1)\.
- \[25\]E\. Haber and L\. Ruthotto\(2017\)Stable architectures for deep neural networks\.Inverse Problems34\(1\)\.External Links:[Document](https://dx.doi.org/10.1088/1361-6420/aa9a90)Cited by:[§2\.14](https://arxiv.org/html/2607.28989#S2.SS14.p1.8)\.
- \[26\]K\. He, X\. Zhang, S\. Ren, and J\. Sun\(2016\)Deep residual learning for image recognition\.InProceedings of the IEEE Conference on Computer Vision and Pattern Recognition \(CVPR\),pp\. 770–778\.External Links:[Document](https://dx.doi.org/10.1109/CVPR.2016.90),[Link](https://openaccess.thecvf.com/content_cvpr_2016/html/He_Deep_Residual_Learning_CVPR_2016_paper.html)Cited by:[§1](https://arxiv.org/html/2607.28989#S1.p2.1)\.
- \[27\]J\. J\. Hopfield\(1982\)Neural networks and physical systems with emergent collective computational abilities\.Proceedings of the National Academy of Sciences79\(8\),pp\. 2554–2558\.External Links:[Document](https://dx.doi.org/10.1073/pnas.79.8.2554)Cited by:[§1](https://arxiv.org/html/2607.28989#S1.p2.1),[§1](https://arxiv.org/html/2607.28989#S1.p5.1),[§2\.1](https://arxiv.org/html/2607.28989#S2.SS1.SSS0.Px2.p1.1),[§2\.8](https://arxiv.org/html/2607.28989#S2.SS8.p1.1),[§4](https://arxiv.org/html/2607.28989#S4.p4.1)\.
- \[28\]S\. Ioffe and C\. Szegedy\(2015\)Batch normalization: accelerating deep network training by reducing internal covariate shift\.InProceedings of the 32nd International Conference on Machine Learning,Proceedings of Machine Learning Research, Vol\.37,pp\. 448–456\.External Links:[Link](https://proceedings.mlr.press/v37/ioffe15.html)Cited by:[§2\.4](https://arxiv.org/html/2607.28989#S2.SS4.p1.6)\.
- \[29\]J\. J\. Irwin, T\. Sterling, M\. M\. Mysinger, E\. S\. Bolstad, and R\. G\. Coleman\(2012\)ZINC: a free tool to discover chemistry for biology\.Journal of Chemical Information and Modeling52\(7\),pp\. 1757–1768\.External Links:[Document](https://dx.doi.org/10.1021/ci3001277)Cited by:[§1](https://arxiv.org/html/2607.28989#S1.p6.1),[§2\.1](https://arxiv.org/html/2607.28989#S2.SS1.p2.27),[§3](https://arxiv.org/html/2607.28989#S3.SS0.SSS0.Px4.p1.2),[4th item](https://arxiv.org/html/2607.28989#Sx4.I1.i4.p1.1)\.
- \[30\]A\. Jaegle, S\. Borgeaud, J\. Alayrac, C\. Doersch, C\. Ionescu, D\. Ding, S\. Koppula, D\. Zoran, A\. Brock, E\. Shelhamer, O\. Hénaff, M\. M\. Botvinick, A\. Zisserman, O\. Vinyals, and J\. Carreira\(2022\)Perceiver IO: a general architecture for structured inputs & outputs\.InInternational Conference on Learning Representations \(ICLR\),External Links:[Link](https://openreview.net/forum?id=fILj7WpI-g)Cited by:[§1](https://arxiv.org/html/2607.28989#S1.p4.1),[§3](https://arxiv.org/html/2607.28989#S3.SS0.SSS0.Px8.p1.1)\.
- \[31\]I\. T\. Jolliffe\(2002\)Principal component analysis\.2 edition,Springer Series in Statistics,Springer\.External Links:[Document](https://dx.doi.org/10.1007/b98835)Cited by:[§3](https://arxiv.org/html/2607.28989#S3.SS0.SSS0.Px2.p1.4)\.
- \[32\]M\. Khona and I\. R\. Fiete\(2022\)Attractor and integrator networks in the brain\.Nature Reviews Neuroscience23,pp\. 744–766\.External Links:[Document](https://dx.doi.org/10.1038/s41583-022-00642-0)Cited by:[§1](https://arxiv.org/html/2607.28989#S1.p5.1)\.
- \[33\]D\. P\. Kingma and J\. Ba\(2015\)Adam: a method for stochastic optimization\.InInternational Conference on Learning Representations \(ICLR\),External Links:[Link](https://arxiv.org/abs/1412.6980)Cited by:[§2\.12](https://arxiv.org/html/2607.28989#S2.SS12.p1.10)\.
- \[34\]T\. N\. Kipf and M\. Welling\(2017\)Semi\-supervised classification with graph convolutional networks\.InInternational Conference on Learning Representations \(ICLR\),External Links:[Link](https://openreview.net/forum?id=SJU4ayYgl)Cited by:[§1](https://arxiv.org/html/2607.28989#S1.p4.1),[§2\.3](https://arxiv.org/html/2607.28989#S2.SS3.p1.12),[§3](https://arxiv.org/html/2607.28989#S3.SS0.SSS0.Px4.p1.2),[§3](https://arxiv.org/html/2607.28989#S3.SS0.SSS0.Px5.p1.1),[Table 1](https://arxiv.org/html/2607.28989#S3.T1.7.3.2),[Table 2](https://arxiv.org/html/2607.28989#S3.T2.14.11.1.2),[§4](https://arxiv.org/html/2607.28989#S4.p4.1)\.
- \[35\]J\. Koutník, K\. Greff, F\. Gomez, and J\. Schmidhuber\(2014\)A clockwork rnn\.InProceedings of the 31st International Conference on Machine Learning,Proceedings of Machine Learning Research, Vol\.32,pp\. 1863–1871\.External Links:[Link](https://proceedings.mlr.press/v32/koutnik14.html)Cited by:[§2\.1](https://arxiv.org/html/2607.28989#S2.SS1.SSS0.Px1.p2.38),[§3](https://arxiv.org/html/2607.28989#S3.SS0.SSS0.Px7.p1.1)\.
- \[36\]A\. Krizhevsky and G\. Hinton\(2009\)Learning multiple layers of features from tiny images\.Technical reportUniversity of Toronto\.External Links:[Link](https://www.cs.toronto.edu/%CB%9Ckriz/learning-features-2009-TR.pdf)Cited by:[§2\.1](https://arxiv.org/html/2607.28989#S2.SS1.p1.18),[§2\.1](https://arxiv.org/html/2607.28989#S2.SS1.p2.7),[§3](https://arxiv.org/html/2607.28989#S3.SS0.SSS0.Px3.p1.4),[4th item](https://arxiv.org/html/2607.28989#Sx4.I1.i4.p1.1)\.
- \[37\]Y\. LeCun, L\. Bottou, Y\. Bengio, and P\. Haffner\(1998\)Gradient\-based learning applied to document recognition\.Proceedings of the IEEE86\(11\),pp\. 2278–2324\.External Links:[Document](https://dx.doi.org/10.1109/5.726791)Cited by:[§1](https://arxiv.org/html/2607.28989#S1.p10.3),[§2\.1](https://arxiv.org/html/2607.28989#S2.SS1.p1.18),[§2\.1](https://arxiv.org/html/2607.28989#S2.SS1.p2.7),[§3](https://arxiv.org/html/2607.28989#S3.SS0.SSS0.Px1.p1.3),[4th item](https://arxiv.org/html/2607.28989#Sx4.I1.i4.p1.1)\.
- \[38\]Y\. LeCun, S\. Chopra, R\. Hadsell, M\. Ranzato, and F\. J\. Huang\(2006\)A tutorial on energy\-based learning\.InPredicting Structured Data,External Links:[Link](http://yann.lecun.com/exdb/publis/pdf/lecun-06.pdf)Cited by:[§1](https://arxiv.org/html/2607.28989#S1.p2.1),[§2\.8](https://arxiv.org/html/2607.28989#S2.SS8.p1.1)\.
- \[39\]Q\. Li, Z\. Han, and X\. Wu\(2018\)Deeper insights into graph convolutional networks for semi\-supervised learning\.InProceedings of the AAAI Conference on Artificial Intelligence,Vol\.32,pp\. 3538–3545\.External Links:[Document](https://dx.doi.org/10.1609/aaai.v32i1.11604)Cited by:[§3](https://arxiv.org/html/2607.28989#S3.SS0.SSS0.Px7.p3.8)\.
- \[40\]S\. McCallum, K\. Arora, and J\. Foster\(2025\)Reversible deep equilibrium models\.arXiv preprint arXiv:2509\.12917\.External Links:[Document](https://dx.doi.org/10.48550/arXiv.2509.12917),[Link](https://arxiv.org/abs/2509.12917)Cited by:[§4](https://arxiv.org/html/2607.28989#S4.p5.1)\.
- \[41\]T\. Miyato, T\. Kataoka, M\. Koyama, and Y\. Yoshida\(2018\)Spectral normalization for generative adversarial networks\.InInternational Conference on Learning Representations \(ICLR\),External Links:[Link](https://openreview.net/forum?id=B1QRgziT-)Cited by:[§2\.14](https://arxiv.org/html/2607.28989#S2.SS14.p1.8)\.
- \[42\]R\. Pascanu, T\. Mikolov, and Y\. Bengio\(2013\)On the difficulty of training recurrent neural networks\.InProceedings of the 30th International Conference on Machine Learning,Proceedings of Machine Learning Research, Vol\.28,pp\. 1310–1318\.External Links:[Link](https://proceedings.mlr.press/v28/pascanu13.html)Cited by:[§2\.12](https://arxiv.org/html/2607.28989#S2.SS12.p1.10),[§2\.7](https://arxiv.org/html/2607.28989#S2.SS7.p1.4)\.
- \[43\]B\. Poole, S\. Lahiri, M\. Raghu, J\. Sohl\-Dickstein, and S\. Ganguli\(2016\)Exponential expressivity in deep neural networks through transient chaos\.InAdvances in Neural Information Processing Systems \(NeurIPS\),Vol\.29\.External Links:[Link](https://papers.nips.cc/paper/6322-exponential-expressivity-in-deep-neural-networks-through-transient-chaos)Cited by:[§2\.10](https://arxiv.org/html/2607.28989#S2.SS10.SSS0.Px1.p4.4)\.
- \[44\]R\. P\. N\. Rao and D\. H\. Ballard\(1999\)Predictive coding in the visual cortex: a functional interpretation of some extra\-classical receptive\-field effects\.Nature Neuroscience2,pp\. 79–87\.External Links:[Document](https://dx.doi.org/10.1038/4580)Cited by:[§1](https://arxiv.org/html/2607.28989#S1.p5.1),[§2\.1](https://arxiv.org/html/2607.28989#S2.SS1.SSS0.Px2.p1.1)\.
- \[45\]P\. J\. Rousseeuw\(1987\)Silhouettes: a graphical aid to the interpretation and validation of cluster analysis\.Journal of Computational and Applied Mathematics20,pp\. 53–65\.External Links:[Document](https://dx.doi.org/10.1016/0377-0427%2887%2990125-7)Cited by:[§3](https://arxiv.org/html/2607.28989#S3.SS0.SSS0.Px2.p1.4)\.
- \[46\]D\. E\. Rumelhart, G\. E\. Hinton, and R\. J\. Williams\(1986\)Learning representations by back\-propagating errors\.Nature323,pp\. 533–536\.External Links:[Document](https://dx.doi.org/10.1038/323533a0)Cited by:[§1](https://arxiv.org/html/2607.28989#S1.p2.1)\.
- \[47\]Y\. Saad and M\. H\. Schultz\(1986\)GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems\.SIAM Journal on Scientific and Statistical Computing7\(3\),pp\. 856–869\.External Links:[Document](https://dx.doi.org/10.1137/0907058)Cited by:[§1](https://arxiv.org/html/2607.28989#S1.p9.1),[§2\.1](https://arxiv.org/html/2607.28989#S2.SS1.SSS0.Px3.p1.3)\.
- \[48\]F\. Scarselli, M\. Gori, A\. C\. Tsoi, M\. Hagenbuchner, and G\. Monfardini\(2009\)The graph neural network model\.IEEE Transactions on Neural Networks20\(1\),pp\. 61–80\.External Links:[Document](https://dx.doi.org/10.1109/TNN.2008.2005605)Cited by:[§1](https://arxiv.org/html/2607.28989#S1.p2.1)\.
- \[49\]P\. Sen, G\. Namata, M\. Bilgic, L\. Getoor, B\. Galligher, and T\. Eliassi\-Rad\(2008\)Collective classification in network data\.AI Magazine29\(3\),pp\. 93–106\.External Links:[Document](https://dx.doi.org/10.1609/aimag.v29i3.2157)Cited by:[§1](https://arxiv.org/html/2607.28989#S1.p6.1),[§2\.1](https://arxiv.org/html/2607.28989#S2.SS1.p1.18),[§2\.1](https://arxiv.org/html/2607.28989#S2.SS1.p2.27),[§3](https://arxiv.org/html/2607.28989#S3.SS0.SSS0.Px5.p1.1),[4th item](https://arxiv.org/html/2607.28989#Sx4.I1.i4.p1.1)\.
- \[50\]N\. Srivastava, G\. Hinton, A\. Krizhevsky, I\. Sutskever, and R\. Salakhutdinov\(2014\)Dropout: a simple way to prevent neural networks from overfitting\.Journal of Machine Learning Research15\(56\),pp\. 1929–1958\.External Links:[Link](https://jmlr.org/papers/v15/srivastava14a.html)Cited by:[§2\.4](https://arxiv.org/html/2607.28989#S2.SS4.p1.6)\.
- \[51\]L\. van der Maaten and G\. Hinton\(2008\)Visualizing data using t\-SNE\.Journal of Machine Learning Research9\(86\),pp\. 2579–2605\.External Links:[Link](https://www.jmlr.org/papers/v9/vandermaaten08a.html)Cited by:[§3](https://arxiv.org/html/2607.28989#S3.SS0.SSS0.Px2.p1.4)\.
- \[52\]A\. Vaswani, N\. Shazeer, N\. Parmar, J\. Uszkoreit, L\. Jones, A\. N\. Gomez, L\. Kaiser, and I\. Polosukhin\(2017\)Attention is all you need\.InAdvances in Neural Information Processing Systems \(NeurIPS\),Vol\.30,pp\. 5998–6008\.External Links:[Link](https://papers.nips.cc/paper/7181-attention-is-all-you-need)Cited by:[§1](https://arxiv.org/html/2607.28989#S1.p1.1),[§1](https://arxiv.org/html/2607.28989#S1.p4.1),[§2\.6](https://arxiv.org/html/2607.28989#S2.SS6.p1.30),[§4](https://arxiv.org/html/2607.28989#S4.p4.1)\.
- \[53\]P\. Veličković, G\. Cucurull, A\. Casanova, A\. Romero, P\. Liò, and Y\. Bengio\(2018\)Graph attention networks\.InInternational Conference on Learning Representations \(ICLR\),External Links:[Link](https://openreview.net/forum?id=rJXMpikCZ)Cited by:[§1](https://arxiv.org/html/2607.28989#S1.p4.1),[§2\.1](https://arxiv.org/html/2607.28989#S2.SS1.p2.7),[§2\.3](https://arxiv.org/html/2607.28989#S2.SS3.p1.2),[§3](https://arxiv.org/html/2607.28989#S3.SS0.SSS0.Px4.p1.2),[§3](https://arxiv.org/html/2607.28989#S3.SS0.SSS0.Px5.p1.1),[Table 1](https://arxiv.org/html/2607.28989#S3.T1.8.4.2),[Table 2](https://arxiv.org/html/2607.28989#S3.T2),[Table 2](https://arxiv.org/html/2607.28989#S3.T2.14.11.1.3),[§4](https://arxiv.org/html/2607.28989#S4.p4.1)\.
- \[54\]H\. F\. Walker and P\. Ni\(2011\)Anderson acceleration for fixed\-point iterations\.SIAM Journal on Numerical Analysis49\(4\),pp\. 1715–1735\.External Links:[Document](https://dx.doi.org/10.1137/10078356X)Cited by:[§2\.1](https://arxiv.org/html/2607.28989#S2.SS1.SSS0.Px3.p1.3),[§2\.9](https://arxiv.org/html/2607.28989#S2.SS9.p1.3)\.
- \[55\]Y\. Wang, Y\. Sun, Z\. Liu, S\. E\. Sarma, M\. M\. Bronstein, and J\. M\. Solomon\(2019\)Dynamic graph CNN for learning on point clouds\.ACM Transactions on Graphics38\(5\),pp\. 146\.External Links:[Document](https://dx.doi.org/10.1145/3326362)Cited by:[§2\.5](https://arxiv.org/html/2607.28989#S2.SS5.p1.7)\.
- \[56\]H\. R\. Wilson and J\. D\. Cowan\(1972\)Excitatory and inhibitory interactions in localized populations of model neurons\.Biophysical Journal12\(1\),pp\. 1–24\.External Links:[Document](https://dx.doi.org/10.1016/S0006-3495%2872%2986068-5)Cited by:[§1](https://arxiv.org/html/2607.28989#S1.p5.1),[§2\.1](https://arxiv.org/html/2607.28989#S2.SS1.SSS0.Px2.p1.1)\.
- \[57\]K\. Xu, W\. Hu, J\. Leskovec, and S\. Jegelka\(2019\)How powerful are graph neural networks?\.InInternational Conference on Learning Representations \(ICLR\),External Links:[Link](https://openreview.net/forum?id=ryGs6iA5Km)Cited by:[§3](https://arxiv.org/html/2607.28989#S3.SS0.SSS0.Px4.p1.2),[Table 1](https://arxiv.org/html/2607.28989#S3.T1.9.5.2)\.
- \[58\]Z\. Yang, W\. W\. Cohen, and R\. Salakhutdinov\(2016\)Revisiting semi\-supervised learning with graph embeddings\.InProceedings of the 33rd International Conference on Machine Learning,Proceedings of Machine Learning Research, Vol\.48,pp\. 40–48\.External Links:[Link](https://proceedings.mlr.press/v48/yanga16.html)Cited by:[§1](https://arxiv.org/html/2607.28989#S1.p6.1),[§2\.1](https://arxiv.org/html/2607.28989#S2.SS1.p1.18),[§2\.1](https://arxiv.org/html/2607.28989#S2.SS1.p2.27),[§3](https://arxiv.org/html/2607.28989#S3.SS0.SSS0.Px5.p1.1),[Table 2](https://arxiv.org/html/2607.28989#S3.T2),[4th item](https://arxiv.org/html/2607.28989#Sx4.I1.i4.p1.1)\.