Paper deep dive
Multi-Attractor GNNs: Set-Valued Expressivity Beyond Unique Equilibria
Jialin Liu
Intelligence
Status: not_run | Model: - | Prompt: - | Confidence: 0%
Entities (0)
Relation Signals (0)
No relation signals yet.
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Recurrent and equilibrium graph neural networks (GNNs) often enforce a unique fixed point or use one training target per graph. Yet many combinatorial and scientific problems admit multiple valid solutions, with no preferred one. A designated target can then impose an arbitrary selection rule. For tasks invariant to node relabeling, a symmetric graph may have a symmetric solution set but no symmetric solution. We show that multiple equilibria enable one weight-tied message-passing GNN to represent set-valued equivariant maps: different initializations approach different valid solutions. Under stated regularity assumptions, we first construct globally Lipschitz, permutation-equivariant dynamics that converge almost surely to valid solutions and reach every solution branch with positive probability. We then establish approximate realization by recurrent message passing with continuous component maps, with arbitrarily small update and limiting errors and arbitrarily high probability. This goes beyond standard universality arguments: although message passing alone cannot distinguish symmetric nodes, the evolving state keeps nodes distinguishable at every finite step without auxiliary node identifiers. Such dynamics can be learned without solution labels using problem-specific energies. On Ising ground states, structural module detection in protein graphs, and chemical reaction steady states, the learned updates produce multiple high-quality predictions with high numerical convergence rates. They achieve better average solution quality than the tested unique-equilibrium, single-target, and feedforward baselines, while remaining competitive with much larger diffusion-based solvers.
Tags
Links
- Source: https://arxiv.org/abs/2609.35274v1
- Canonical: https://arxiv.org/abs/2609.35274v1
Trouble viewing inline? Open PDF directly →
Full Text
219,297 characters extracted from source content.
Expand or collapse full text
Multi-Attractor GNNs: Set-Valued Expressivity Beyond Unique Equilibria Jialin Liu Address: (JL) School of Data, Mathematical, and Statistical Sciences, University of Central Florida, Orlando, FL 32826. Email address: jialin.liu@ucf.edu Date: September 28, 2026 Abstract. Recurrent and equilibrium graph neural networks (GNNs) are often designed to have a unique fixed point or trained against a single target per graph. Yet many combinatorial and scientific problems are inherently set-valued: several valid solutions coexist for the same graph, and nothing in the problem singles one out. Training against a single designated target can then make the model fit an arbitrary selection rule rather than the full solution set. For example, when a task is invariant to node relabeling, a symmetric graph has a symmetric solution set, yet may admit no symmetric solution; isolating a single target then imposes an arbitrary symmetry-breaking choice. In this paper, we show that multiple equilibria of recurrent GNNs are a resource rather than a defect: one weight-tied message-passing GNN can represent such set-valued equivariant maps through its attractor landscape, where different initializations approach different valid solutions. Under the stated regularity assumptions, we establish this expressive power in two steps. First, we prove the existence of globally Lipschitz, permutation-equivariant multi-attractor dynamics that converge almost surely to valid solutions while reaching every solution branch with positive probability. Second, we establish their approximate realization by recurrent message passing with continuous component maps, with arbitrarily small update and limiting errors and arbitrarily high probability. This goes beyond standard universality arguments, since message passing cannot by itself distinguish symmetric nodes: we show that the evolving state keeps nodes distinguishable at every finite step, without auxiliary node identifiers. In practice, such dynamics can be learned without solution labels from problem-specific energies. On three scientific tasks (ground states of Ising models, structural module detection in protein graphs, and steady states of chemical reaction networks), the learned updates produce multiple high-quality predictions and high numerical convergence rates, with better average solution quality than the tested unique-equilibrium, single-target, and feedforward baselines, while remaining competitive with far larger diffusion-based solvers. 1. Introduction Recurrent graph neural networks (GNNs) repeatedly apply a shared message-passing update, allowing computation to deepen without introducing new parameters at every step [35]. This iterative structure provides a natural mechanism for capturing long-range dependencies, while equilibrium-based variants support implicit differentiation without storing the full sequence of intermediate states [24, 6, 67]. For an n-node graph G=(A,X)G=(A,X), A encodes connectivity and edge attributes, and XiX_i contains node i’s features. A prediction y∈ℝn×qy ^n× q assigns each node a q-dimensional output yiy_i. A feedforward GNN computes y=GNNθ(G)y=GNN_θ(G) directly; a recurrent predictor updates a state yt∈ℝn×qy_t ^n× q: yt+1=GNNθ(G,yt),t=0,1,2,….y_t+1=GNN_θ(G,y_t), t=0,1,2,…. From a deterministic or random y0y_0, each step receives (Xi,(yt)i)(X_i,(y_t)_i) at node i and exchanges messages along edges to produce (yt+1)i(y_t+1)_i. The graph and parameters θ remain fixed. The usual goal is a limit y∞y_∞ approximating a prescribed target f(G)f(G). An equilibrium, or fixed point, is a state unchanged by the update. Implicit GNNs often seek well-posedness: existence and uniqueness of this state for each input [24, 39, 47]. A standard sufficient condition is global contraction: each update shrinks the distance between any two states by a common factor κ<1κ<1. Iteration then converges to the same fixed point from every initialization [53]. Strong monotonicity of an associated operator offers another route to uniqueness [63, 6]; uniqueness alone need not ensure convergence of direct iteration. Single-target training instead encourages one designated output without guaranteeing uniqueness. Contraction rules out initialization-dependent solutions regardless of the training objective. Many graph or combinatorial tasks admit multiple valid solutions for the same input. Selecting one target can discard meaningful alternatives or conflict with input symmetries. For example: • Two coloring. Consider proper two-coloring of the four-node cycle 1−2−3−4−11-2-3-4-1, with identical node features: each node receives a color in 0,1\0,1\, and adjacent nodes must have different colors. There are exactly two valid assignments, (0,1,0,1)(0,1,0,1) and (1,0,1,0)(1,0,1,0), which are exchanged by a one-step rotation that leaves the input graph unchanged. A vanilla GNN receiving only the graph must assign identical outputs to all nodes and therefore cannot produce either coloring. Symmetry breaking can distinguish these nodes and enable a valid coloring, but single-target training still arbitrarily favors one equally valid solution. • Nonlinear equilibria. Chemical reaction networks are represented as bipartite graphs connecting species to participating reactions, with reaction parameters as attributes; the task is to predict stationary species concentrations. Even one dynamic species can admit multiple solutions for the same graph and parameters: a Schlögl model has steady-state equation 0=8−14c+7c2−c30=8-14c+7c^2-c^3 and three positive solutions c∈1,2,4c∈\1,2,4\. All three satisfy the same physical equation, which selects no preferred target. Discovering multiple steady states can reveal coexisting stable states and help explain switching between them. These examples motivate recovering multiple valid solutions per graph without a preferred target. Recurrent GNNs provide a natural mechanism: initializations y0y_0 and y0′y_0 can converge to different limits y∞y_∞ and y∞′y_∞ under the same learned update (Figure 1). Only the initial node states change between runs; G and θ remain fixed. We seek convergence to a valid solution within each run and access to the graph’s different solutions across runs. Figure 1. Different initializations can lead to different equilibria under the same recurrent GNN update, with the graph and model parameters held fixed. Building on the set-valued dynamical viewpoint of Jore and Liu [30], we study permutation-equivariant representation and message-passing realization with recurrent GNNs. Our contributions address three questions. • Can one shared update represent many solutions? Under the stated assumptions, we prove that a single recurrent message-passing GNN can converge arbitrarily close to the solution set, with arbitrarily high probability over graphs and initializations (Theorems 1 and 2). • Is a unique equilibrium merely less diverse? Under suitable conditions, every continuous, equivariant, contractive update misses all valid solutions with positive probability, even on almost-surely asymmetric graphs. This applies to MP-GNNs of arbitrary finite depth and width, regardless of training (Theorem 3). • Can such dynamics be learned without solution labels, and do they pay off? On Ising ground states, protein module detection, and chemical steady states, energy or residual training from random starts yields multiple high-quality predictions and high numerical convergence rates. Average solution quality exceeds the tested single-target, unique-equilibrium, and feedforward baselines and remains competitive with much larger diffusion-based solvers. 2. Main results Existing expressivity theory studies single-valued invariant or equivariant maps [66, 33, 1]. Can one GNN instead represent solution sets through initialization-dependent limits? We establish existence, illustrate the mechanism, and identify an obstruction under contraction. 2.1. Can Recurrent GNNs Represent Set-Valued Equivariant Maps? To make the question precise, let Yn=ℝn×qY_n=R^n× q be the space of node-level outputs. Let G be a family of n-node graphs that is closed under node permutations. A permutation π acts on a graph by relabeling its nodes, and on YnY_n by permuting rows. A finite solution map assigns to each G∈G a set (G)=f1(G),…,fM(G)⊆Yn,S(G)=\f_1(G),…,f_M(G)\ Y_n, where the branches faf_a may coincide on some G. We require the map to be set-valued equivariant, (πG)=π(G).S(π G)=π\,S(G). Individual branches need not be equivariant. In the two-coloring example (Section 1), a one-node rotation preserves the graph but exchanges (0,1,0,1)(0,1,0,1) and (1,0,1,0)(1,0,1,0). The solution set is unchanged; selecting either coloring alone violates equivariance. A recurrent GNN represents solutions through the limits of its update yt+1=GNNθ(G,yt)y_t+1=GNN_θ(G,y_t) with y0∼μny_0 _n, where μn _n is an initialization distribution on YnY_n with a density and full support, such as i.i.d. standard Gaussian entries. We collect the limits reached with positive probability in the attractor set θ(G):=s∈Yn:ℙy0∼μn[limt→∞GNNθt(G,y0)=s]>0,A_θ(G):= \s∈ Y_n:P_y_0 _n [ _t→∞GNN_θ^\,t(G,y_0)=s ]>0 \, where GNNθt(G,⋅)GNN_θ^\,t(G,·) denotes t repeated applications of the block. Ideally, one shared update satisfies, for every G: (i) almost every trajectory ytt\y_t\_t converges to a limit in (G)S(G); and (i) every s∈(G)s (G) is reached with positive probability. Together, these ensure θ(G)=(G),A_θ(G)=S(G), with no probability mass on nonconvergent trajectories or invalid limits. Theorem 1 constructs a globally Lipschitz, equivariant operator satisfying (i)–(i) on suitably regular graphs. Theorem 2 realizes its dynamics approximately by recurrent message passing: limits approach the solution set with arbitrarily high joint probability, and each branch with positive joint probability over graphs and initializations. Theorem 1 (Existence of equivariant multi-attractor dynamics). Let G be a permutation-invariant family of n-node graphs, and let (G)=f1(G),…,fM(G)⊆Yn=ℝn×qS(G)=\f_1(G),…,f_M(G)\ Y_n=R^n× q be set-valued equivariant, with branch labels as in Definition 1. Let ⋆G_ consist of the graphs where the branches are locally Lipschitz and which branches coincide does not change under sufficiently small graph perturbations. Then there exists a single operator T:(G,y)↦TG(y)T:(G,y) T_G(y) that is jointly globally Lipschitz and permutation equivariant, TπG(πy)=πTG(y)T_π G(π y)=π T_G(y). For every G∈⋆G _ , the iteration yt+1=TG(yt)y_t+1=T_G(y_t), initialized from any distribution μn _n with a density and full support on YnY_n, satisfies yt→y∞∈(G) y_t→ y_∞ (G) almost surely, surely, ℙy0∼μn(y∞=s)>0 _y_0 _n(y_∞=s)>0 for every s∈(G). every s (G). Thus Lipschitz regularity, which bounds sensitivity without requiring contraction, permits multiple solutions and convergent trajectories (Appendix A). We next realize this behavior by message passing. Message-passing GNNs. For a graph G=(A,X)G=(A,X) and state y∈Yny∈ Y_n, define a GNN block GNNθ(G,y)GNN_θ(G,y) with finitely many layers and finite hidden dimensions: (1) hi(0) h_i^(0) =enc(Xi,yi), =enc(X_i,y_i), hi(ℓ+1) h_i^( +1) =Uℓ(hi(ℓ),∑j=1nMℓ(hi(ℓ),hj(ℓ),Aji)), =U_ \! (h_i^( ), _j=1^nM_ (h_i^( ),h_j^( ),A_ji) ), [GNNθ(G,y)]i [GNN_θ(G,y)]_i =O(hi(L),∑j=1nRL(hj(L))), =O\! (h_i^(L), _j=1^nR_L(h_j^(L)) ), where ℓ=0,…,L−1 =0,…,L-1 and i=1,…,ni=1,…,n. The encoder and component maps Mℓ,Uℓ,RL,OM_ ,U_ ,R_L,O are arbitrary continuous maps between finite-dimensional Euclidean spaces, shared across nodes. The edge features AjiA_ji encode edge attributes and edge presence; messages can be chosen to vanish on absent edges. Our theory uses arbitrary continuous component maps, rather than a fixed MLP parameterization; its guarantee is representational, not a guarantee for training. Theorem 2 (Approximate realization by a recurrent MP-GNN). Let S and ⋆G_ be as in Theorem 1, and assume q≥2q≥ 2. Let G follow a tight Borel probability on ⋆G_ , and independently draw y0∼μny_0 _n, where μn _n has a density and full support on YnY_n. Let TGT_G be the ideal operator in Theorem 1. For any ε,ρ>0 ,ρ>0 and δ∈(0,1)δ∈(0,1), there exists a single GNN of the form (1) such that, with probability at least 1−δ1-δ over (G,y0)(G,y_0), the iteration yt+1=GNNθ(G,yt)y_t+1=GNN_θ(G,y_t) converges to a limit y∞y_∞ satisfying supt≥0‖GNNθ(G,yt)−TG(yt)‖<ρanddist(y∞,(G))<ε. _t≥ 0 \|GNN_θ(G,y_t)-T_G(y_t) \|<ρ~~~~and~~~~dist(y_∞,S(G))< . Moreover, every solution branch is reached within ε with positive joint probability: ℙ(yt→y∞,∥y∞−fa(G)∥<ε)>0,a=1,…,M.P\! (y_t→ y_∞,\ \|y_∞-f_a(G)\|< )>0, a=1,…,M. The branch guarantee is joint over graphs and initializations, not per graph. Vanilla GNNs lack universality for continuous equivariant operators [66]; random initialization can provide unique node identifiers to overcome this [1]. But subsequent states depend on the graph and each other, even from i.i.d. y0y_0. Do they still distinguish nodes? Our proof establishes that independence is unnecessary: the ideal trajectories zt:=TGt(y0)z_t:=T_G^t(y_0) retain distinct node states at every finite time, ℙ(zt,i≠zt,j for all i≠j and all t∈ℕ0)=1.P\! (z_t,i≠ z_t,j for all i≠ j and all t _0 )=1. These evolving identifiers enable recurrent message passing to approximate the ideal dynamics (Appendix B). 2.2. Toy Example: Single-Source Graph Diffusion On a problem with exactly three solutions, we test whether one trained recurrent GNN approaches different solutions from different random initial states and recovers all three across runs. We compare against contraction and single-target supervision. Problem. The input is a three-node undirected graph G=(A,0)G=(A,0) with independent edge weights w12,w23,w31∼Unif[0.5,1.5]w_12,w_23,w_31 [0.5,1.5] and no distinguishing static node features. Here A is the symmetric weighted adjacency matrix with zero diagonal. Define LG=diag(A)−AL_G=diag(A1)-A and HG=I3+LGH_G=I_3+L_G. We seek a source indicator b∈0,13b∈\0,1\^3 selecting exactly one node, ⊤b=11 b=1, and its diffusion field u∈ℝ3u ^3 satisfying HGu=bH_Gu=b. Since HGH_G is positive definite, each source choice determines one diffusion field, giving exactly three solutions: (2) (G)=[HG−1ekek]:k=1,2,3⊂ℝ3×2,S(G)= \[H_G^-1e_k\ \ e_k]:k=1,2,3 \ ^3× 2, where eke_k is the kkth standard basis vector. The source is an output to discover, not an input supplied to the model. Figure 2 illustrates one graph and its three solutions. Model and training. The recurrent state yt=[utbt]∈ℝ3×2y_t=[u_t\ b_t] ^3× 2 stores two values per node. From six independent standard Gaussian entries, we iterate yt+1=GNNθ(G,yt)y_t+1=GNN_θ(G,y_t) using a shared block with two message-passing layers, following (1). For each fixed graph and model, we seek convergent trajectories reaching each solution in (G)S(G) with positive probability across initializations. We train using the energy function (3) EG(u,b)=13‖HGu−b‖22+∑i=13bi(1−bi)+(⊤b−1)2,b∈[0,1]3.E_G(u,b)= 13\|H_Gu-b\|_2^2+ _i=1^3b_i(1-b_i)+(1 b-1)^2, b∈[0,1]^3. The three terms are all nonnegative and EG(u,b)=0E_G(u,b)=0 holds exactly on (G)S(G). Training uses eight random starts per graph and minimizes the terminal energy G,y0[EG(uT,bT)]E_G,y_0[E_G(u_T,b_T)] with T=10T=10. This objective permits different starts to approach different solutions without prescribing branch labels. Baselines. We compare two controls with the same architecture, changing either the update constraint or training target. The contractive GNN uses the same energy objective but projects its parameters throughout training to enforce Lipy(GNNθ(G,⋅))≤0.9Lip_y(GNN_θ(G,·))≤ 0.9, ensuring a unique fixed point for each graph. The single-target-supervised GNN has no contraction constraint. For each training graph, we uniformly select one of its three exact solutions and retain it across all initializations and training visits. Squared-error training asks all starts to predict that target without guaranteeing a unique equilibrium. Details are in Appendix D. Results. Figure 3 tracks the source coordinates btb_t of 128 trajectories on one held-out graph, with identical initial states across all three models. The multi-attractor model forms three clusters near the valid source indicators e1,e2,e3e_1,e_2,e_3, whereas the supervised and contractive controls each concentrate near a single point, visibly separated from all three. The diffusion coordinates utu_t exhibit the same qualitative pattern (Appendix D, Figure 4). Thus energy training can produce multiple attractors near different solutions, consistent with Theorems 1 and 2, while contraction or single-target supervision can concentrate trajectories far from every valid solution. Figure 2. One graph, three valid solutions to single-source diffusion. Left: the input graph, with edge weights labeled. Right: each choice of source gives a solution fk(G)=[HG−1ekek]f_k(G)=[H_G^-1e_k\ \ e_k]. Node labels show (ui,bi)(u_i,b_i), with diffusion values rounded to four decimal places; orange marks the selected source (bi=1b_i=1). All three outputs are valid for the same input, which specifies no source. Figure 3. Different initial states can approach different solutions under the same learned update. Each blue point represents one of 128 trajectories, plotted in source coordinates ((bt)1,(bt)2,(bt)3)((b_t)_1,(b_t)_2,(b_t)_3); magenta circles mark the three valid source indicators e1,e2,e3e_1,e_2,e_3. Columns show recurrent steps t=0,2,4,6,8t=0,2,4,6,8. The multi-attractor model (top) forms three clusters near these targets, whereas the single-target-supervised (middle) and contractive (bottom) controls each concentrate away from them. All rows use the same held-out graph (test graph 0), training seed 1, and initial states. 2.3. Can a Unique Equilibrium Miss Every Valid Solution? Do these failures reflect training or a structural limitation? A unique equilibrium cannot represent multiple solutions through different initializations, but can it always recover one valid solution? The following theorem answers no for continuous, permutation-equivariant, globally contractive updates: even with independently sampled edge weights, their equilibrium remains separated from every valid solution with positive probability, regardless of training or initialization. Theorem 3. Let △G_ be the family of weighted triangles in (2), with edge weights in [w−,w+][w_-,w_+], where 0<w−<w+<∞0<w_-<w_+<∞. Let ν be the graph law induced by sampling the three edge weights independently and uniformly from [w−,w+][w_-,w_+]. Consider any jointly continuous, permutation-equivariant operator F:△×Y3→Y3F:G_ × Y_3→ Y_3, where Y3=ℝ3×2Y_3=R^3× 2, satisfying (4) ‖FG(y)−FG(z)‖∗≤κ‖y−z‖∗(G∈△,y,z∈Y3)\|F_G(y)-F_G(z)\|_*≤κ\|y-z\|_* (G _ ,\ y,z∈ Y_3) for some fixed norm ∥⋅∥∗\|·\|_* and κ∈[0,1)κ∈[0,1). This includes any MP-GNN of (1) with continuous component maps satisfying (4), regardless of its finite depth and width. For every 0<ε<2/30< < 2/3, there exists δF,ε>0 _F, >0 such that, under any joint law of (G,Y0)(G,Y_0) with graph marginal ν, the recurrence Yt+1=FG(Yt)Y_t+1=F_G(Y_t) converges to a limit Y∞Y_∞ and ℙ(dist(Y∞,(G))>ε)≥δF,ε,P\! (dist (Y_∞,S(G) )> )≥ _F, , where distance is measured in the Euclidean norm. In particular, no such operator converges to a valid solution almost surely over independently weighted graphs. Random initialization distinguishes nodes even on symmetric graphs; Figure 3 shows its influence in the early iterations. Contraction eventually erases these differences: all trajectories approach one equilibrium. On the symmetric triangle, equivariance forces this equilibrium to preserve the graph’s symmetries, preventing single-source selection. For independently weighted graphs, the equilibrium need not be symmetric, but its continuous dependence on the graph extends the failure to nearby graphs with positive probability. The failure probability depends on the operator; no uniform lower bound across models is claimed. See Appendix C for the proof and extension. 3. Scientific applications We test whether energy or residual training learns multiple high-quality predictions and numerically convergent dynamics. These empirical tests, including scalar states outside Theorem 2’s q≥2q≥ 2 setting, do not certify asymptotic convergence or complete coverage. Baselines use the reported budgets; inference computation is not matched. 3.1. Ising Models Problem. The Ising model is a classical model of interacting magnetic spins and is a standard discrete optimization problem. In the zero-field setting, each node of the weighted graph G=(V,E,J)G=(V,E,J) carries a spin si∈−1,+1s_i∈\-1,+1\, with energy and ground states (5) (G)=argmins∈−1,+1|V|EG(s),EG(s)=−∑i,j∈EJijsisj.S(G)= arg\,min_s∈\-1,+1\^|V|E_G(s), E_G(s)=- _\i,j\∈ EJ_ijs_is_j. We use open triangular lattices with antiferromagnetic couplings (Ji,j<0J_i,j<0), which favor opposite neighboring spins. These preferences conflict on triangles, creating competing low-energy configurations; global spin reversal also preserves energy. This multiplicity naturally tests set-valued prediction: we seek diverse low-energy configurations per graph rather than one arbitrarily selected target. Models and baselines. Our multi-attractor model (Ours in the tables) uses yt+1=GNNθ(G,yt)y_t+1=GNN_θ(G,y_t) from y0∼(0,I)y_0 (0,I), refining continuous spin estimates and applying the sign function to yTy_T to predict discrete spins s s. Training minimizes [EG(yT)]E[E_G(y_T)] with T=10T=10 shared updates, encouraging low energy without prescribing a solution for each initialization. Testing uses 200 updates. We compare four baseline groups. (A) Single-target supervision minimizes MSE to one Gurobi solution per graph using the same recurrent architecture, testing the effect of a designated target on accuracy and diversity. (B) Unique-equilibrium models are energy-trained IGNNs [24] with contraction constraints and increasing hidden widths, comparing quality and multiplicity under enforced uniqueness at several capacities. (C) Feedforward models stack 1, 3, or 10 independently parameterized copies of our block with the same energy objective, comparing untied depth with shared recurrence. (D) Distributional models learn sampling distributions or generative processes, whereas ours learns a shared update with initialization-dependent solutions. We compare quality and diversity with Mean-field GNN, Annealed Bernoulli GNN, DIFUSCO, and VAG (VAG-CO) [32, 58, 59, 51] as alternative approaches to set-valued prediction. Appendix E provides all details. Table 1. Ising results on 500 test graphs with 20 candidates per graph (mean ± standard deviation over three training seeds). e¯ e: mean energy per node (lower is better); Hit10: percentage of candidates within 10% of the certified minimum energy; C200C_200: percentage of trajectories converged by step 200; D20D_20: mean number of distinct optimal configurations found per graph, counting global spin flips separately; G20G_20: percentage of graphs with at least two such configurations, up to numerical tolerance. Details in Appendix E. Method e¯↓ e Hit10 ↑ C200↑C_200 D20↑D_20 G20↑G_20 Ours −1.0794±0.0028-1.0794± 0.0028 83.98±0.9383.98± 0.93 99.96±0.0299.96± 0.02 2.16±0.232.16± 0.23 68.73±5.1968.73± 5.19 Single-target Supervision 0.2077±0.58780.2077± 0.5878 9.45±4.379.45± 4.37 10.11±3.3810.11± 3.38 0.10±0.050.10± 0.05 0.80±0.990.80± 0.99 IGNN adapter (H16) −0.4581±0.2538-0.4581± 0.2538 5.20±3.315.20± 3.31 100.00±0.00100.00± 0.00 0.03±0.020.03± 0.02 0.00±0.000.00± 0.00 IGNN adapter (H32) −0.7516±0.0185-0.7516± 0.0185 8.00±0.718.00± 0.71 100.00±0.00100.00± 0.00 0.06±0.000.06± 0.00 0.00±0.000.00± 0.00 IGNN adapter (H64) −0.7954±0.0143-0.7954± 0.0143 13.20±2.1413.20± 2.14 100.00±0.00100.00± 0.00 0.08±0.010.08± 0.01 0.00±0.000.00± 0.00 Feedforward (×1× 1) −0.9720±0.0035-0.9720± 0.0035 36.13±0.7536.13± 0.75 – 1.31±0.031.31± 0.03 36.67±1.2336.67± 1.23 Feedforward (×3× 3) −1.0515±0.0016-1.0515± 0.0016 69.07±0.6669.07± 0.66 – 0.95±0.950.95± 0.95 22.40±31.6822.40± 31.68 Feedforward (×10× 10) −1.0538±0.0012-1.0538± 0.0012 70.93±1.5270.93± 1.52 – 0.28±0.010.28± 0.01 0.00±0.000.00± 0.00 Mean-field GNN −0.9736±0.0026-0.9736± 0.0026 36.71±0.7836.71± 0.78 – 1.34±0.061.34± 0.06 37.80±1.3137.80± 1.31 Annealed Bernoulli GNN −0.9766±0.0011-0.9766± 0.0011 37.71±0.8637.71± 0.86 – 1.37±0.031.37± 0.03 38.27±0.6638.27± 0.66 DIFUSCO (matched) −0.9033±0.0023-0.9033± 0.0023 37.64±0.2337.64± 0.23 – 1.56±0.031.56± 0.03 45.87±0.0945.87± 0.09 DIFUSCO (large) −1.0362±0.0024-1.0362± 0.0024 69.82±1.6469.82± 1.64 – 2.60±0.102.60± 0.10 68.00±2.9568.00± 2.95 VAG (matched) −0.8893±0.0075-0.8893± 0.0075 17.18±2.0217.18± 2.02 – 0.65±0.110.65± 0.11 17.40±4.9817.40± 4.98 VAG (large) −0.8814±0.0547-0.8814± 0.0547 17.34±5.6917.34± 5.69 – 0.56±0.170.56± 0.17 14.00±5.7414.00± 5.74 Quality and multiple solutions. Our model achieves the lowest mean energy and highest Hit10 (Table 1): 83.98%83.98\% of candidates are within 10% of optimum. It finds 2.162.16 distinct optima per graph on average and at least two on 68.73%68.73\% of graphs. Numerical convergence reaches 99.96%99.96\% by step 200 and 100% by step 500 for the same trajectories (Appendix E). Single-target supervision has substantially poorer quality and diversity; the IGNN adapters converge reliably but return one configuration per graph with low Hit10. Deeper untied models improve quality but lag in Hit10 and distinct-optimum discovery. DIFUSCO (large) finds more optima (D20=2.60D_20=2.60), but has worse mean energy despite approximately 606 times as many parameters. Consistent with our theory, one shared update trained on energy can converge to multiple high-quality solutions, offering a favorable combination of convergence, quality, and diversity over single-target training or enforced uniqueness in these experiments. 3.2. Structural module detection in protein graphs Problem. We seek structural modules in protein graphs: groups more internally connected than expected from their node degrees. We replace graph classification on the 1,113 PROTEINS graphs [9, 42] with instance-wise modularity maximization [46]. For an undirected graph G=(V,E)G=(V,E) with node degree did_i and node assignments zi∈1,…,Kz_i∈\1,…,K\, the objective and solution set are QG(z)=1|E|∑i,j∈E[zi=zj]−∑k=1K(∑idi[zi=k]2|E|)2,(G)=argmaxz∈1,…,K|V|QG(z).Q_G(z)= 1|E| _\i,j\∈ E1[z_i=z_j]- _k=1^K ( _id_i1[z_i=k]2|E| )^2,~~S(G)= arg\,max_z∈\1,…,K\^|V|Q_G(z). We allow at most four modules (K=4K=4, with empty modules permitted). Partitions can have similar energies even modulo module-label permutations, motivating distinct low-energy predictions without one selected target. Models and baselines. Our shared recurrent update Pt+1=GNNθ(G,Pt)P_t+1=GNN_θ(G,P_t) refines a randomly initialized soft partition, where Pt,i∈ℝKP_t,i ^K contains node i’s assignment weights over K modules. After T iterations, each node selects its highest-weight module, zi=argmaxk(PT,i)kz_i= _k(P_T,i)_k. Training minimizes [eG(PT)]E[e_G(P_T)] with T=5T=5 over training graphs and initializations, where eGe_G is a differentiable soft-partition energy derived from modularity. At test time, we apply the learned update for 200 steps from 20 independent initializations per graph. The four baseline groups (Section 3.1) test training objectives, equilibrium uniqueness, weight sharing, and alternative solution generation. Appendix F details the soft energy, architectures, training, and evaluation protocols. Results. Our model has the highest mean modularity (0.55600.5560) and Hit10 (77.49%77.49\%) in Table 2, finding 5.835.83 distinct qualified partitions per graph on average and at least two on 80.78%80.78\% of graphs. Single-target supervision and unique-equilibrium IGNNs have poorer quality and diversity. Untied depth improves performance, but five independent blocks still lag in mean modularity, Hit10, and partition discovery. DIFUSCO (large) finds more qualified partitions (D20.10=11.02D^.10_20=11.02) but has lower mean modularity with approximately 155155 times as many parameters. Together with the convergence results in Appendix F (Table 10), these findings support our theoretical premise: energy training without prescribed targets lets one shared update recover multiple high-quality solutions from different initializations. Table 2. Structural module detection on 111 PROTEINS test graphs with 20 candidates per graph. Q¯G Q_G averages modularity over candidates and then graphs. For optimal modularity QG⋆Q_G , Hit10 is the percentage of candidates with relative gap (QG⋆−QG(z))/|QG⋆|≤0.10(Q_G -Q_G(z))/|Q_G |≤ 0.10. D20.10D^.10_20: mean number of distinct partitions per graph with absolute gap QG⋆−QG(z)≤0.10Q_G -Q_G(z)≤ 0.10; G20.10G^.10_20: percentage of graphs with at least two such partitions; C200C_200: percentage of trajectories numerically converged at t=200t=200. Results: mean ± standard deviation over training seeds 0/1/2 on the same test set. Full protocol in Appendix F. Method Q¯G↑ Q_G Hit10 ↑ C200↑C_200 D20.10↑D^.10_20 G20.10↑G^.10_20 Ours 0.5560±0.00240.5560± 0.0024 77.49±2.0677.49± 2.06 88.05±2.2388.05± 2.23 5.83±0.385.83± 0.38 80.78±0.4280.78± 0.42 Single-target Supervision 0.4911±0.01140.4911± 0.0114 41.28±7.1641.28± 7.16 50.81±10.8450.81± 10.84 1.11±0.081.11± 0.08 33.03±2.2533.03± 2.25 IGNN (H16) 0.4147±0.00790.4147± 0.0079 8.11±2.658.11± 2.65 100.00±0.00100.00± 0.00 0.20±0.040.20± 0.04 0.00±0.000.00± 0.00 IGNN (H32) 0.4196±0.00440.4196± 0.0044 11.71±2.2111.71± 2.21 100.00±0.00100.00± 0.00 0.22±0.030.22± 0.03 0.00±0.000.00± 0.00 IGNN (H64) 0.4235±0.01210.4235± 0.0121 8.41±2.368.41± 2.36 100.00±0.00100.00± 0.00 0.22±0.060.22± 0.06 0.00±0.000.00± 0.00 Feedforward (×1× 1) 0.4436±0.00480.4436± 0.0048 16.52±2.2216.52± 2.22 – 1.35±0.231.35± 0.23 30.93±6.2630.93± 6.26 Feedforward (×2× 2) 0.4846±0.01560.4846± 0.0156 29.50±3.6129.50± 3.61 – 2.09±0.262.09± 0.26 53.15±6.6253.15± 6.62 Feedforward (×3× 3) 0.5088±0.01080.5088± 0.0108 39.92±5.1739.92± 5.17 – 2.95±0.622.95± 0.62 71.47±9.7171.47± 9.71 Feedforward (×4× 4) 0.5285±0.00220.5285± 0.0022 51.35±1.6851.35± 1.68 – 4.24±0.514.24± 0.51 82.88±6.6282.88± 6.62 Feedforward (×5× 5) 0.5305±0.00440.5305± 0.0044 54.97±3.5954.97± 3.59 – 4.32±0.474.32± 0.47 78.68±4.4978.68± 4.49 Mean-field GNN 0.4020±0.01220.4020± 0.0122 8.53±2.018.53± 2.01 – 0.69±0.050.69± 0.05 16.22±0.0016.22± 0.00 Annealed Bernoulli GNN 0.4131±0.01030.4131± 0.0103 9.91±2.359.91± 2.35 – 0.80±0.010.80± 0.01 18.32±1.5318.32± 1.53 DIFUSCO (matched) 0.5115±0.01270.5115± 0.0127 41.85±7.9741.85± 7.97 – 8.47±1.208.47± 1.20 97.30±1.2797.30± 1.27 DIFUSCO (large) 0.5515±0.00140.5515± 0.0014 74.49±1.6674.49± 1.66 – 11.02±0.0311.02± 0.03 97.00±0.4297.00± 0.42 VAG (matched) 0.3952±0.00240.3952± 0.0024 8.05±0.148.05± 0.14 – 0.57±0.010.57± 0.01 12.01±0.8512.01± 0.85 VAG (large) 0.2577±0.02370.2577± 0.0237 0.44±0.310.44± 0.31 – 0.35±0.150.35± 0.15 7.21±1.477.21± 1.47 3.3. Chemical reaction-network multistationarity Problem. A chemical reaction network has s species i and R reactions r. The matrices α,β∈ℤ≥0R×sα,β _≥ 0^R× s give the numbers of molecules of species i consumed (αri _ri) and produced (βri _ri) by reaction r; κ∈ℝ>0Rκ _>0^R contains rate constants. The graph G=(α,β,κ)G=(α,β,κ) has one node per species and reaction, connected bidirectionally when αri+βri>0 _ri+ _ri>0. Edges carry (αri,βri)( _ri, _ri); reaction nodes encode logκr _r. Under mass-action kinetics, the species concentrations c∈ℝ>0sc _>0^s evolve according to c˙=fG(c)=Nv(c),Nir=βri−αri,vr(c)=κr∏i=1sciαri. c=f_G(c)=Nv(c), N_ir= _ri- _ri, v_r(c)= _r _i=1^sc_i _ri. We seek stationary solutions (G)=c∈ℝ>0s:fG(c)=0S(G)=\c _>0^s:f_G(c)=0\. Multiple positive solutions for one G constitute multistationarity, a set-valued prediction task. The solver seeks roots of fGf_G without requiring physical stability under the chemical ODE. Our kinetic instances use published network structures [68], with two numerically verified stationary witnesses per instance. Models and baselines. Species log concentrations follow the shared update zt+1=GNNθ(G,zt)z_t+1=GNN_θ(G,z_t), with output cT=expzTc_T= z_T. Training minimizes equation residuals without prescribed targets, using randomized unrolls of at most 32 steps; the main comparison uses T=200T=200. The four baseline groups (Section 3.1) are single-target MSE, physics-trained contractive IGNNs, untied feedforward depth controls, and distributional models. For continuous states, the latter comprise DDPM, DDIM, and DIS with matched and larger upstream-recipe graph adaptations; details are in Appendix G. Table 3. Chemistry results on 236 test graphs with 20 candidates per graph. E¯res E_ res: mean Huber penalty on the mass-action residual (lower is better); Hit10−3_10^-3: percentage of final states with normalized physical residual R≤10−3R≤ 10^-3; C: percentage of recurrent trajectories with next-step log-concentration RMS change below 10−310^-3 at T=200T=200; D20rawD_20 raw: mean number of distinct accepted prediction clusters per graph; G20rawG_20 raw: percentage of graphs with at least two. These count separated accepted predictions, not certified distinct roots. Entries: mean ± population standard deviation over training seeds 0/1/2. Method E¯res↓ E_ res Hit10−3↑_10^-3 C↑C D20raw↑D_20 raw G20raw↑G_20 raw Ours 0.0188 ± 0.0046 83.77 ± 3.44 91.99 ± 0.35 3.506 ± 0.314 41.67 ± 17.01 Single-target Supervis. 0.0373 ± 0.0059 18.78 ± 3.88 88.48 ± 4.29 1.371 ± 0.137 11.30 ± 0.72 IGNN (H16) 0.1214 ± 0.0823 0.00 ± 0.00 100.00 ± 0.00 0.000 ± 0.000 0.00 ± 0.00 IGNN (H32) 0.1180 ± 0.0709 0.00 ± 0.00 100.00 ± 0.00 0.000 ± 0.000 0.00 ± 0.00 IGNN (H64) 0.0960 ± 0.0056 0.00 ± 0.00 100.00 ± 0.00 0.000 ± 0.000 0.00 ± 0.00 Feedforward (×1× 1) 0.1082 ± 0.0081 0.00 ± 0.00 — 0.000 ± 0.000 0.00 ± 0.00 Feedforward (×2× 2) 0.0586 ± 0.0050 0.01 ± 0.02 — 0.003 ± 0.004 0.14 ± 0.20 Feedforward (×3× 3) 0.0441 ± 0.0065 0.11 ± 0.10 — 0.021 ± 0.019 0.28 ± 0.20 Feedforward (×4× 4) 0.0295 ± 0.0016 0.42 ± 0.35 — 0.076 ± 0.061 0.42 ± 0.35 Feedforward (×5× 5) 0.0255 ± 0.0022 0.58 ± 0.41 — 0.110 ± 0.078 0.71 ± 0.53 DDPM / matched 0.0576 ± 0.0047 0.03 ± 0.01 — 0.006 ± 0.002 0.00 ± 0.00 DDIM / matched 0.0732 ± 0.0055 0.01 ± 0.01 — 0.003 ± 0.002 0.00 ± 0.00 DIS / matched 0.9315 ± 0.0007 0.00 ± 0.00 — 0.000 ± 0.000 0.00 ± 0.00 DDPM / large 0.1197 ± 0.0040 0.01 ± 0.01 — 0.003 ± 0.002 0.00 ± 0.00 DDIM / large 0.1317 ± 0.0030 0.00 ± 0.00 — 0.000 ± 0.000 0.00 ± 0.00 DIS / large 0.0463 ± 0.0003 0.02 ± 0.02 — 0.004 ± 0.003 0.00 ± 0.00 Results. Our model achieves 83.77%83.77\% candidate hits and 3.5063.506 accepted prediction clusters per graph on average, with at least two on 41.67%41.67\% of graphs (Table 3). Numerical settling increases from 91.99%91.99\% at T=200T=200 to 96.15%96.15\% at T=5,000T=5,000; remaining oscillations and unresolved cases appear in Table 13. Single-target supervision reaches only 18.78%18.78\% hits; IGNNs settle in their hidden-state operators but yield no accepted predictions, illustrating that convergence alone does not ensure validity. Untied depth improves accuracy, but five independent stages reach only 0.58%0.58\% hits; tested distributional models also have low hit rates. These results demonstrate multiple residual-qualified predictions and high numerical settling rates. They do not certify distinct exact roots or complete root coverage. Baselines differ in objectives and inference budgets, so the comparison does not isolate recurrence at equal computational cost. 4. Conclusions Multi-attractor dynamics provide a mechanism for representing set-valued graph tasks with recurrent, weight-tied GNNs. Under the stated assumptions, one shared update can converge arbitrarily close to different valid solutions, while continuous, equivariant contraction can obstruct recovery of even one. Across three scientific tasks, energy or residual training produces multiple high-quality predictions and high numerical convergence rates without solution labels. These findings motivate preserving access to multiple equilibria when tasks admit multiple answers. This general framework opens opportunities across graph applications, with success depending on problem-specific objectives and optimization. Understanding how training shapes convergence and solution coverage remains an important direction. References [1] R. Abboud, I. I. Ceylan, M. Grohe, and T. Lukasiewicz (2020) The surprising power of graph neural networks with random node initialization. arXiv preprint arXiv:2010.01179. Cited by: Appendix H, Appendix H, §2.1, §2. [2] L. A. Adamic and E. Adar (2003) Friends and neighbors on the web. Social networks 25 (3), p. 211–230. Cited by: Appendix F. [3] C. Anil, A. Pokle, K. Liang, J. Treutlein, Y. Wu, S. Bai, J. Z. Kolter, and R. B. Grosse (2022) Path independent equilibrium models can better exploit test-time computation. Advances in Neural Information Processing Systems 35, p. 7796–7809. Cited by: Appendix H. [4] W. Azizian and marc lelarge (2021) Expressive power of invariant and equivariant graph neural networks. In International Conference on Learning Representations, External Links: Link Cited by: §B.2, Appendix H. [5] S. Bai, J. Z. Kolter, and V. Koltun (2019) Deep equilibrium models. Advances in neural information processing systems 32. Cited by: Appendix H. [6] J. Baker, Q. Wang, C. D. Hauck, and B. Wang (2023) Implicit graph neural networks: a monotone operator viewpoint. In International conference on machine learning, p. 1521–1548. Cited by: Appendix H, §1, §1. [7] A. Bansal, A. Schwarzschild, E. Borgnia, Z. Emam, F. Huang, M. Goldblum, and T. Goldstein (2022) End-to-end algorithm synthesis with recurrent networks: extrapolation without overthinking. Advances in Neural Information Processing Systems 35, p. 20232–20242. Cited by: Appendix H. [8] J. Berner, L. Richter, and K. Ullrich (2024) An optimal control perspective on diffusion-based generative modeling. Transactions on Machine Learning Research. External Links: ISSN 2835-8856, Link Cited by: 4th item, Appendix H. [9] K. M. Borgwardt, C. S. Ong, S. Schönauer, S. V. N. Vishwanathan, A. J. Smola, and H. Kriegel (2005) Protein function prediction via graph kernels. Bioinformatics. Cited by: §3.2. [10] P. Breiding and S. Timme (2018) HomotopyContinuation.jl: a package for homotopy continuation in Julia. In Mathematical Software – ICMS 2018, Lecture Notes in Computer Science, Vol. 10931. Cited by: Appendix H. [11] S. Brody, U. Alon, and E. Yahav (2022) How attentive are graph attention networks?. In International Conference on Learning Representations, External Links: Link Cited by: Appendix G. [12] C. Cai and Y. Wang (2018) A simple yet effective baseline for non-attributed graph classification. arXiv preprint arXiv:1811.03508. Cited by: Appendix F. [13] Q. Chen, T. Zhang, L. Yang, Q. Han, A. Wang, R. Sun, X. Luo, and T. Chang (2024) Symilo: a symmetry-aware learning framework for integer linear optimization. Advances in Neural Information Processing Systems 37, p. 24411–24434. Cited by: Appendix H. [14] Z. Chen, S. Villar, L. Chen, and J. Bruna (2019) On the equivalence between graph isomorphism testing and function approximation with gnns. Advances in neural information processing systems 32. Cited by: Appendix H. [15] K. Cho, B. Van Merriënboer, Ç. Gulçehre, D. Bahdanau, F. Bougares, H. Schwenk, and Y. Bengio (2014) Learning phrase representations using rnn encoder–decoder for statistical machine translation. In Proceedings of the 2014 conference on empirical methods in natural language processing (EMNLP), p. 1724–1734. Cited by: Appendix G. [16] G. Corso, L. Cavalleri, D. Beaini, P. Liò, and P. Veličković (2020) Principal neighbourhood aggregation for graph nets. Advances in neural information processing systems 33, p. 13260–13271. Cited by: Appendix F. [17] H. Dai, Z. Kozareva, B. Dai, A. Smola, and L. Song (2018) Learning steady-states of iterative algorithms over graphs. In International conference on machine learning, p. 1106–1114. Cited by: Appendix H. [18] V. P. Dwivedi, A. T. Luu, T. Laurent, Y. Bengio, and X. Bresson (2022) Graph neural networks with learnable structural and positional representations. In International Conference on Learning Representations, External Links: Link Cited by: Appendix F, Appendix H. [19] N. Dym, H. Lawrence, and J. W. Siegel (2024) Equivariant frames and the impossibility of continuous canonicalization. In International Conference on Machine Learning, Cited by: Appendix H. [20] P. E. Farrell, Á. Birkisson, and S. W. Funke (2015) Deflation techniques for finding distinct solutions of nonlinear partial differential equations. SIAM Journal on Scientific Computing 37 (4), p. A2026–A2045. Cited by: Appendix H. [21] F. Geerts and J. L. Reutter (2022) Expressiveness and approximation properties of graph neural networks. arXiv preprint arXiv:2204.04661. Cited by: §B.2, Appendix H. [22] A. Goel, D. Lim, H. Lawrence, S. Jegelka, and N. Huang (2026) Any-subgroup equivariant networks via symmetry breaking. arXiv preprint arXiv:2603.19486. Cited by: Appendix H. [23] L. Gonon, T. Meyer-Brandis, and N. Weber (2026) Universality and approximation rates of graph neural networks with random features. arXiv preprint arXiv:2607.26699. Cited by: Appendix H. [24] F. Gu, H. Chang, W. Zhu, S. Sojoudi, and L. El Ghaoui (2020) Implicit graph neural networks. Advances in neural information processing systems 33, p. 11984–11995. Cited by: 2nd item, 2nd item, Appendix H, §1, §1, §3.1. [25] Gurobi Optimization, LLC (2026) Gurobi Optimizer Reference Manual. External Links: Link Cited by: Appendix E. [26] A. Hagberg, P. J. Swart, and D. A. Schult (2007) Exploring network structure, dynamics, and function using networkx. Technical report Los Alamos National Laboratory (LANL). Cited by: Appendix F, Appendix F. [27] F. Hendriks, O. Rokoš, M. Doškář, M. G. Geers, and V. Menkovski (2025) Equivariant flow matching for symmetry-breaking bifurcation problems. arXiv preprint arXiv:2509.03340. Cited by: Appendix H. [28] J. Ho, A. Jain, and P. Abbeel (2020) Denoising diffusion probabilistic models. Advances in neural information processing systems 33, p. 6840–6851. Cited by: 4th item, Appendix H. [29] J. J. Hopfield (1982) Neural networks and physical systems with emergent collective computational abilities.. Proceedings of the national academy of sciences 79 (8), p. 2554–2558. Cited by: Appendix H. [30] C. Jore and J. Liu (2026) Bifurcation models: learning set-valued solution maps with weight-tied dynamics. arXiv preprint arXiv:2605.07277. Cited by: §A.3, Appendix E, Appendix H, §1. [31] S. Kaba and S. Ravanbakhsh (2023) Symmetry breaking and equivariant neural networks. arXiv preprint arXiv:2312.09016. Cited by: Appendix H. [32] N. Karalias and A. Loukas (2020) Erdos goes neural: an unsupervised learning framework for combinatorial optimization on graphs. Advances in Neural Information Processing Systems 33, p. 6659–6672. Cited by: 4th item, Appendix H, §3.1. [33] N. Keriven and G. Peyré (2019) Universal invariant and equivariant graph neural networks. Advances in neural information processing systems 32. Cited by: §B.2, Appendix H, §2. [34] H. Lawrence, V. Portilheiro, Y. Zhang, and S. Kaba (2025) Improving equivariant networks with probabilistic symmetry breaking. In International Conference on Learning Representations, Cited by: Appendix H. [35] Y. Li, D. Tarlow, M. Brockschmidt, and R. Zemel (2015) Gated graph sequence neural networks. arXiv preprint arXiv:1511.05493. Cited by: Appendix H, §1. [36] Z. Li, Q. Chen, and V. Koltun (2018) Combinatorial optimization with graph convolutional networks and guided tree search. Advances in neural information processing systems 31. Cited by: Appendix H. [37] D. Liben-Nowell and J. Kleinberg (2003) The link prediction problem for social networks. In Proceedings of the twelfth international conference on Information and knowledge management, p. 556–559. Cited by: Appendix F. [38] J. Liu, L. Ding, S. J. Osher, and W. Yin (2026) Expressive power of implicit models: rich equilibria and test-time scaling. In International Conference on Learning Representations, Vol. 2026, p. 17164–17212. Cited by: Appendix H, Lemma 3. [39] J. Liu, K. Kawaguchi, B. Hooi, Y. Wang, and X. Xiao (2021) Eignn: efficient infinite-depth graph neural networks. Advances in Neural Information Processing Systems 34, p. 18762–18773. Cited by: Appendix H, §1. [40] A. Loukas (2019) What graph neural networks cannot learn: depth vs width. arXiv preprint arXiv:1907.03199. Cited by: Appendix H. [41] H. Maron, H. Ben-Hamu, H. Serviansky, and Y. Lipman (2019) Provably powerful graph networks. Advances in neural information processing systems 32. Cited by: Appendix H. [42] C. Morris, N. M. Kriege, F. Bause, K. Kersting, P. Mutzel, and M. Neumann (2020) Tudataset: a collection of benchmark datasets for learning with graphs. arXiv preprint arXiv:2007.08663. Cited by: §3.2. [43] C. Morris, M. Ritzert, M. Fey, W. L. Hamilton, J. E. Lenssen, G. Rattan, and M. Grohe (2019) Weisfeiler and leman go neural: higher-order graph neural networks. In Proceedings of the AAAI conference on artificial intelligence, Cited by: Appendix H. [44] J. R. Munkres (2014) Topology. 2 edition, Pearson Education Limited, Harlow, England. External Links: ISBN 978-1-292-02362-5 Cited by: §B.2. [45] Y. Nandwani, D. Jindal, P. Singla, et al. (2020) Neural learning of one-of-many solutions for combinatorial problems in structured output spaces. arXiv preprint arXiv:2008.11990. Cited by: Appendix H. [46] M. E. Newman (2006) Modularity and community structure in networks. Proceedings of the national academy of sciences 103 (23), p. 8577–8582. Cited by: §3.2. [47] J. Park, J. Choo, and J. Park (2021) Convergent graph solvers. arXiv preprint arXiv:2106.01680. Cited by: Appendix H, §1. [48] H. Ramsauer, B. Schäfl, J. Lehner, P. Seidl, M. Widrich, T. Adler, L. Gruber, M. Holzleitner, M. Pavlović, G. K. Sandve, et al. (2020) Hopfield networks is all you need. arXiv preprint arXiv:2008.02217. Cited by: Appendix H. [49] F. Ricci-Tersenghi (2012) The bethe approximation for solving the inverse ising problem: a comparison with other inference methods. Journal of Statistical Mechanics: Theory and Experiment 2012 (08), p. P08015. Cited by: Appendix E. [50] E. Rosenbluth and M. Grohe (2026) Repetition makes perfect: recurrent graph neural networks match message passing limit. In Proceedings of the AAAI Conference on Artificial Intelligence, Cited by: Appendix H. [51] S. Sanokowski, W. Berghammer, S. Hochreiter, and S. Lehner (2023) Variational annealing on graphs for combinatorial optimization. Advances in Neural Information Processing Systems 36, p. 63907–63930. Cited by: 4th item, Appendix H, §3.1. [52] R. Sato, M. Yamada, and H. Kashima (2021) Random features strengthen graph neural networks. In Proceedings of the 2021 SIAM international conference on data mining (SDM), p. 333–341. Cited by: Appendix H, Appendix H. [53] F. Scarselli, M. Gori, A. C. Tsoi, M. Hagenbuchner, and G. Monfardini (2008) The graph neural network model. IEEE transactions on neural networks 20 (1), p. 61–80. Cited by: Appendix H, §1. [54] M. J. Schuetz, J. K. Brubaker, and H. G. Katzgraber (2022) Combinatorial optimization with physics-inspired graph neural networks. Nature Machine Intelligence 4 (4), p. 367–377. Cited by: Appendix H. [55] D. Selsam, M. Lamm, B. Bünz, P. Liang, L. de Moura, and D. L. Dill (2018) Learning a sat solver from single-bit supervision. arXiv preprint arXiv:1802.03685. Cited by: Appendix H. [56] T. E. Smidt, M. Geiger, and B. K. Miller (2021) Finding symmetry breaking order parameters with Euclidean neural networks. Physical Review Research 3, p. L012002. Cited by: Appendix H. [57] J. Song, C. Meng, and S. Ermon (2021) Denoising diffusion implicit models. In International Conference on Learning Representations, External Links: Link Cited by: 4th item, Appendix H. [58] H. Sun, E. K. Guha, and H. Dai (2022) Annealed training for combinatorial optimization on graphs. arXiv preprint arXiv:2207.11542. Cited by: 4th item, Appendix H, §3.1. [59] Z. Sun and Y. Yang (2023) Difusco: graph-based diffusion solvers for combinatorial optimization. Advances in neural information processing systems 36, p. 3706–3731. Cited by: 4th item, Appendix H, §3.1. [60] H. Tang, Z. Huang, J. Gu, B. Lu, and H. Su (2020) Towards scale-invariant graph-related problem solving by iterative homogeneous gnns. Advances in Neural Information Processing Systems 33, p. 15811–15822. Cited by: Appendix H. [61] J. Toenshoff, M. Ritzert, H. Wolf, and M. Grohe (2021) Graph neural networks for maximum constraint satisfaction. Frontiers in artificial intelligence 3, p. 580607. Cited by: Appendix H. [62] P. Veličković, R. Ying, M. Padovano, R. Hadsell, and C. Blundell (2019) Neural execution of graph algorithms. arXiv preprint arXiv:1910.10593. Cited by: Appendix H. [63] E. Winston and J. Z. Kolter (2020) Monotone operator equilibrium networks. Advances in neural information processing systems 33, p. 10718–10728. Cited by: Appendix H, §1. [64] D. Wu, L. Wang, and P. Zhang (2019) Solving statistical mechanics using variational autoregressive networks. Physical review letters 122 (8), p. 080602. Cited by: Appendix H. [65] Y. Xie and T. Smidt (2024) Equivariant symmetry breaking sets. Transactions on Machine Learning Research. Cited by: Appendix H. [66] K. Xu, W. Hu, J. Leskovec, and S. Jegelka (2018) How powerful are graph neural networks?. arXiv preprint arXiv:1810.00826. Cited by: §B.2, Appendix H, §2.1, §2. [67] Y. Yang, T. Liu, Y. Wang, Z. Huang, and D. Wipf (2025) Implicit vs unfolded graph neural networks. Journal of Machine Learning Research 26 (82), p. 1–46. Cited by: Appendix H, §1. [68] S. Yao, A. Sadeghimanesh, and M. England (2025) Understanding multistationarity of fully open reaction networks: s. yao, ah sadeghimanesh, m. england. Bulletin of Mathematical Biology 87 (12), p. 176. Cited by: Appendix G, §3.3. [69] D. Zhang, H. Dai, N. Malkin, A. C. Courville, Y. Bengio, and L. Pan (2023) Let the flows tell: solving graph combinatorial problems with gflownets. Advances in neural information processing systems 36, p. 11952–11969. Cited by: Appendix H. [70] T. Zhou, L. Lü, and Y. Zhang (2009) Predicting missing links via local information. The European Physical Journal B 71 (4), p. 623–630. Cited by: Appendix F. Appendix Contents Appendix A Ideal Operator Existence In the appendices, S(G)S(G) and nY_n denote (G)S(G) and YnY_n from the main text. Fix the number of nodes n. We represent a graph instance as G=(A,X)∈ℰn:=ℝn×n×de×ℝn×dv,G=(A,X) _n:=R^n× n× d_e×R^n× d_v, where A=(A(1),…,A(de))A=(A^(1),…,A^(d_e)) contains the edge features and X contains the node features. The output space is n=ℝn×q≅ℝm.Y_n=R^n× q ^m. Both spaces are equipped with their Euclidean norms. For a permutation π of the n nodes, let PπP_π denote its permutation matrix, with Pπei=eπ(i)P_πe_i=e_π(i). Relabeling the nodes of a graph and an output is defined by πG:=((PπA(r)Pπ⊤)r=1de,PπX),πy:=Pπy.π G:= ( (P_πA^(r)P_π )_r=1^d_e,P_πX ), π y:=P_πy. Since permutation matrices are orthogonal, these relabelings preserve Euclidean distances. Throughout, ⊆ℰnG _n is closed under node relabeling: G∈⟹πG∈G π G for every node permutation π. The domain G need not be bounded, closed, connected, or compact. A.1. Basic concepts and main result Definition 1 (Permutation-equivariant finite-branch solution map). A finite-branch solution map on G is represented, for some integer M≥1M≥ 1, by labeled branches f1,…,fM:→nf_1,…,f_M:G _n through S(G):=f1(G),…,fM(G).S(G):=\f_1(G),…,f_M(G)\. Distinct branch indices are allowed to have coincident values. Since S(G)S(G) is a set, if fi(G)=fj(G)f_i(G)=f_j(G) for i≠ji≠ j, their common value appears only once in S(G)S(G). We say that this representation is permutation equivariant if, for every node permutation π, there exists a permutation τπ _π of the branch indices 1,…,M\1,…,M\ such that fτπ(i)(πG)=πfi(G)for every G∈ and i∈1,…,M.f_ _π(i)(π G)=π f_i(G) every G and i∈\1,…,M\. Thus, relabeling the nodes may relabel the solution branches, but it does not change the underlying solution set. In particular, S(πG)=πS(G):=πs:s∈S(G).S(π G)=π S(G):=\π s:s∈ S(G)\. Definition 2 (Regular and stable graph instances). For each branch, let Di:=G∈:fi is not locally Lipschitz at G,D:=⋃i=1MDi.D_i:=\G :f_i is not locally Lipschitz at G\, D:= _i=1^MD_i. A graph G0G_0 has a stable collapse pattern if some relatively open neighborhood V of G0G_0 admits a partition 1,…,M=C1∪˙⋯∪˙Cr\1,…,M\=C_1 ∪·s ∪C_r such that, for every G∈VG∈ V, fi(G)=fj(G)⟺i,j∈Cα for some α.f_i(G)=f_j(G) i,j∈ C_α for some α. Let U be the set of graphs without a stable collapse pattern and define ⋆:=∖(D∪U).G_ :=G (D∪ U). We call ⋆G_ the regular stable domain. Both D and U are relatively closed in G. No emptiness, finiteness, or measure-zero assumption is imposed on D. Remark (Equivariant branch labels). At the level of the underlying set-valued map, the fixed label permutations in Definition 1 impose no additional restriction. Indeed, given any finite representation S(G)=f1(G),…,fM(G)S(G)=\f_1(G),…,f_M(G)\ satisfying S(πG)=πS(G)S(π G)=π S(G), define f~(i,σ)(G):=σfi(σ−1G),τ~π(i,σ):=(i,πσ),σ∈n. f_(i,σ)(G):=σ f_i(σ^-1G), τ_π(i,σ):=(i,πσ), σ∈ S_n. These Mn!Mn! labels represent the same solution set and satisfy Definition 1. If D and ⋆G_ are computed from the original labels, the symmetrized representation has D~=⋃σ∈nσD,~⋆=⋂σ∈nσ⋆. D= _σ∈ S_nσ D, G_ = _σ∈ S_n _ . For the second identity, a stable partition restricts to each σ-block; conversely, local continuity and separation of the distinct solution values make matches between these blocks locally constant. Thus symmetrization may shrink the regular stable domain, and probability assumptions on that domain must be rechecked. If the original representation already satisfies Definition 1, both domains are invariant by Lemma 1 and remain unchanged. Definition 3 (Switching set and strict basins). For G∈⋆G _ , the reduced switching set is ΣG∘:=y∈n:#argmins∈S(G)‖y−s‖≥2, _G := \y _n:\# *argmin_s∈ S(G)\|y-s\|≥ 2 \, where coincident branch values are counted only once. For s∈S(G)s∈ S(G), its strict Voronoi basin is Vs(G):=y∈n:‖y−s‖<‖y−r‖ for every r∈S(G)∖s.V_s(G):=\y _n:\|y-s\|<\|y-r\| for every r∈ S(G) \s\\. Theorem 4 (Jointly globally Lipschitz equivariant multi-attractor dynamics). Let S be an equivariant finite-branch solution map in the sense of Definition 1, with regular stable domain and switching geometry given by Definitions 2 and 3. Then there exists T:×n→n,(G,y)↦TG(y),T:G×Y_n _n, (G,y) T_G(y), with the following properties. (1) Regularity and symmetry. The map T is jointly globally Lipschitz and permutation equivariant. In particular, for some LT<∞L_T<∞, (6) ‖TG(y)−TH(z)‖≤LT‖G−H‖2+‖y−z‖2,TπG(πy)=πTG(y).\|T_G(y)-T_H(z)\|≤ L_T \|G-H\|^2+\|y-z\|^2, T_π G(π y)=π T_G(y). (2) Exceptional configurations. One may choose TG(y)=yT_G(y)=y whenever G∈D∪UG∈ D∪ U, or whenever G∈⋆G _ and y∈ΣG∘y∈ _G . (3) Basin-wise convergence. Fix G∈⋆G _ and y0∉ΣG∘y_0∉ _G , and let s:=argminr∈S(G)‖y0−r‖.s:= *argmin_r∈ S(G)\|y_0-r\|. Then yt+1=TG(yt)y_t+1=T_G(y_t) converges to s, and there is a trajectory-dependent constant c(G,y0)∈(0,1)c(G,y_0)∈(0,1) such that (7) ‖yt−s‖≤(1−c(G,y0))t‖y0−s‖.\|y_t-s\|≤ (1-c(G,y_0) )^t\|y_0-s\|. No uniform positive lower bound on c(G,y0)c(G,y_0) is claimed. (4) Almost-everywhere recovery. For every G∈⋆G _ , ℒm(ΣG∘)=0L^m( _G )=0. Let μ be a Borel probability measure on nY_n with μ≪ℒmμ ^m, and draw y0∼μy_0 μ. Then the iteration yty_t converges to S(G)S(G) for μ-almost every initialization. If in addition μ(Vs(G))>0μ(V_s(G))>0 for every s∈S(G)s∈ S(G), then (8) suppℒ(y∞∣G)=S(G),suppL(y_∞ G)=S(G), where y∞y_∞ denotes the almost-sure iteration limit. A.2. Bad set and displacement functions Before the proof, we state some definitions and lemmas. Definition 4 (Joint bad set). Let Z:=×nZ:=G×Y_n carry the Euclidean product metric, and define ∘:=(G,y):G∈⋆,y∈ΣG∘. :=\(G,y):G _ ,\ y∈ _G \. The joint bad set and its complement are (9) B:=((D∪U)×n)∪∘¯Z,Ω:=Z∖B.B:= ((D∪ U)×Y_n )∪ ^\,Z, :=Z B. Here, B contains singular and unstable graph instances and the closed joint switching graph. Definition 5 (Selector and normalized displacement). Define d:Z→[0,1]d:Z→[0,1] by d(G,y):=min1,dist((G,y),B),B≠∅,1,B=∅.d(G,y):= cases \1,dist((G,y),B)\,&B≠ ,\\ 1,&B= . cases For (G,y)∈Ω(G,y)∈ , the nearest distinct solution is unique; define the branch selector P and normalized displacement H as follows: (10) P(G,y):=argmins∈S(G)‖y−s‖,H(G,y):=P(G,y)−yd(G,y).P(G,y):= *argmin_s∈ S(G)\|y-s\|, H(G,y):= P(G,y)-yd(G,y). First, let’s consider the invariance/equivariance properties of the above definitions. Lemma 1. The bad set B and its complement Ω are both permutation invariant. Proof. Fix a node permutation π. Its action on G is an isometric bijection and therefore maps relatively open neighborhoods to relatively open neighborhoods. Definition 1 gives fτπ(i)=π∘fi∘π−1.f_ _π(i)=π f_i π^-1. Consequently, fif_i is locally Lipschitz at G if and only if fτπ(i)f_ _π(i) is locally Lipschitz at πGπ G. Since τπ _π is a bijection of the labels, πD=Dπ D=D. If V and Cα\C_α\ give a stable collapse pattern at G0G_0, then πVπ V and τπ(Cα)\ _π(C_α)\ give one at πG0π G_0. Applying the same argument to π−1π^-1 yields πU=Uπ U=U. In particular, π⋆=⋆ _ =G_ , and π((D∪U)×n)=(D∪U)×n.π((D∪ U)×Y_n)=(D∪ U)×Y_n. Moreover, solution-set equivariance and distance preservation imply π(∘)=∘π( )= . Indeed, let (G,y)∈∘(G,y)∈ . Then there exist distinct s1,s2∈S(G)s_1,s_2∈ S(G) such that ‖y−s1‖=‖y−s2‖=dist(y,S(G)).\|y-s_1\|=\|y-s_2\|=dist(y,S(G)). Since S(πG)=πS(G)S(π G)=π S(G) and the permutation action is an isometry, πs1π s_1 and πs2π s_2 are distinct nearest elements of S(πG)S(π G) to πyπ y. Since πG∈⋆π G _ , it follows that π(G,y)=(πG,πy)∈∘.π(G,y)=(π G,π y)∈ . Thus π(∘)⊆∘.π( ) . Applying the same argument to π−1π^-1 gives the reverse inclusion and therefore π(∘)=∘.π( )= . Since the permutation mapping is a homeomorphism, π(∘¯Z)=π(∘)¯Z=∘¯Z.π ( ^\,Z )= π( )^\,Z= ^\,Z. Using B=((D∪U)×n)∪∘¯ZB=((D∪ U)×Y_n)∪ ^\,Z now gives πB=Bπ B=B. In addition, πΩ=π(Z∖B)=πZ∖πB=Z∖B=Ω.π =π(Z B)=π Z π B=Z B= . which shows the invariance of Ω and finishes the proof. ∎ Lemma 2. Function d is permutation invariant, and P and H are permutation equivariant. Proof. Since the permutation action is an isometry and πB=Bπ B=B, d(πG,πy)=d(G,y).d(π G,π y)=d(G,y). Let s⋆=P(G,y)s_ =P(G,y). By solution-set equivariance, S(πG)=πS(G).S(π G)=π S(G). Thus every s~∈S(πG) s∈ S(π G) has the form s~=πs s=π s for some s∈S(G)s∈ S(G). Since permutations preserve distances, ‖πy−s~‖=‖y−s‖≥‖y−s⋆‖=‖πy−πs⋆‖.\|π y- s\|=\|y-s\|≥\|y-s_ \|=\|π y-π s_ \|. Hence πs⋆π s_ is a nearest solution to πyπ y. Its uniqueness on Ω gives P(πG,πy)=πP(G,y).P(π G,π y)=π P(G,y). Finally, by Definition 5, H(πG,πy)=P(πG,πy)−πyd(πG,πy)=πP(G,y)−yd(G,y)=πH(G,y),H(π G,π y)= P(π G,π y)-π yd(π G,π y)=π P(G,y)-yd(G,y)=π H(G,y), which finishes the proof. ∎ To further establish the properties of displacement H, we state and use the following lemma. Lemma 3 ([38, Theorem A.4]). Let X be an arbitrary subset of a Euclidean space, equipped with the induced metric, and let h:X→ℝh:X be locally Lipschitz intrinsically on X. Then there exists a function ε:X→(0,1) :X→(0,1) such that both ε and εh h are globally Lipschitz on X. Here intrinsic local Lipschitzness means Lipschitzness on a relative neighborhood X∩B(x,r)X∩ B(x,r) of each x∈Xx∈ X. The cited theorem requires no openness, closedness, boundedness, or local compactness of X, so it applies to X=ΩX= . We first extend this result to the vector-valued and equivariant form required below. Consider a general group setting where Γ is a finite group with orthogonal representations: ρX:Γ→O(p),ρY:Γ→O(m). _X: (p), _Y: (m). For brevity, given γ∈Γγ∈ , x∈ℝpx ^p, and v∈ℝmv ^m, we write: γx:=ρX(γ)x,γv:=ρY(γ)v.γ x:= _X(γ)x, γ v:= _Y(γ)v. The previously defined permutation π is an example of such an action; it applies to both G and y, acting as an orthogonal transformation on each domain. Lemma 4 (Equivariant bounded flattening). Let X⊆ℝpX ^p be Γ -invariant. Suppose that H:X→ℝmH:X ^m is locally Lipschitz and Γ -equivariant. Then there exists a Γ -invariant function β:X→(0,1)β:X→(0,1) such that both β and βHβ H are bounded and globally Lipschitz on X. Moreover, βHβ H is Γ -equivariant. Proof. Write H=(H1,…,Hm).H=(H_1,…,H_m). For every coordinate HjH_j, apply Lemma 3 to obtain a function εj:X→(0,1) _j:X→(0,1) such that εj _j and uj:=εjHju_j:= _jH_j are globally Lipschitz. Define aj(x):=εj(x)1+|uj(x)|.a_j(x):= _j(x)1+|u_j(x)|. The scalar map r↦(1+|r|)−1r (1+|r|)^-1 is bounded and globally Lipschitz. Since εj _j is also bounded and globally Lipschitz, aja_j is bounded and globally Lipschitz, and 0<aj<1.0<a_j<1. Moreover, ajHj=uj1+|uj|.a_jH_j= u_j1+|u_j|. Because r↦r/(1+|r|)r r/(1+|r|) is bounded and globally Lipschitz, ajHja_jH_j is bounded and globally Lipschitz as well. Set a0(x):=∏j=1maj(x).a_0(x):= _j=1^ma_j(x). A finite product of bounded globally Lipschitz functions is bounded and globally Lipschitz. For every coordinate j, a0Hj=(ajHj)∏k≠jak,a_0H_j=(a_jH_j) _k≠ ja_k, which is again bounded and globally Lipschitz. Hence a0Ha_0H is bounded and globally Lipschitz. Finally, define β(x):=∏γ∈Γa0(γx).β(x):= _γ∈ a_0(γ x). Since X is Γ -invariant, each a0∘γa_0 γ is well defined on X. Since the representation ρX _X is orthogonal, Lip(a0∘γ)≤Lip(a0),‖a0∘γ‖∞≤‖a0‖∞.Lip(a_0 γ) (a_0), \|a_0 γ\|_∞≤\|a_0\|_∞. The finiteness of Γ therefore implies that β is bounded and globally Lipschitz. For every γ0∈Γ _0∈ , β(γ0x)=∏γ∈Γa0(γγ0x)=∏γ∈Γa0(γx)=β(x),β( _0x)= _γ∈ a_0(γ _0x)= _γ∈ a_0(γ x)=β(x), where the second equality follows because right multiplication by γ0 _0 permutes the elements of Γ . Thus β is Γ -invariant. Furthermore, β(x)H(x)=a0(x)H(x)∏γ∈Γγ≠ea0(γx),β(x)H(x)=a_0(x)H(x) _ subarraycγ∈ \\ γ≠ e subarraya_0(γ x), which is a finite product of bounded globally Lipschitz factors. Thus βHβ H is bounded and globally Lipschitz. Moreover, the equivariance of βHβ H follows from the invariance of β and equivariance of H. This proves the lemma. ∎ A.3. Proof of the main result Proof of Theorem 4. We construct T and verify the four asserted properties. Step 1: The set B is closed and permutation-invariant. For each branch, the locus on which it is locally Lipschitz is relatively open. Therefore D is relatively closed in G. The set of points having a stable collapse pattern is also relatively open by definition, so U is relatively closed. It follows that (D∪U)×n(D∪ U)×Y_n is relatively closed in Z. The set ∘¯Z ^\,Z is closed by definition. Hence B is closed in Z. In addition, Lemma 1 provides the permutation-invariance of the bad set B and its complement. Step 2: The nearest-solution selector is locally Lipschitz on Z∖BZ B. Fix (G0,y0)∈Ω(G_0,y_0)∈ . Then G0∉D∪UG_0∉ D∪ U, the collapse pattern is stable near G0G_0, and y0y_0 has a unique nearest distinct branch value. Let C⋆C_ be the stable cluster realizing this nearest value, and choose any representative i⋆∈C⋆i_ ∈ C_ . Branches in C⋆C_ agree identically throughout some neighborhood of G0G_0. If C⋆C_ is the only stable cluster, then all branch values agree throughout that neighborhood, and hence P(G,y)=fi⋆(G)P(G,y)=f_i_ (G) there; the local Lipschitz conclusion follows immediately. Otherwise, choose one representative iαi_α from every other stable cluster. Since the nearest branch value at (G0,y0)(G_0,y_0) is unique, the margin δ:=minα≠⋆(∥y0−fiα(G0)∥−∥y0−fi⋆(G0)∥)δ:= _α≠ (\|y_0-f_i_α(G_0)\|-\|y_0-f_i_ (G_0)\| ) is strictly positive. All regular branches are continuous near G0G_0. After shrinking the neighborhood of (G0,y0)(G_0,y_0) if necessary, C⋆C_ remains the unique nearest cluster throughout that neighborhood. Hence P(G,y)=fi⋆(G)P(G,y)=f_i_ (G) there. Because G0∉DG_0∉ D, the representative branch fi⋆f_i_ is locally Lipschitz near G0G_0. Thus P is jointly locally Lipschitz near (G0,y0)(G_0,y_0). Since (G0,y0)(G_0,y_0) was arbitrary, P:Ω→nP: _n is jointly locally Lipschitz. This is the same local mechanism used in [30, Lemma A.12], but no boundedness or compactness of the graph domain is required here. Step 3: Construction of a damping factor that removes the singularities. Since B is relatively closed in Z, for every z∈Z∖Bz∈ Z B we have dist(z,B)>0dist(z,B)>0; otherwise z∈B¯Z=Bz∈ B^\,Z=B. Hence it holds that d(z)>0for z∈Ω,d(z)=0for z∈B.d(z)>0 z∈ , d(z)=0 z∈ B. Moreover, d is bounded, 11-Lipschitz, and permutation-invariant. Since d is 11-Lipschitz and strictly positive on Ω , it is locally bounded away from zero there. Hence d−1d^-1 is locally Lipschitz on Ω . Together with the local Lipschitzness of P, this implies that H=(P−y)/dH=(P-y)/d is locally Lipschitz. In addition, thanks to Lemma 2, H is permutation equivariant. Apply Lemma 4 on Ω and permutation action group. There exists a permutation-invariant function β:Ω→(0,1)β: →(0,1) such that A:=βHA:=β H is bounded and globally Lipschitz on Ω . Define q(G,y):=d(G,y)A(G,y),(G,y)∈Ω,0,(G,y)∈B.q(G,y):= casesd(G,y)A(G,y),&(G,y)∈ ,\\ 0,&(G,y)∈ B. cases For (G,y)∈Ω(G,y)∈ , q(G,y)=β(G,y)(P(G,y)−y).q(G,y)=β(G,y) (P(G,y)-y ). Finally, define TG(y):=y+q(G,y).T_G(y):=y+q(G,y). Equivalently, (11) TG(y)=(1−β(G,y))y+β(G,y)PG(y),(G,y)∉B,y,(G,y)∈B.T_G(y)= cases(1-β(G,y))y+β(G,y)P_G(y),&(G,y)∉ B,\\[2.84526pt] y,&(G,y)∈ B. cases The freezing property on B is immediate from this definition. Step 4: Joint global Lipschitz continuity of T. Let ‖A‖∞≤MA,Lip(A)≤LA.\|A\|_∞≤ M_A, (A)≤ L_A. For any z,z′∈Ωz,z ∈ , ‖q(z)−q(z′)‖ \|q(z)-q(z )\| =‖d(z)A(z)−d(z′)A(z′)‖ =\|d(z)A(z)-d(z )A(z )\| ≤|d(z)−d(z′)|‖A(z)‖+d(z′)‖A(z)−A(z′)‖ ≤|d(z)-d(z )|\,\|A(z)\|+d(z )\|A(z)-A(z )\| ≤(MA+LA)‖z−z′‖, ≤(M_A+L_A)\|z-z \|, where we used 0≤d≤10≤ d≤ 1 (This follows from Definition 5) and the 11-Lipschitz property of d. If z∈Ωz∈ and z′∈Bz ∈ B, then ‖q(z)−q(z′)‖=|q(z)|=d(z)‖A(z)‖≤MAdist(z,B)≤MA‖z−z′‖.\|q(z)-q(z )\|=\|q(z)\|=d(z)\|A(z)\|≤ M_Adist(z,B)≤ M_A\|z-z \|. If both points lie in B, the difference is zero. Thus q is globally Lipschitz on all of Z. Since the coordinate projection (G,y)↦y(G,y) y is 11-Lipschitz and TG(y)=y+q(G,y),T_G(y)=y+q(G,y), T is jointly globally Lipschitz. In particular, one may take LT≤1+MA+LA.L_T≤ 1+M_A+L_A. Step 5: Permutation equivariance. The functions d and β are invariant, whereas H is equivariant. Therefore A(πG,πy)=πA(G,y)A(π G,π y)=π A(G,y) on Ω . Since B is invariant, it follows that q(πG,πy)=πq(G,y)q(π G,π y)=π q(G,y) on all of Z. Consequently, TπG(πy)=πy+q(πG,πy)=πy+πq(G,y)=πTG(y).T_π G(π y)=π y+q(π G,π y)=π y+π q(G,y)=π T_G(y). Step 6: On a regular stable graph, taking the closure does not enlarge the switching section. Fix G∈⋆G _ . We claim that (12) y:(G,y)∈B=ΣG∘.\y:(G,y)∈ B\= _G . The inclusion “⊇ ” follows immediately from the definition of the closure. For the reverse inclusion, it suffices to show that any point in the G-section of ∘¯Z ^\,Z belongs to ΣG∘ _G . Suppose that (Gk,yk)∈∘,(Gk,yk)→(G,y).(G_k,y_k)∈ , (G_k,y_k)→(G,y). For every k, choose two labels whose distinct branch values attain the common minimum distance from yky_k. Because there are only finitely many pairs of labels, after passing to a subsequence we may assume that these labels are fixed, say i≠ji≠ j. Since G is a stable-collapse point, if i and j belonged to the same stable cluster, then the two branches would agree identically throughout a neighborhood of G. For all sufficiently large k, this would contradict the fact that they represent distinct branch values at GkG_k. Therefore i and j belong to different stable clusters, and in particular fi(G)≠fj(G).f_i(G)≠ f_j(G). Because G∉DG∉ D, every branch is continuous at G. Passing to the limit in the equal-minimum relations gives ‖y−fi(G)‖=‖y−fj(G)‖=minℓ‖y−fℓ(G)‖.\|y-f_i(G)\|=\|y-f_j(G)\|= _ \|y-f_ (G)\|. Thus y∈ΣG∘y∈ _G , proving the reverse inclusion. Step 7: Invariance of the strict Voronoi segment. Fix G∈⋆,y0∉ΣG∘,G _ , y_0∉ _G , and let s=PG(y0).s=P_G(y_0). For any other distinct solution r∈S(G)∖sr∈ S(G) \s\, ‖y0−s‖<‖y0−r‖.\|y_0-s\|<\|y_0-r\|. For a∈[0,1]a∈[0,1], define za=(1−a)y0+as.z_a=(1-a)y_0+as. Then ‖za−s‖=(1−a)‖y0−s‖.\|z_a-s\|=(1-a)\|y_0-s\|. On the other hand, the reverse triangle inequality yields ‖za−r‖ \|z_a-r\| ≥‖y0−r‖−‖za−y0‖ ≥\|y_0-r\|-\|z_a-y_0\| =‖y0−r‖−a‖y0−s‖ =\|y_0-r\|-a\|y_0-s\| >(1−a)‖y0−s‖ >(1-a)\|y_0-s\| =‖za−s‖. =\|z_a-s\|. Therefore PG(za)=sandza∉ΣG∘P_G(z_a)=s z_a∉ _G for every a∈[0,1]a∈[0,1]. Hence PG(z)=sfor every z∈[y0,s],P_G(z)=s every z∈[y_0,s], and Step 6 implies G×[y0,s]⊆Ω.\G\×[y_0,s] . Step 8: Convergence. The set G×[y0,s]\G\×[y_0,s] is a compact subset of Ω . Since β is continuous and strictly positive on Ω , c(G,y0):=minz∈[y0,s]β(G,z)>0.c(G,y_0):= _z∈[y_0,s]β(G,z)>0. Furthermore, 0<β<10<β<1. If yt∈[y0,s]y_t∈[y_0,s], then Step 7 gives PG(yt)=sP_G(y_t)=s, and therefore yt+1=(1−β(G,yt))yt+β(G,yt)s∈[yt,s]⊆[y0,s].y_t+1=(1-β(G,y_t))y_t+β(G,y_t)s∈[y_t,s] [y_0,s]. By induction, every iterate remains on [y0,s][y_0,s]. Moreover, (13) yt+1−s=(1−β(G,yt))(yt−s).y_t+1-s= (1-β(G,y_t) )(y_t-s). It follows that ‖yt+1−s‖≤(1−c(G,y0))‖yt−s‖.\|y_t+1-s\|≤ (1-c(G,y_0) )\|y_t-s\|. Iterating this estimate gives ‖yt−s‖≤(1−c(G,y0))t‖y0−s‖.\|y_t-s\|≤ (1-c(G,y_0) )^t\|y_0-s\|. Hence yt→sy_t→ s. Step 9: Almost-everywhere convergence and full support. Fix G∈⋆G _ . For any two distinct solutions s,r∈S(G)s,r∈ S(G), the equidistance set y:‖y−s‖=‖y−r‖\y:\|y-s\|=\|y-r\|\ is an affine hyperplane, because the defining equality is equivalent to 2⟨y,r−s⟩=‖r‖2−‖s‖2.2 y,r-s =\|r\|^2-\|s\|^2. Therefore ΣG∘⊆⋃s,r∈S(G)s≠ry:‖y−s‖=‖y−r‖ _G _ subarraycs,r∈ S(G)\\ s≠ r subarray\y:\|y-s\|=\|y-r\|\ is contained in a finite union of Lebesgue-null affine hyperplanes. Thus ℒm(ΣG∘)=0.L^m( _G )=0. If μ≪ℒmμ ^m, then μ(ΣG∘)=0.μ( _G )=0. Step 8 therefore shows that the iteration converges to a point of S(G)S(G) for μ-almost every initialization. Finally, if every strict Voronoi basin has positive μ-mass, then for every s∈S(G)s∈ S(G), Pr(y∞=s∣G)=μ(Vs(G))>0. (y_∞=s G)=μ(V_s(G))>0. Since S(G)S(G) is finite and y∞∈S(G)y_∞∈ S(G) almost surely, suppℒ(y∞∣G)=S(G).suppL(y_∞ G)=S(G). This proves all the conclusions. ∎ Appendix B Family-wide asymptotic realization by a recurrent MP-GNN We now realize the ideal dynamics of Appendix A using one message-passing GNN block that is reused at every recurrent step: (14) yt+1=GNN(G,yt).y_t+1=GNN(G,y_t). We suppress the parameter subscript θ in this appendix. The random initialization supplies both the initial node addresses and the choice of solution basin. At later steps, the block receives only the graph and the current state. Throughout this appendix, write dy:=q≥2d_y:=q≥ 2 and m=ndym=nd_y, using the graph and output spaces of Appendix A. B.1. Architecture and main result The recurrent block. For one call, write y+=GNN(G,y)y^+=GNN(G,y). Choose a finite number of hidden layers L, finite hidden widths, and a continuous encoder enc:ℝdv×ℝdy⟶ℝd0.enc:R^d_v×R^d_y ^d_0. The block has the form (15) hi(0) h_i^(0) =enc(Xi,yi), =enc(X_i,y_i), mi(ℓ) m_i^( ) =∑j=1nMℓ(hi(ℓ),hj(ℓ),Aji), = _j=1^nM_ \! (h_i^( ),h_j^( ),A_ji ), hi(ℓ+1) h_i^( +1) =Uℓ(hi(ℓ),mi(ℓ)),0≤ℓ<L, =U_ \! (h_i^( ),m_i^( ) ), 0≤ <L, g(L) g^(L) =∑j=1nRL(hj(L)), = _j=1^nR_L\! (h_j^(L) ), yi+ y_i^+ =O(hi(L),g(L)). =O\! (h_i^(L),g^(L) ). Here the aggregation runs over all n senders, including j=ij=i, since the edge tensor A assigns a feature vector to every ordered pair. The component maps, Mℓ,Uℓ=0L−1,RL,O\M_ ,U_ \_ =0^L-1,R_L,O, are arbitrary continuous maps between finite-dimensional Euclidean spaces and are shared across nodes. The only global operation is the final invariant sum: no global aggregate is formed or broadcast in the hidden layers. The entire L-layer block, including its encoder and readout, is denoted as GNN(G,⋅)GNN(G,·) and is reused with the same parameters at every recurrent time. Hidden representations are recomputed within each call; only yty_t is carried to the next call. Remark (Sparse graphs). For an edge set E(G)⊆[n]×[n]E(G) [n]×[n], encode Aji=(eji,aji)A_ji=(e_ji,a_ji), where eji=(j,i)∈E(G)e_ji=1_\(j,i)∈ E(G)\ and aji=0a_ji=0 whenever eji=0e_ji=0. The class of such tensors is closed under relabeling, so Appendix A applies to any permutation-invariant domain of these graphs. The block constructed in Lemma 5 has M0=ci⊗cj⊗AjiM_0=c_i c_j A_ji, which vanishes on every absent pair, even outside the compact realization family. Hence its all-pairs message sum agrees exactly with the neighbor sum over (j,i)∈E(G)(j,i)∈ E(G). On the compact realization family, ℛedgeR_ edge records (1,aji)(1,a_ji) in present-edge slots and zero in absent-edge slots. Consequently, the block in Theorem 5 can use neighbor-only hidden message passing followed by the readout in (15). It’s straightforward that reindexing the message and readout sums gives the equivariance: (16) GNN(πG,πy)=πGNN(G,y).GNN(π G,π y)= (G,y). Facts about the ideal update. Fix the operator T constructed in Theorem 4. Let Ω and P be as in Definitions 4 and 5. The proof in Appendix A gives (17) TG(y)=(1−β(G,y))y+β(G,y)P(G,y),0<β(G,y)<1,(G,y)∈Ω,T_G(y)=(1-β(G,y))y+β(G,y)P(G,y), 0<β(G,y)<1, (G,y)∈ , where β is continuous. The selector is continuous on Ω and is constant along the segment from y to P(G,y)P(G,y). We will also use the bound (18) ‖TG(y)−y‖≤dist(y,S(G))(G∈⋆,y∈n).\|T_G(y)-y\| (y,S(G)) (G _ ,\ y _n). Indeed, off the switching set this follows from (17); on the switching set TG(y)=yT_G(y)=y. Theorem 5 (Family-wide asymptotic realization from random initialization). Assume dy≥2d_y≥ 2. Let (G,Y0)(G,Y_0) have a tight Borel probability law ℙP on ×nG×Y_n, and let μG _G be a regular conditional distribution of Y0Y_0 given G. Such a Borel kernel exists by disintegrating the joint law on the ambient Euclidean product and restricting the graph variable to G. Assume: (1) G∈⋆G _ almost surely; (2) conditionally on G, the law μG _G of Y0Y_0 satisfies μG≪ℒm _G ^m almost surely; (3) for almost every G and every s∈S(G)s∈ S(G), μG(Vs(G))>0. _G(V_s(G))>0. All solution-dependent probability events are understood on the full-measure set G∈⋆G _ . Fix a limit tolerance ε>0 >0, an operator tolerance ρ>0ρ>0, and a probability loss δ∈(0,1)δ∈(0,1). Then there exist a compact permutation-invariant event Eε,ρ,δ⊆×nE_ ,ρ,δ ×Y_n with (19) ℙ(Eε,ρ,δ)≥1−δP(E_ ,ρ,δ)≥ 1-δ and one deterministic, autonomous, weight-tied MP-GNN block such that, for every (G,y0)∈Eε,ρ,δ(G,y_0)∈ E_ ,ρ,δ, the trajectory y^0=y0,y^t+1=GNN(G,y^t),t≥0, y_0=y_0, y_t+1=GNN(G, y_t), t≥ 0, converges to a limit y^∞ y_∞ satisfying (20) dist(y^∞,S(G)) ( y_∞,S(G)) <ε, < , supt≥0‖GNN(G,y^t)−TG(y^t)‖ _t≥ 0 \|GNN(G, y_t)-T_G( y_t) \| <ρ. <ρ. The block receives only (G,y^t)(G, y_t) at every recurrent time. The graph projection ε,ρ,δ:=projEε,ρ,δK_ ,ρ,δ:=proj_GE_ ,ρ,δ is compact, and the same fixed collection of continuous component maps works for every graph in this projection on its corresponding retained initialization fiber. The block is permutation equivariant as in (16). In particular, ℙ(y^t converges and dist(y^∞,S(G))<ε)≥1−δ.P\! ( y_t converges and dist( y_∞,S(G))< )≥ 1-δ. Moreover, the good event may be chosen so that every labeled branch has positive joint capture probability: (21) ℙ(Y0∈Vfa(G)(G),y^t converges to y^∞,∥y^∞−fa(G)∥<ε)>0,a=1,…,M.P\! (Y_0∈ V_f_a(G)(G), y_t converges to y_∞, \| y_∞-f_a(G)\|< )>0, a=1,…,M. Coincident labels refer to the same distinct target. The learned limit is required to be ε -close to the solution set; it need not belong exactly to S(G)S(G). The probability law need not be permutation invariant. Equivariance is a property of the block, and the retained event is made invariant by its construction. The branch-capture guarantee is joint over graphs and initializations; it does not assert positive retained mass in every basin of every individual graph. B.2. Exact realization using the current node states In this subsection, we establish conditions under which a continuous equivariant operator can be exactly realized by a message-passing GNN in our setting. We begin with the relevant symmetry-breaking condition. Without suitable symmetry-breaking features, standard message-passing GNNs may fail to distinguish graph inputs that the target operator distinguishes, limiting their ability to represent arbitrary continuous equivariant operators (see [66]). We therefore introduce a condition under which the current node states provide distinct node addresses, and prove that these addresses enable exact realization on a compact family. This result provides the basis for constructing a message-passing GNN that follows the ideal dynamics introduced in Appendix A. Definition 6 (Robust state-based node symmetry breaking). Let K⊆×nK ×Y_n be compact and permutation invariant. We say that K admits robust state-based node symmetry breaking if there exist finitely many compact address regions 1,…,J⊆ℝdyA_1,…,A_J ^d_y and η>0η>0 such that dist(j,k)≥η(j≠k),dist(A_j,A_k)≥η (j≠ k), and, for every (G,y)∈K(G,y)∈ K, the n node arguments y1,…,yny_1,…,y_n belong to n different address regions. Because there are finitely many compact address regions, pairwise disjointness already guarantees a positive separation margin; η simply denotes a common margin. The regions 1,…,JA_1,…,A_J and their separation margin η are fixed for the whole family K, while the assignment of nodes to regions may vary between inputs. Since the node states y1,…,yny_1,…,y_n belong to different address regions, each yiy_i serves as a unique node identifier within the current input (G,y)(G,y). For example, take dy=2d_y=2 and a fixed two-node graph G0G_0 that is unchanged by swapping its nodes. Write uθ=(cosθ,sinθ)u_θ=( θ, θ) and define two opposite arcs: 1=uθ:θ∈[−π/4,π/4],2=−1.A_1=\u_θ:θ∈[-π/4,π/4]\, _2=-A_1. These compact regions lie on the right and left sides of the unit circle, respectively, with dist(1,2)=2dist(A_1,A_2)= 2. The compact permutation-invariant family K=(G0,(u,−u)),(G0,(−u,u)):u∈1K= \(G_0,(u,-u)),\ (G_0,(-u,u)):u _1 \ therefore satisfies Definition 6 (i.e., admits robust state-based node symmetry breaking) with J=2J=2 and η=2η= 2: one node lies in each region, although their assignments can be swapped. In contrast, the following full-circle family does not satisfy Definition 6 Kcircle=(G0,(uθ,−uθ)):θ∈[0,2π].K_ circle= \ (G_0,(u_θ,-u_θ) ):θ∈[0,2π] \. The two node states are still always distance 22 apart, but each now ranges over the entire circle. There is no finite collection of compact regions that simultaneously (i) are positively separated, (i) cover the entire unit circle, and (i) place uθu_θ and −uθ-u_θ in different regions for every θ. Indeed, because the circle is connected, conditions (i) and (i) force it to lie entirely in one region, contradicting (i). In the remainder of this subsection, we show how continuous equivariant operators can be exactly realized by message-passing GNNs under robust state-based node symmetry breaking. Subsections B.3 and B.4 then establish this condition for the compact families used in our construction. Lemma 5 (Exact MP-GNN realization on a compact family). Let K⊆×nK ×Y_n be compact, permutation invariant, and admit robust state-based node symmetry breaking. Fix integers u,v≥0u,v≥ 0 with u+v≥1u+v≥ 1 and, when v>0v>0, positive integers p1,…,pvp_1,…,p_v. For ℓ=1,…,u =1,…,u, let φℓ:K→ℝ _ :K be continuous and permutation invariant. For b=1,…,vb=1,…,v, let Fb:K→ℝn×pbF_b:K ^n× p_b be continuous and permutation equivariant. Then there exists a MP-GNN, with output width u+∑b=1vpbu+ _b=1^vp_b and the structural form (15) realizes all these maps simultaneously and exactly on K: its output row at node i is τ(G,y,i):=(φ1(G,y),…,φu(G,y),[F1(G,y)]i,…,[Fv(G,y)]i).τ(G,y,i):= ( _1(G,y),…, _u(G,y),[F_1(G,y)]_i,…,[F_v(G,y)]_i ). Empty output blocks are omitted. The component maps and all widths are chosen once for the whole compact family. Proof. We construct one message-passing layer that records the graph and state in node-code order, and then define a continuous readout on that record. Step 1: Construct the message-passing representation. Fix one collection of compact address regions 1,…,JA_1,…,A_J satisfying Definition 6 simultaneously for all (G,y)∈K(G,y)∈ K. These regions remain fixed throughout the construction; only the region containing each node state may vary with the input. The map that equals the jjth standard basis vector ej∈ℝJe_j ^J on jA_j is continuous on their closed union. The Tietze extension theorem [44, Theorem 35.1(b)], applied coordinatewise, gives a continuous map c:ℝdy⟶ℝJ,c(w)=ej(w∈j).c:R^d_y ^J, c(w)=e_j (w _j). For (G,y)∈K(G,y)∈ K, the codes ci:=c(yi)c_i:=c(y_i) are distinct standard basis vectors. Set d0=dv+dy+Jd_0=d_v+d_y+J and define enc(Xi,yi):=(Xi,yi,c(yi)).enc(X_i,y_i):=(X_i,y_i,c(y_i)). For the single local layer, take M0(hi(0),hj(0),Aji) M_0\! (h_i^(0),h_j^(0),A_ji ) :=ci⊗cj⊗Aji, :=c_i c_j A_ji, U0(hi(0),mi(0)) U_0\! (h_i^(0),m_i^(0) ) :=(ci,ci⊗(1,Xi,yi),mi(0)). := (c_i,\ c_i (1,X_i,y_i),\ m_i^(0) ). Thus hi(1)=(ci,Ni,Ei)h_i^(1)=(c_i,N_i,E_i), where Ni:=ci⊗(1,Xi,yi),Ei:=∑j=1nci⊗cj⊗Aji.N_i:=c_i (1,X_i,y_i), E_i:= _j=1^nc_i c_j A_ji. Set R1(ci,Ni,Ei):=(Ni,Ei)R_1(c_i,N_i,E_i):=(N_i,E_i). The final global sum is g(1)=(ℛnode,ℛedge),ℛnode:=∑ici⊗(1,Xi,yi),ℛedge:=∑i,jci⊗cj⊗Aji.g^(1)=(R_ node,R_ edge), _ node:= _ic_i (1,X_i,y_i), _ edge:= _i,jc_i c_j A_ji. All tensors are vectorized before concatenation. The displayed component maps are continuous and shared across nodes. Step 2: Recover the graph up to relabeling. For a standard basis vector eke_k, the tensor ek⊗ze_k z places z in its kkth block and zero in the remaining blocks. Hence ℛnodeR_ node records (1,Xi,yi)(1,X_i,y_i) in the slot indexed by cic_i, and ℛedgeR_ edge records AjiA_ji in the slot indexed by the ordered pair (ci,cj)(c_i,c_j). Note that, the roles of the constant-one coordinates are therefore as follows: The 1 in (1,Xi,yi)(1,X_i,y_i) is a node-occupancy indicator. Without it, a node satisfying Xi=0X_i=0 and yi=0y_i=0 would produce a zero block, indistinguishable from an unused address slot. This ensures that the global tensors reconstruct nodes and edges even when their features are zero; they do not carry solution information or participate in symmetry breaking. Suppose (G,y)(G,y) and (G′,y′)(G ,y ) produce the same two global tensors. Both inputs have the same number n of nodes, each occupying a distinct address slot marked by the constant-one coordinate. Matching these occupied slots therefore defines a unique permutation π∈nπ∈ S_n such that cπ(i)′=ci,i∈[n].c _π(i)=c_i, i∈[n]. Equality of the corresponding blocks then implies Xπ(i)′=Xi,yπ(i)′=yi,Aπ(j),π(i)′=Aji.X _π(i)=X_i, y _π(i)=y_i, A _π(j),π(i)=A_ji. Thus (G′,y′)=(πG,πy)(G ,y )=(π G,π y). If the receiver codes also agree, ci′=cic _i =c_i, their distinctness gives i′=π(i)i =π(i). Therefore Γ(G,y,i):=(ci,ℛnode,ℛedge) (G,y,i):=(c_i,R_ node,R_ edge) determines the input and its receiving node up to simultaneous relabeling. Step 3: The target is determined by this representation. If Γ(G,y,i)=Γ(G′,y′,i′) (G,y,i)= (G ,y ,i ), Step 2 supplies a permutation matching both inputs and receivers. Invariance of each φℓ _ and equivariance of each FbF_b imply τ(G,y,i)=τ(G′,y′,i′).τ(G,y,i)=τ(G ,y ,i ). Consequently, the prescription D0(Γ(G,y,i)):=τ(G,y,i)D_0 ( (G,y,i) ):=τ(G,y,i) defines a single-valued map on the reconstruction image ℐ:=Γ(K×[n])I:= (K×[n]). Step 4: Construct a continuous output map. Both Γ and τ are continuous on the compact space K×[n]K×[n], where [n][n] has the discrete topology. The induced map D0D_0 is continuous. For completeness, if wk→w_k→ w in ℐI, choose preimages zk∈K×[n]z_k∈ K×[n]. Every subsequence of these preimages has a further convergent subsequence, say zkj→z_k_j→ z. Continuity gives Γ(z)=w (z)=w and D0(wkj)=τ(zkj)⟶τ(z)=D0(w).D_0(w_k_j)=τ(z_k_j) τ(z)=D_0(w). If D0(wk)D_0(w_k) did not converge to D0(w)D_0(w), a subsequence staying a positive distance away would contradict this argument. Thus D0D_0 is continuous. This is the compact-domain quotient argument in explicit form. The image ℐI is compact and hence closed in its Euclidean ambient space. Coordinatewise Tietze extension gives a continuous ambient map D~ D agreeing with D0D_0 on ℐI. Define O((ci,Ni,Ei),g(1)):=D~(ci,ℛnode,ℛedge).O ((c_i,N_i,E_i),g^(1) ):= D(c_i,R_ node,R_ edge). On K, this output equals τ(G,y,i)τ(G,y,i), as required. ∎ Lemma 5 shows that, once robust node-state symmetry breaking is available, any prescribed finite collection of continuous permutation-invariant scalar maps and continuous permutation-equivariant node maps can be realized simultaneously and exactly on K. This proof differs from the Stone–Weierstrass density arguments commonly used in GNN universality theory, e.g., [33, 4, 21]. It constructs a finite-dimensional representation (ci,ℛnode,ℛedge)(c_i,R_ node,R_ edge) that reconstructs (G,y,i)(G,y,i) up to simultaneous node relabeling. The target therefore factors continuously through this representation, and the Tietze extension theorem provides an ambient continuous readout. Because the present MP-GNN class allows arbitrary continuous component maps, this construction yields exact realization. By contrast, Stone–Weierstrass uses a function algebra that contains the constants and separates the relevant points or equivalence classes to establish density in in the relevant continuous-function space, yielding uniform approximation rather than exact equality. B.3. Node separation and a common collection of boxes To apply Lemma 5, we need the node states to provide robust addresses (symmetry breaking). The next two lemmas establish this condition: the first derives node separation from random initialization, and the second places the separated states in common address boxes, allowing arbitrarily small probability losses. Lemma 6 (Finite-time node separation for the ideal dynamics). Assume dy≥2d_y≥ 2, G∈⋆G _ almost surely, and μG≪ℒm _G ^m almost surely. Put Yt:=TGt(Y0)Y_t:=T_G^t(Y_0). Then (22) ℙ((Yt)i≠(Yt)j for all t∈ℕ0 and i≠j|G)=1for almost every G.P\! ((Y_t)_i≠(Y_t)_j for all t _0 and i≠ j\, |\,G )=1 almost every G. In particular, this event has joint probability one. Let E⊆⋆×nE _ ×Y_n be compact and permutation invariant, let H∈ℕ0H _0, and let γ>0γ>0. There exist η>0η>0 and a compact invariant E′⊆E E such that ℙ(E∖E′) (E E ) <γ, <γ, ‖[TGt(y0)]i−[TGt(y0)]j‖ \|[T_G^t(y_0)]_i-[T_G^t(y_0)]_j\| ≥η((G,y0)∈E′, 0≤t≤H,i≠j). ≥η ((G,y_0)∈ E ,\ 0≤ t≤ H,\ i≠ j). If finitely many measurable subsets Ca⊆EC_a E have positive probability, E′E may also be chosen with ℙ(Ca∩E′)>0P(C_a∩ E )>0 for every a. Proof. For n=1n=1 there are no node pairs; take E′=E =E and any η>0η>0. Assume n≥2n≥ 2, and fix a graph G∈⋆G _ for which the conditional initialization law is absolutely continuous. The key observation is that a collision forces an initial node state difference onto a fixed line determined by a solution. In two or more dimensions, such an event has zero measure. For i<ji<j, the node-state-difference map Dij:n→ℝdy,Dij(y)=yi−yj,D_ij:Y_n ^d_y, D_ij(y)=y_i-y_j, is linear and surjective. For any s∈S(G)s∈ S(G), define Ns,ij:=Dij−1(spansi−sj),N_s,ij:=D_ij^-1 (span\s_i-s_j\ ), which consists of states whose node difference yi−yjy_i-y_j is a scalar multiple of the target difference si−sjs_i-s_j. The span on the right has dimension at most one. Write W=spansi−sjW=span\s_i-s_j\. Since DijD_ij is surjective, the rank–nullity theorem gives dimkerDij=m−dy D_ij=m-d_y. Its restriction to Ns,ij=Dij−1(W)N_s,ij=D_ij^-1(W) is a surjective linear map onto W with the same kernel. Applying rank–nullity to this restriction yields dimNs,ij=dimkerDij+dimW=m−dy+dimspansi−sj. N_s,ij= D_ij+ W=m-d_y+ \s_i-s_j\. The span has dimension zero if si=sjs_i=s_j and one otherwise. Thus, since dy≥2d_y≥ 2, dimNs,ij≤m−dy+1≤m−1<m. N_s,ij≤ m-d_y+1≤ m-1<m. Thus Ns,ijN_s,ij is a proper linear subspace and has zero Lebesgue measure. Define set NG:=ΣG∘⋃(⋃s∈S(G)⋃i<jNs,ij)N_G:= _G ( _s∈ S(G)\ _i<jN_s,ij ) where ΣG∘⊆n _G _n is the switching set for the fixed graph G, as defined in Definition 3. The set NGN_G collects exceptional initializations excluded to ensure a unique selected branch from f1(G),⋯,fM(G)\f_1(G),·s,f_M(G)\ and rule out node-state collisions at every finite time. The set NGN_G is Lebesgue null: the switching set is Lebesgue null by Theorem 4, and the remaining union is finite. Absolute continuity gives μG(NG)=0 _G(N_G)=0. For y0∉NGy_0∉ N_G, let s=P(G,y0)s=P(G,y_0). Basin preservation and (17) give (23) yt=s+at(y0−s),at:=∏k=0t−1(1−β(G,yk))>0,a0:=1.y_t=s+a_t(y_0-s), a_t:= _k=0^t-1(1-β(G,y_k))>0, a_0:=1. If (yt)i=(yt)j(y_t)_i=(y_t)_j at a finite time, then 0=at((y0)i−(y0)j)+(1−at)(si−sj).0=a_t ((y_0)_i-(y_0)_j )+(1-a_t)(s_i-s_j). Since at>0a_t>0, the initial row difference belongs to spansi−sjspan\s_i-s_j\, contradicting y0∉NGy_0∉ N_G. The same null set works for all finite times. This proves the conditional claim (22), and integration over G proves the joint claim. For the compact restriction, define ΔH(G,y0):=min0≤t≤Hi<j‖[TGt(y0)]i−[TGt(y0)]j‖. _H(G,y_0):= _ subarrayc0≤ t≤ H\\ i<j subarray\|[T_G^t(y_0)]_i-[T_G^t(y_0)]_j\|. Finite iterates of T are jointly continuous and equivariant, so ΔH _H is continuous and invariant. The probability-one statement implies ℙ(E∩ΔH=0)=0P(E∩\ _H=0\)=0, and hence ℙ(E∩ΔH<η)⟶0(η↓0).P(E∩\ _H<η\) 0 (η 0). Choose η>0η>0 so that this loss is below γ and, when subsets CaC_a are prescribed, below 12minaℙ(Ca) 12 _aP(C_a). Then E′:=E∩ΔH≥ηE :=E∩\ _H≥η\ is compact, invariant, and has all the required properties. ∎ Lemma 7 (Separated boxes for finitely many state maps). Let E be a compact metric space with a finite Borel measure μ. Let T be a nonempty finite index set, and let Vt=(vt,1,…,vt,n):E⟶ℝn×dy,t∈,V_t=(v_t,1,…,v_t,n):E ^n× d_y, t , be continuous. (Here i denotes the node index.) Suppose, for some η>0η>0, ‖vt,i(x)−vt,j(x)‖≥η(x∈E,t∈,i≠j).\|v_t,i(x)-v_t,j(x)\|≥η (x∈ E,\ t ,\ i≠ j). For any ω,γ>0ω,γ>0, there exist a compact E′⊆E E and finitely many pairwise positively separated compact axis-aligned boxes Q1,…,QJ⊂ℝdyQ_1,…,Q_J ^d_y, with centers b1,…,bJb_1,…,b_J, such that: (1) μ(E∖E′)<γμ(E E )<γ; (2) every vt,i(x)v_t,i(x) lies in exactly one box, and the n items of each Vt(x)V_t(x) lie in different boxes; (3) replacing every item in Vt(x)V_t(x) by its box center changes the full state by less than ω: ‖Vt(x)−b(Vt(x))‖<ω(x∈E′,t∈).\|V_t(x)-b(V_t(x))\|<ω (x∈ E ,\ t ). If finitely many measurable subsets Ca⊆EC_a E have positive measure, E′E can preserve positive measure in every CaC_a. If E carries a continuous node-permutation action and Vt(πx)=πVt(x)V_t(π x)=π V_t(x), then E′E can also be chosen permutation invariant. Proof. We construct one common grid in the node-state space ℝdyR^d_y. First, we choose its mesh small enough to separate different nodes and control the error from replacing states by cell centers. Second, we shift the grid so that the set of inputs placing any node state on a grid boundary has measure zero. Third, we remove inputs whose node states lie in a sufficiently thin neighborhood of these boundaries, losing arbitrarily little measure. The remaining node states then lie in finitely many separated compact boxes. Step 1: Choose a sufficiently small mesh. Since E is compact and there are only finitely many continuous maps vt,iv_t,i, the union of all their ranges is compact and hence bounded. Choose a grid mesh ℓ satisfying 0<ℓ<min2ωndy,ηdy.0< < \ 2ω nd_y, η d_y \. For a shift vector a∈[0,ℓ)dya∈[0, )^d_y, the grid cells are a+ℓz+[0,ℓ]dy,z∈ℤdy.a+ z+[0, ]^d_y, z ^d_y. Each cell has diameter dyℓ<η d_y <η, so it cannot contain two node states from the same Vt(x)V_t(x), whose distance is at least η. The other bound on ℓ will control the center-replacement error in Step 4. Step 2: Shift the grid to give its boundaries zero measure. Fix a coordinate s∈[dy]s∈[d_y]. We first identify the coordinate values that occur on a set of inputs of positive measure: Bs:=b∈ℝ:μx∈E:(vt,i(x))s=b>0 for some t∈,i∈[n].B_s:= \b :μ\x∈ E:(v_t,i(x))_s=b\>0 for some t ,\ i∈[n] \. This set is at most countable. Indeed, for any fixed (t,i)(t,i) and positive integer p, only finitely many values b can satisfy μx∈E:(vt,i(x))s=b≥1p.μ\x∈ E:(v_t,i(x))_s=b\≥ 1p. The corresponding level sets are disjoint, and infinitely many such sets would contradict μ(E)<∞μ(E)<∞. Every positive-measure level has measure at least 1/p1/p for some p, so there are at most countably many such levels. Taking the finite union over (t,i)(t,i) shows that BsB_s is countable. The grid boundaries in coordinate s occur at the values as+kℓa_s+k , k∈ℤk . We want none of these values to belong to BsB_s. A shift that places some b∈Bsb∈ B_s on a boundary must satisfy as=b−kℓa_s=b-k for some integer k. Thus all such shifts belong to the countable set b−kℓ:b∈Bs,k∈ℤ.\b-k :b∈ B_s,\ k \. Since the interval [0,ℓ)[0, ) is uncountable, it contains a shift asa_s outside this set. For this choice, every boundary value as+kℓa_s+k lies outside BsB_s. By the definition of BsB_s, μx∈E:(vt,i(x))s=as+kℓ=0(t∈,i∈[n],k∈ℤ).μ\x∈ E:(v_t,i(x))_s=a_s+k \=0 (t ,\ i∈[n],\ k ). Choose each coordinate shift in this way and fix a=(a1,…,ady)a=(a_1,…,a_d_y) for the entire family. The union of all grid-boundary hyperplanes is ℋ:=⋃s=1dy⋃k∈ℤv∈ℝdy:vs=as+kℓ.H:= _s=1^d_y\ _k \v ^d_y:v_s=a_s+k \. By the countable union bound, the set of inputs for which any vt,i(x)v_t,i(x) lies on ℋH has measure zero. Step 3: Remove a thin neighborhood of the boundaries. For 0<τ<ℓ/20<τ< /2, remove the inputs Nτ:=x∈E:mint∈,i∈[n]dist(vt,i(x),ℋ)<τ.N_τ:= \x∈ E: _t ,\,i∈[n]dist(v_t,i(x),H)<τ \. The set ℋH is closed: for each coordinate s, the lattice as+ℓℤa_s+ is closed, and there are only finitely many coordinates. The displayed minimum is continuous, so NτN_τ is relatively open in E. As τ↓0τ 0, these sets decrease to the zero-measure event that some node state lies on ℋH. Since μ is finite, continuity of measure from above gives μ(Nτ)→0μ(N_τ)→ 0. Choose τ sufficiently small that μ(Nτ)<γμ(N_τ)<γ and, if the subsets CaC_a are prescribed, also μ(Nτ)<12minaμ(Ca).μ(N_τ)< 12 _aμ(C_a). Then E′:=E∖NτE :=E N_τ is compact, loses less than γ measure, and retains positive measure in every prescribed CaC_a. Node permutations reorder the node states without changing their internal coordinates. Since every node is tested against the same fixed grid, equivariance gives mint,idist(vt,i(πx),ℋ)=mint,idist(vt,π−1(i)(x),ℋ)=mint,jdist(vt,j(x),ℋ). _t,idist(v_t,i(π x),H)= _t,idist(v_t,π^-1(i)(x),H)= _t,jdist(v_t,j(x),H). Thus NτN_τ and E′=E∖NτE =E N_τ are permutation invariant, regardless of whether μ is invariant. Step 4: Obtain separated boxes and bound the center error. Every retained node state lies at distance at least τ from every grid boundary, and therefore belongs to one of the shrunken boxes Qz:=a+ℓz+[τ,ℓ−τ]dy.Q_z:=a+ z+[τ, -τ]^d_y. Keep only boxes containing a node state vt,i(x)v_t,i(x) for some x∈E′x∈ E , and enumerate them as Q1,…,QJQ_1,…,Q_J. Only finitely many are needed because all node-state ranges are bounded. These boxes are compact, and distinct boxes have distance at least 2τ>02τ>0. Their diameters satisfy diam(Qj)=dy(ℓ−2τ)<η.diam(Q_j)= d_y( -2τ)<η. Consequently, the n node states of each retained Vt(x)V_t(x) belong to different boxes. Finally, each node state is within dy(ℓ−2τ)/2 d_y( -2τ)/2 of its box center. Summing the squared errors over the n nodes gives ‖Vt(x)−b(Vt(x))‖≤ndy2(ℓ−2τ)<ω.\|V_t(x)-b(V_t(x))\|≤ nd_y2( -2τ)<ω. This proves all the claims. ∎ In the application below, E′E consists of retained graph–initialization pairs, and Vt(G,y0)=TGt(y0)V_t(G,y_0)=T_G^t(y_0) for finitely many times. The boxes therefore provide one fixed collection of address regions for all retained graphs, initializations, and these times. Subsection B.4 then constructs terminal dynamics that keep each node in its box and contract toward its center, preserving the addressing condition while ensuring convergence near the selected solution. B.4. Proof of the main result Given the preparations, we are now ready to prove Theorem 5. Proof of Theorem 5. Choose positive probability budgets satisfying δloc+δstop+δsep+δgrid<δ. _ loc+ _ stop+ _ sep+ _ grid<δ. These budgets allow us to exclude graph–initialization pairs with total probability less than δ, while preserving positive probability for each labeled solution branch. The remaining pairs form a compact family whose ideal trajectories all approach their selected solutions within a common finite number of steps. During these steps, different nodes have uniformly separated states, which can be placed in small, separated boxes from one fixed collection shared across all retained inputs. We then construct an update rule that follows TGT_G exactly until the state is sufficiently close to its selected solution. From that point onward, each node moves halfway toward the center of its current box at every step. This ensures convergence while keeping different nodes distinguishable. Switching close to the solution and using small boxes controls both the final solution error and the deviation from the ideal update. Lemma 5 then realizes this rule with a single fixed GNN block, whose repeated application generates all retained trajectories. Step 1: Localize the initializations and the ideal trajectories. For G∈⋆G _ , the proof of Theorem 4 identifies the G-section of the joint bad set with ΣG∘ _G , even after taking the joint closure. Its switching-set nullity and the first two probability assumptions therefore give ℙ(Ω)=1P( )=1. For each label a, define Ba:=(G,y0)∈Ω:P(G,y0)=fa(G).B_a:=\(G,y_0)∈ :P(G,y_0)=f_a(G)\. These sets are Borel, since P and the regular branches are continuous on their domains. For G∈⋆G _ , the G-section of BaB_a is Vfa(G)(G)V_f_a(G)(G), and the section is empty otherwise. Integration of the indicator of BaB_a against the conditional kernel is measurable in G, so disintegration and the positive-basin assumption give, with ℙGP_G denoting the graph marginal, ℙ(Ba)=∫⋆μG(Vfa(G)(G))dℙG(G)>0.P(B_a)= _G_ _G(V_f_a(G)(G))\,dP_G(G)>0. A tight finite Borel measure on a metric space is inner regular by compact sets: restrict first to a compact set of arbitrarily high probability and use Borel regularity there. Hence choose compact sets C0⊆ΩC_0 and Ca⊆BaC_a B_a such that ℙ(C0)>1−δloc,ℙ(Ca)>0(a=1,…,M).P(C_0)>1- _ loc, (C_a)>0 (a=1,…,M). Finite symmetrization gives the compact invariant event E0:=⋃π∈nπ(C0∪⋃a=1MCa)⊆Ω,ℙ(E0)>1−δloc.E^0:= _π∈ S_nπ (C_0∪ _a=1^MC_a ) , (E^0)>1- _ loc. The compact sets CaC_a are the branch subsets whose mass we preserve. Write s(G,y0):=P(G,y0)s(G,y_0):=P(G,y_0). The union of the full ideal segments, :=(G,(1−a)y0+as(G,y0)):(G,y0)∈E0,a∈[0,1],C:= \ (G,(1-a)y_0+a\,s(G,y_0) ):(G,y_0)∈ E^0,\ a∈[0,1] \, is a compact subset of Ω : it is the continuous image of E0×[0,1]E^0×[0,1], and Appendix A shows that each segment stays in its initial strict basin. Set b∗:=min(G,y)∈β(G,y)∈(0,1),R0:=max(G,y0)∈E0‖y0−s(G,y0)‖<∞.b_*:= _(G,y) β(G,y)∈(0,1), R_0:= _(G,y_0)∈ E^0\|y_0-s(G,y_0)\|<∞. Equation (17) then gives the uniform bound (24) ‖TGt(y0)−s(G,y0)‖≤(1−b∗)tR0((G,y0)∈E0,t≥0).\|T_G^t(y_0)-s(G,y_0)\|≤(1-b_*)^tR_0 ((G,y_0)∈ E^0,\ t≥ 0). Step 2: Choose an entrance threshold with a positive margin. We use distance to the solution set to decide when to switch from the ideal update TG(⋅)T_G(·) to the final convergent update. By excluding a small-probability set of graph–initialization pairs, we ensure that this distance stays uniformly away from the chosen threshold at every time up to a common finite horizon. The resulting gap separates the states where the two updates are used, allowing us to combine them into a continuous rule. For x=(G,y0)∈E0x=(G,y_0)∈ E^0, define yt(x):=TGt(y0),dt(x):=dist(yt(x),S(G)).y_t(x):=T_G^t(y_0), d_t(x):=dist(y_t(x),S(G)). Recall that E0E^0 is a compact subset of Ω , so its graph projection is a compact subset of ⋆G_ . Each branch faf_a is locally Lipschitz on ⋆G_ and hence continuous on this projection. Moreover, y0(x)y_0(x) is a coordinate projection, and the recurrence yt+1(x)=TG(yt(x))y_t+1(x)=T_G(y_t(x)) uses the jointly continuous map (G,y)↦TG(y)(G,y) T_G(y). Induction therefore shows that yt(x)y_t(x) is continuous in x=(G,y0)x=(G,y_0) for every finite t. Consequently, dt(x)=mina=1,…,M‖yt(x)−fa(G)‖d_t(x)= _a=1,…,M\|y_t(x)-f_a(G)\| is continuous as the minimum of finitely many continuous functions. The distance dtd_t is also permutation invariant. Indeed, equivariance of T gives yt(πx)=πyt(x)y_t(π x)=π y_t(x), while S(πG)=πS(G)S(π G)=π S(G). Since node permutations preserve Euclidean distances, dt(πx)=dist(πyt(x),πS(G))=dist(yt(x),S(G))=dt(x).d_t(π x)=dist(π y_t(x),π S(G))=dist(y_t(x),S(G))=d_t(x). Finally, Appendix A shows that the ideal trajectory remains in the strict basin of its selected solution s(x)s(x). Thus s(x)s(x) remains the unique nearest solution, and dt(x)=‖yt(x)−s(x)‖.d_t(x)=\|y_t(x)-s(x)\|. Choose the radius before choosing a time horizon. For each fixed t∈ℕ0t _0, there are at most countably many distance values b≥0b≥ 0 for which ℙx∈E0:dt(x)=b>0,P\x∈ E^0:d_t(x)=b\>0, by the same countability argument used in Step 2 of the proof of Lemma 7. Since there are only countably many times t, the collection of all such values across all finite times is still countable. The interval (0,minε,ρ/2)(0, \ ,ρ/2\) is uncountable, so we can choose a radius r in this interval outside that collection. This ensures ℙx∈E0:dt(x)=r=0for every t∈ℕ0.P\x∈ E^0:d_t(x)=r\=0 every t _0. Now choose H∈ℕ0H _0 with (1−b∗)HR0<r/4.(1-b_*)^HR_0<r/4. In particular, every trajectory from E0E^0 is within r/4r/4 of its selected solution by time H. For ζ>0ζ>0, the events Nstop(ζ):=x∈E0:|dt(x)−r|<ζ for some 0≤t≤HN_ stop(ζ):=\x∈ E^0:|d_t(x)-r|<ζ for some 0≤ t≤ H\ decrease to a null event as ζ↓0ζ 0. Choose 0<ζ<r/20<ζ<r/2 such that ℙ(Nstop(ζ))<minδstop,12minaℙ(Ca).P(N_ stop(ζ))< \ _ stop, 12 _aP(C_a) \. The event Estop:=E0∖Nstop(ζ)E_ stop:=E^0 N_ stop(ζ) is compact and permutation invariant, with ℙ(Estop)>1−δloc−δstop,ℙ(Ca∩Estop)>12ℙ(Ca)>0(a=1,…,M).P(E_ stop)>1- _ loc- _ stop, (C_a∩ E_ stop)> 12P(C_a)>0 (a=1,…,M). On EstopE_ stop, every ideal trajectory stays at least ζ away from the distance threshold r through time H: |dt(x)−r|≥ζ(x∈Estop, 0≤t≤H).|d_t(x)-r|≥ζ (x∈ E_ stop,\ 0≤ t≤ H). Thus, at each of these times, the distance to the solution set is either at most r−ζr-ζ or at least r+ζr+ζ. On EstopE_ stop, define the entrance time h(x):=min0≤t≤H:dt(x)≤r−ζ.h(x):= \0≤ t≤ H:d_t(x)≤ r-ζ\. It exists because dH<r/4<r−ζd_H<r/4<r-ζ. Every retained trajectory satisfies (25) dt(x)≥r+ζ(0≤t<h(x)),dh(x)(x)≤r−ζ.d_t(x)≥ r+ζ (0≤ t<h(x)), d_h(x)(x)≤ r-ζ. Thus it first visits the inner neighborhood of radius r−ζr-ζ after remaining outside the larger neighborhood of radius r+ζr+ζ. For each k∈0,…,Hk∈\0,…,H\, its entrance-time level set is Lk:=Estop∩dk≤r−ζ∩⋂t<kdt≥r+ζ.L_k:=E_ stop∩\d_k≤ r-ζ\∩ _t<k\d_t≥ r+ζ\. Thus LkL_k consists of graph–initialization pairs, (G,y0)(G,y_0), whose ideal trajectories first enter the (r−ζ)(r-ζ)-neighborhood of the solution set at time k. This is compact and invariant. When k=0k=0, there are no pre-entrance inequalities. The sets LkL_k form a finite partition of EstopE_ stop. The horizon H is a uniform upper bound on this retained family, whereas h(x)h(x) may be smaller and depends on the trajectory. Step 3: Obtain one collection of distinct node-address boxes. Apply Lemma 6 to EstopE_ stop, horizon H, loss δsep _ sep, and the subsets Ca∩EstopC_a∩ E_ stop. It gives a compact invariant Esep⊆EstopE_ sep E_ stop and η>0η>0 such that ℙ(Estop∖Esep)<δsep,ℙ(Ca∩Esep)>0,P(E_ stop E_ sep)< _ sep, (C_a∩ E_ sep)>0, and, uniformly over all retained inputs and times up to H, distinct nodes have states at least η apart: ∥[TGt(y0)]i−[TGt(y0)]j∥≥η((G,y0)∈Esep, 0≤t≤H,i≠j).\|[T_G^t(y_0)]_i-[T_G^t(y_0)]_j\|≥η ((G,y_0)∈ E_ sep,\ 0≤ t≤ H,\ i≠ j ). Then choose any 0<ω<ζ0<ω<ζ, apply Lemma 7 on EsepE_ sep with the restricted measure ℙP, tolerance ω, loss δgrid _ grid, subsets Ca∩EsepC_a∩ E_ sep, and state maps Vt(G,y0):=TGt(y0),t=0,…,H.V_t(G,y_0):=T_G^t(y_0), t=0,…,H. The maps are continuous and equivariant, and the row-separation hypothesis was just established. The lemma gives boxes Q1,…,QJQ_1,…,Q_J and a compact invariant final event (26) E:=Eε,ρ,δ⊆Esep,ℙ(E)>1−δ,ℙ(Ca∩E)>0(a=1,…,M).E:=E_ ,ρ,δ E_ sep, (E)>1-δ, (C_a∩ E)>0 (a=1,…,M). The probability bound follows by adding the four losses. For every x∈Ex∈ E and t∈0,…,Ht∈\0,…,H\, there exist pairwise distinct indices j1,…,jn∈[J]j_1,…,j_n∈[J] such that (27) [yt(x)]i∈Qji(i=1,…,n).[y_t(x)]_i∈ Q_j_i (i=1,…,n). Replacing each node state by the center of its box changes the full state by less than ω: ‖yt(x)−(bj1,…,bjn)‖<ω, \|y_t(x)-(b_j_1,…,b_j_n) \|<ω, where bjb_j is the center of QjQ_j. The indices jij_i may depend on (x,t)(x,t), while the boxes and their centers are fixed for the whole family. Using the final retained event E, define two families of graph–state pairs: Ktr K_ tr :=⋃k=0H⋃t=0k−1(G,yt(x)):x=(G,y0)∈E∩Lk, := _k=0^H\ _t=0^k-1\(G,y_t(x)):x=(G,y_0)∈ E∩ L_k\, Kent K_ ent :=⋃k=0H(G,yk(x)):x=(G,y0)∈E∩Lk. := _k=0^H\(G,y_k(x)):x=(G,y_0)∈ E∩ L_k\. Here KtrK_ tr collects the states of each retained trajectory before its first visit to the (r−ζ)(r-ζ)-neighborhood of S(G)S(G); these are the transient states. The family KentK_ ent collects the state at that first visit for each trajectory; these are the entrance states. By the definition of LkL_k, dist(y,S(G))≥r+ζon Ktr,dist(y,S(G))≤r−ζon Kent.dist(y,S(G))≥ r+ζ K_ tr, (y,S(G))≤ r-ζ K_ ent. Both families are compact, being finite unions of continuous images of the compact sets E∩LkE∩ L_k. They are permutation invariant because these sets are invariant and the maps yty_t are equivariant. In every state from either family, different nodes belong to different boxes from the same collection Q1,…,QJQ_1,…,Q_J—i.e., see (27). Hence, KtrK_ tr and KentK_ ent admit robust node-state symmetry breaking (Definition 6). Step 4: Construct the terminal dynamics and realize both regimes. Let Q:=⋃j=1JQjQ:= _j=1^JQ_j, and define b(v):=bj(v∈Qj)b(v):=b_j~(v∈ Q_j), where bjb_j is the center of QjQ_j. This map is continuous on Q because it is constant on each box and the boxes are positively separated. For a state y with every node state yi∈Qy_i∈ Q, define b(y):=(b(y1),…,b(yn)),Ψ(y):=b(y)+12(y−b(y)).b(y):=(b(y_1),…,b(y_n)), (y):=b(y)+ 12(y-b(y)). The update Ψ moves each node state halfway toward the center of its current box. Since the boxes are convex, each node remains in that same box, so b(Ψ(y))=b(y).b( (y))=b(y). Repeated updates therefore converge to the fixed center pattern b(y)b(y), halving each node’s distance to its box center at every step. For each entrance state, include the entire line segment joining it to its box-center pattern. These segments form the terminal family: (28) Kterm:=(G,b(y)+a(y−b(y))):(G,y)∈Kent,a∈[0,1].K_ term:= \ (G,b(y)+a(y-b(y)) ):(G,y)∈ K_ ent,\ a∈[0,1] \. Here a=1a=1 gives the entrance state y, while a=0a=0 gives the center pattern b(y)b(y) toward which the terminal updates converge. Fix an entrance state y. Its node states belong to distinct boxes: there are pairwise distinct indices j1,…,jnj_1,…,j_n such that yi∈Qjiy_i∈ Q_j_i. For any state z=b(y)+a(y−b(y))z=b(y)+a(y-b(y)) on the segment, zi=bji+a(yi−bji)∈Qji,a∈[0,1],z_i=b_j_i+a(y_i-b_j_i)∈ Q_j_i, a∈[0,1], because QjiQ_j_i is convex and contains both yiy_i and its center bjib_j_i. Thus each node remains in the box it occupied at entrance, and different nodes continue to occupy different boxes—i.e., KtermK_ term admits robust node-state symmetry breaking (Definition 6). The family KtermK_ term is compact as a continuous image of Kent×[0,1]K_ ent×[0,1], and it is permutation invariant. Moreover, Ψ(z)=b(y)+a2(y−b(y))∈Kterm(G), (z)=b(y)+ a2(y-b(y))∈ K_ term(G), where Kterm(G):=w:(G,w)∈KtermK_ term(G):=\w:(G,w)∈ K_ term\. Thus every terminal update stays in this family, which contains all terminal iterates and their limits. For a terminal state z=b(y)+a(y−b(y))z=b(y)+a(y-b(y)) arising from an entrance state y=yh(x)(x)y=y_h(x)(x), let s=s(x)s=s(x) be that initialization’s selected solution. The grid bound and (25) imply (29) ‖z−y‖<ω,‖z−b(y)‖<ω,‖z−s‖≤r−ζ+ω<r.\|z-y\|<ω, \|z-b(y)\|<ω, \|z-s\|≤ r-ζ+ω<r. These estimates hold also at both endpoints of the segment. In particular, the solution-set distance on KtermK_ term is below r, whereas it is at least r+ζr+ζ on KtrK_ tr. The two compact families, KtermK_ term and KtrK_ tr, are disjoint and, when both are nonempty, positively separated. Set K:=Ktr∪KtermK:=K_ tr∪ K_ term and prescribe (30) F∗(G,y):=TG(y),(G,y)∈Ktr,Ψ(y),(G,y)∈Kterm.F^*(G,y):= casesT_G(y),&(G,y)∈ K_ tr,\\[2.84526pt] (y),&(G,y)∈ K_ term. cases The two pieces are continuous and equivariant, and their domains are disjoint compact subsets of K. Thus F∗F^* is a well-defined continuous equivariant map on K. Since KtrK_ tr and KtermK_ term admit robust node-state symmetry breaking, K admits robust state-based symmetry breaking, with address regions j=QjA_j=Q_j. Apply Lemma 5 with u=0u=0, v=1v=1, and F1=F∗F_1=F^*. It gives one MP-GNN satisfying (31) GNN(G,y)=F∗(G,y)((G,y)∈K).GNN(G,y)=F^*(G,y) ((G,y)∈ K). The block is defined on the entire input domain by the extension steps in Lemma 5. Its proof extends the address encoder to all of ℝdyR^d_y and uses the Tietze extension theorem to extend the decoder from the compact reconstruction image to its entire Euclidean ambient space. The resulting GNN is therefore defined for every (G,y)∈×n(G,y) ×Y_n, while its agreement with F⋆F is guaranteed on K. In particular, it already realizes Ψ on KtermK_ term, so no separate extension of Ψ is needed. Both regimes are implemented by the same fixed block using only the current input (G,y)(G,y), without an additional time input or stopping flag. Step 5: Verify the full trajectory and branch capture. Fix x=(G,y0)∈Ex=(G,y_0)∈ E and write h=h(x)h=h(x). For t<ht<h, each ideal update from the transient family leads either to its next transient state or to the entrance state. From the entrance state onward, Ψ stays in KtermK_ term. Thus exact realization and induction give (32) y^t=TGt(y0)(0≤t≤h),y^h+k=b(y^h)+2−k(y^h−b(y^h))(k≥0). y_t=T_G^t(y_0) (0≤ t≤ h), y_h+k=b( y_h)+2^-k ( y_h-b( y_h) ) (k≥ 0). This includes the case h=0h=0. The trajectory converges to y^∞=b(y^h) y_∞=b( y_h), and (29) gives ‖y^∞−s(G,y0)‖≤r−ζ+ω<r<ε.\| y_∞-s(G,y_0)\|≤ r-ζ+ω<r< . The learned update equals TGT_G on the transient. On the terminal family, b(z)=b(y)b(z)=b(y) for its generating entrance state, so ‖Ψ(z)−z‖=12‖z−b(z)‖<ω/2.\| (z)-z\|= 12\|z-b(z)\|<ω/2. Using (18) and (29), uniformly on KtermK_ term, ‖GNN(G,z)−TG(z)‖ \|GNN(G,z)-T_G(z)\| ≤‖Ψ(z)−z‖+‖z−TG(z)‖ ≤\| (z)-z\|+\|z-T_G(z)\| ≤12‖z−b(z)‖+‖z−s‖ ≤ 12\|z-b(z)\|+\|z-s\| ≤ω/2+r−ζ+ω<r+ω/2<2r<ρ. ≤ω/2+r-ζ+ω<r+ω/2<2r<ρ. The right-hand bounds are constants independent of the trajectory and time, proving the supremum bound in (20). The preceding argument establishes convergence and the uniform operator-error bound for every (G,y0)∈E(G,y_0)∈ E. In particular, every retained trajectory converges to a limit within ε of its selected solution, and hence of S(G)S(G). Together with (26), this gives ℙ(y^t→y^∞,dist(y^∞,S(G))<ε)≥ℙ(E)≥1−δ.P\! ( y_t→ y_∞,\ dist( y_∞,S(G))< ) (E)≥ 1-δ. To verify branch capture, fix a label a. By construction, Ca⊆BaC_a B_a, so every (G,y0)∈Ca∩E(G,y_0)∈ C_a∩ E satisfies y0∈Vfa(G)(G),s(G,y0)=fa(G).y_0∈ V_f_a(G)(G), s(G,y_0)=f_a(G). The limit estimate therefore gives ‖y^∞−fa(G)‖=‖y^∞−s(G,y0)‖<ε((G,y0)∈Ca∩E).\| y_∞-f_a(G)\|=\| y_∞-s(G,y_0)\|< ((G,y_0)∈ C_a∩ E). Since each restriction preserved positive probability in CaC_a, ℙ(Y0∈Vfa(G)(G),y^t→y^∞,∥y^∞−fa(G)∥<ε)≥ℙ(Ca∩E)>0.P\! (Y_0∈ V_f_a(G)(G),\ y_t→ y_∞,\ \| y_∞-f_a(G)\|< ) (C_a∩ E)>0. Thus every labeled branch has positive capture probability under the joint law of graphs and initializations. Finally, the graph projection Kε,ρ,δ:=projEK_ ,ρ,δ:=proj_GE is compact because it is the continuous image of the compact set E. Lemma 5 gives a single GNN that agrees with F⋆F on all of K. For every retained pair (G,y0)∈E(G,y_0)∈ E, the entire learned trajectory stays in K. Thus the same GNN, with all component maps fixed, generates every retained trajectory by repeated application. Moreover, the shared nodewise maps and sum aggregations commute with node relabeling, giving the permutation equivariance: GNN(πG,πy)=πGNN(G,y)GNN(π G,π y)= (G,y). This concludes the proof. ∎ Appendix C An obstruction to almost-sure validity for contractive equivariant dynamics We use the graph space ℰnE_n, output space n=ℝn×qY_n=R^n× q, Euclidean norms, and permutation actions of Appendix A. In particular, πy=Pπyπ y=P_πy, and the graph domain ⊆ℰnG _n is closed under node relabeling. The Euclidean norm on a matrix is its Frobenius norm. All continuity and neighborhood statements on G refer to its relative topology. For a graph G∈G , define its automorphism group and the corresponding fixed subspace by ΓG:=π∈n:πG=G,Fix(ΓG):=y∈n:πy=y for every π∈ΓG. _G:=\π∈ S_n:π G=G\, ( _G):=\y _n:π y=y for every π∈ _G\. The latter is a closed linear subspace, since Fix(ΓG)=⋂π∈ΓGker(y↦Pπy−y).Fix( _G)= _π∈ _G (y P_πy-y ). Let ν be a Borel probability measure on G. We write G0∈suppνG_0 ν when every relatively open neighborhood of G0G_0 has positive ν-measure. Thus a full-support graph law satisfies this condition at every G0∈G_0 . C.1. The operator-level obstruction Theorem 6 (Positive-probability failure of continuous contractive equivariant dynamics). Let S(G)=f1(G),…,fM(G)S(G)=\f_1(G),…,f_M(G)\ be a permutation-equivariant finite-branch solution map in the sense of Definition 1. For this theorem, assume additionally that every branch fa:→nf_a:G _n is continuous. Suppose there is a graph G0∈suppνG_0 ν such that (33) S(G0)∩Fix(ΓG0)=∅.S(G_0) ( _G_0)= . Let F:×n→n,(G,y)↦FG(y),F:G×Y_n _n, (G,y) F_G(y), be jointly continuous and permutation equivariant. Assume that, for some κ∈[0,1)κ∈[0,1), (34) ‖FG(y)−FG(z)‖≤κ‖y−z‖,FπG(πy)=πFG(y),\|F_G(y)-F_G(z)\|≤κ\|y-z\|, F_π G(π y)=π F_G(y), for all G∈G , y,z∈ny,z _n, and node permutations π. Then the following statements hold. (1) A unique, initialization-independent limit. For every G∈G , there is a unique fixed point pF(G)∈np_F(G) _n. Every trajectory yt+1=FG(yt)y_t+1=F_G(y_t) satisfies (35) ‖yt−pF(G)‖≤κt‖y0−pF(G)‖,t≥1.\|y_t-p_F(G)\|≤κ^t\|y_0-p_F(G)\|, t≥ 1. The map pF:→np_F:G _n is continuous and permutation equivariant. In particular, pF(G)∈Fix(ΓG)(G∈).p_F(G) ( _G) (G ). (2) Failure on a neighborhood of the obstruction graph. The number (36) γ0:=min1≤a≤Mdist(fa(G0),Fix(ΓG0)) _0:= _1≤ a≤ Mdist (f_a(G_0),Fix( _G_0) ) is strictly positive and does not depend on F. For every ε∈(0,γ0) ∈(0, _0), there is a relatively open neighborhood VF,εV_F, of G0G_0 such that (37) dist(pF(G),S(G))>ε(G∈VF,ε).dist (p_F(G),S(G) )> (G∈ V_F, ). (3) Positive-probability failure under random initialization. Let (G,Y0)(G,Y_0) have any Borel probability law ℙP on ×nG×Y_n with graph marginal ν, and define Yt+1=FG(Yt),t≥0.Y_t+1=F_G(Y_t), t≥ 0. For every ε∈(0,γ0) ∈(0, _0), set δF,ε:=ν(VF,ε)>0 _F, :=ν(V_F, )>0. Then YtY_t converges for every (G,Y0)(G,Y_0), and its limit Y∞=pF(G)Y_∞=p_F(G) satisfies (38) ℙ(dist(Y∞,S(G))>ε)≥δF,ε>0.P\! (dist(Y_∞,S(G))> )≥ _F, >0. Consequently, ℙ(Y∞∈S(G))<1.P (Y_∞∈ S(G) )<1. No independence, density, or full-support assumption is needed for the initialization law. Proof. We verify the three conclusions. Step 1: Existence, uniqueness, and convergence for each graph. Fix G∈G and y0∈ny_0 _n. Repeated use of (34) gives ‖yt+1−yt‖≤κt‖y1−y0‖.\|y_t+1-y_t\|≤κ^t\|y_1-y_0\|. For integers r>s≥0r>s≥ 0, it follows that ‖yr−ys‖≤∑t=sr−1κt‖y1−y0‖.\|y_r-y_s\|≤ _t=s^r-1κ^t\|y_1-y_0\|. The right-hand side tends to zero as s→∞s→∞, uniformly in r>sr>s, because the geometric series is summable. Hence (yt)(y_t) is Cauchy. The space nY_n is complete, so yty_t has a limit p. Continuity of FGF_G yields FG(p)=limt→∞FG(yt)=limt→∞yt+1=p.F_G(p)= _t→∞F_G(y_t)= _t→∞y_t+1=p. If p and p~ p are fixed points, then ‖p−p~‖≤κ‖p−p~‖,\|p- p\|≤κ\|p- p\|, which forces p=p~p= p because κ<1κ<1. Write this unique fixed point as pF(G)p_F(G). Its uniqueness makes the limit independent of y0y_0, and iteration of ‖yt+1−pF(G)‖≤κ‖yt−pF(G)‖\|y_t+1-p_F(G)\|≤κ\|y_t-p_F(G)\| proves (35). If κ=0κ=0, the trajectory reaches the fixed point after one step. Step 2: Continuity of the fixed-point map. For any G,H∈G,H , the fixed-point equations and (34) imply ‖pF(G)−pF(H)‖ \|p_F(G)-p_F(H)\| =‖FG(pF(G))−FH(pF(H))‖ =\|F_G(p_F(G))-F_H(p_F(H))\| ≤κ‖pF(G)−pF(H)‖+‖FG(pF(H))−FH(pF(H))‖. ≤κ\|p_F(G)-p_F(H)\|+\|F_G(p_F(H))-F_H(p_F(H))\|. Therefore (39) ‖pF(G)−pF(H)‖≤‖FG(pF(H))−FH(pF(H))‖1−κ.\|p_F(G)-p_F(H)\|≤ \|F_G(p_F(H))-F_H(p_F(H))\|1-κ. Fixing H and letting G→HG→ H, joint continuity of F makes the numerator tend to zero. This proves continuity of pFp_F on G. Step 3: Equivariance and preservation of graph automorphisms. By definition, pF(πG)p_F(π G) is a fixed point of πGπ G. For any node permutation π, FπG(πpF(G))=πFG(pF(G))=πpF(G)F_π G(π p_F(G))=π F_G(p_F(G))=π p_F(G) which implies πpF(G)π p_F(G) is also a fixed point of πGπ G. By uniqueness of the fixed point at πGπ G, pF(πG)=πpF(G).p_F(π G)=π p_F(G). If π∈ΓGπ∈ _G, then πG=Gπ G=G, and consequently πpF(G)=pF(G)π p_F(G)=p_F(G). Thus pF(G)∈Fix(ΓG)p_F(G) ( _G). Step 4: A positive error persists on a graph neighborhood. The fixed subspace Fix(ΓG0)Fix( _G_0) is closed. By (33), each of the finitely many points fa(G0)f_a(G_0) lies outside this subspace. Its distance to the subspace is therefore strictly positive, and the finite minimum in (36) satisfies γ0>0 _0>0. Step 3 gives ‖pF(G0)−fa(G0)‖≥γ0(1≤a≤M).\|p_F(G_0)-f_a(G_0)\|≥ _0 (1≤ a≤ M). Fix ε∈(0,γ0) ∈(0, _0) and put η:=(γ0−ε)/2>0η:=( _0- )/2>0. By continuity of pFp_F and of the finitely many branches, there is a relatively open neighborhood VF,εV_F, of G0G_0 such that ∥pF(G)−pF(G0)∥<η,max1≤a≤M∥fa(G)−fa(G0)∥<η(G∈VF,ε).\|p_F(G)-p_F(G_0)\|<η, _1≤ a≤ M\|f_a(G)-f_a(G_0)\|<η (G∈ V_F, ). For every such G and every a, the reverse triangle inequality gives ‖pF(G)−fa(G)‖ \|p_F(G)-f_a(G)\| ≥‖pF(G0)−fa(G0)‖−‖pF(G)−pF(G0)‖−‖fa(G)−fa(G0)‖ ≥\|p_F(G_0)-f_a(G_0)\|-\|p_F(G)-p_F(G_0)\|-\|f_a(G)-f_a(G_0)\| >γ0−2η=ε. > _0-2η= . Taking the minimum over the finite branch set proves (37). Step 5: The probability statement. Because G0∈suppνG_0 ν, δF,ε=ν(VF,ε)>0 _F, =ν(V_F, )>0. The function G⟼dist(pF(G),S(G))=min1≤a≤M‖pF(G)−fa(G)‖G (p_F(G),S(G))= _1≤ a≤ M\|p_F(G)-f_a(G)\| is continuous, so the events in the statement are measurable. By Step 1, Y∞=pF(G)Y_∞=p_F(G) for every initial state. Consequently, G∈VF,ε⊆dist(Y∞,S(G))>ε.\G∈ V_F, \ \dist(Y_∞,S(G))> \. Taking ℙP-probabilities and using its graph marginal ν proves (38). The last assertion follows because the event on the right contains no valid limit. ∎ Remark (Scope of the assumptions). The graph domain need not be compact or connected, and the graph law need not be permutation invariant. Full support of ν is sufficient but not necessary: the proof only uses G0∈suppνG_0 ν. The neighborhood argument only needs continuity of the branches at G0G_0; continuity on G is imposed above to ensure that all global error events are measurable without additional assumptions. In particular, the branch assumption holds on a domain where the branches are locally Lipschitz, as in the regular stable setting of Appendix A. The failure is uniform over initial states on VF,εV_F, only at the level of their limits; no common finite convergence time over the unbounded space nY_n is asserted. C.2. Application to independently weighted triangles Fix 0<w−<w+<∞0<w_-<w_+<∞ and let W=[w−,w+]3W=[w_-,w_+]^3. For w=(w12,w13,w23)∈Ww=(w_12,w_13,w_23)∈ W, define A(w)=(0w12w13w120w23w13w230),G(w)=(A(w),0).A(w)= pmatrix0&w_12&w_13\\ w_12&0&w_23\\ w_13&w_23&0 pmatrix, G(w)=(A(w),0). Here n=3n=3, de=dv=1d_e=d_v=1, and X=0∈ℝ3×1X=0 ^3× 1. Set △:=G(w):w∈W,LG:=diag(A)−A,HG:=I3+LG.G_ :=\G(w):w∈ W\, L_G:=diag(A1)-A, H_G:=I_3+L_G. We take q=2q=2 and write y=[ub]∈3=ℝ3×2y=[u\ b] _3=R^3× 2, where u,b∈ℝ3u,b ^3 are its two columns. The single-source diffusion solution set is (40) S(G)=f1(G),f2(G),f3(G),fk(G):=[HG−1ekek].S(G)=\f_1(G),f_2(G),f_3(G)\, f_k(G):=[H_G^-1e_k\ \ e_k]. Thus a valid output satisfies HGu=bH_Gu=b, with b a one-hot source indicator. Corollary 1 (Impossibility under independent edge weights). Let the three edge weights be independent and uniform on [w−,w+][w_-,w_+], and let ν be the induced law on △G_ . Suppose F:△×3→3F:G_ ×Y_3 _3 satisfies the continuity, equivariance, and global contraction assumptions of Theorem 6. For every 0<ε<23,0< < 23, there is a number δF,ε>0 _F, >0 such that, under any joint law of (G,Y0)(G,Y_0) with graph marginal ν, the recurrence converges and ℙ(dist(Y∞,S(G))>ε)≥δF,ε.P\! (dist(Y_∞,S(G))> )≥ _F, . In particular, no such operator converges to a valid solution almost surely over independently weighted graphs. The same conclusion holds for any full-support Borel law on △G_ . Proof. We verify the solution-map, support, and symmetry hypotheses. Step 1: The graph family and solution branches. Relabeling the nodes permutes the three edge coordinates, so △G_ is permutation invariant. For every v∈ℝ3v ^3, v⊤HGv=‖v‖22+∑1≤i<j≤3wij(vi−vj)2≥‖v‖22.v H_Gv=\|v\|_2^2+ _1≤ i<j≤ 3w_ij(v_i-v_j)^2≥\|v\|_2^2. Hence HGH_G is invertible and ‖HG−1‖op≤1\|H_G^-1\|_op≤ 1. The map G↦HGG H_G is linear in the edge matrix up to the constant identity. For any G,G~∈△G, G _ , the resolvent identity gives HG−1−HG~−1=HG−1(HG~−HG)HG~−1.H_G^-1-H_ G^-1=H_G^-1(H_ G-H_G)H_ G^-1. Consequently, ‖HG−1−HG~−1‖op≤‖HG~−HG‖op.\|H_G^-1-H_ G^-1\|_op≤\|H_ G-H_G\|_op. Since linear maps between finite-dimensional Euclidean spaces are globally Lipschitz, every branch in (40) is globally Lipschitz. Furthermore, HπG=PπHGPπ⊤,HπG−1=PπHG−1Pπ⊤.H_π G=P_πH_GP_π , H_π G^-1=P_πH_G^-1P_π . Using Pπek=eπ(k)P_πe_k=e_π(k) yields fπ(k)(πG)=[PπHG−1ekPπek]=πfk(G).f_π(k)(π G)=[P_πH_G^-1e_k\ \ P_πe_k]=π f_k(G). Thus Definition 1 holds with τπ(k)=π(k) _π(k)=π(k). The branches are distinct for every graph because their second columns are distinct; indeed, ‖fk(G)−fℓ(G)‖≥‖ek−eℓ‖2=2(k≠ℓ).\|f_k(G)-f_ (G)\|≥\|e_k-e_ \|_2= 2 (k≠ ). Their collapse partition consists of three singleton sets throughout the domain. Hence D=U=∅D=U= and (△)⋆=△(G_ )_ =G_ in the notation of Definition 2. Step 2: The graph law has full support. The map w↦G(w)w G(w) is a homeomorphism from W onto △G_ . In fact, ‖G(w)−G(w~)‖=2‖w−w~‖2.\|G(w)-G( w)\|= 2\,\|w- w\|_2. Every nonempty relatively open subset of the cube W has positive three-dimensional Lebesgue measure: it contains the intersection of W with a sufficiently small Euclidean ball, and this intersection has positive volume. The uniform product law has constant positive density on W. Its pushforward ν therefore has full support on △G_ . Step 3: A symmetric graph provides a fixed positive margin. Choose any c∈(w−,w+)c∈(w_-,w_+) and set G0=G(c,c,c)G_0=G(c,c,c). Every node permutation fixes G0G_0, so ΓG0=3 _G_0= S_3 and Fix(ΓG0)=[αβ]:α,β∈ℝ.Fix( _G_0)=\[ 1\ \ 1]:α,β \. For any k∈1,2,3k∈\1,2,3\ and β∈ℝβ , ‖ek−β‖22=(1−β)2+2β2=3(β−13)2+23≥23.\|e_k- 1\|_2^2=(1-β)^2+2β^2=3 (β- 13 )^2+ 23≥ 23. It follows that S(G0)∩Fix(ΓG0)=∅,γ0≥23.S(G_0) ( _G_0)= , _0≥ 23. Step 2 gives G0∈suppνG_0 ν. Applying Theorem 6 proves the claimed probability bound for every ε∈(0,2/3) ∈(0, 2/3) and for any initialization law. The proof uses only full support of the graph law, so the last assertion follows as well. ∎ C.3. Consequence for recurrent message-passing GNNs Corollary 2 (Contractive recurrent MP-GNNs). Consider a deterministic, autonomous, weight-tied MP-GNN block of the form (15), with continuous component maps and only yty_t carried between recurrent calls. Suppose its complete update satisfies ‖GNN(G,y)−GNN(G,z)‖≤κ‖y−z‖(G∈,y,z∈n)\|GNN(G,y)-GNN(G,z)\|≤κ\|y-z\| (G ,\ y,z _n) for some κ<1κ<1. Then the conclusions of Theorem 6 hold whenever its solution-map and graph-law hypotheses hold. In particular, Corollary 1 applies on △G_ , for every finite number of hidden layers and every choice of finite hidden widths. Proof. Set FG(y):=GNN(G,y)F_G(y):=GNN(G,y). The encoder is continuous in (G,y)(G,y). Inductively, if all hi(ℓ)h_i^( ) are continuous, then each message is continuous by continuity of MℓM_ , its finite sum is continuous, and hi(ℓ+1)h_i^( +1) is continuous by continuity of UℓU_ . The final finite sum and continuous output map preserve continuity. Thus F is jointly continuous. Under the simultaneous relabeling (G,y)↦(πG,πy)(G,y) (π G,π y), the encoder satisfies hπ(i)(0)(πG,πy)=hi(0)(G,y)h_π(i)^(0)(π G,π y)=h_i^(0)(G,y). If this identity holds at layer ℓ , reindexing the sender sum by j↦π(j)j π(j) gives the same identity for messages and then for layer ℓ+1 +1. The final global sum is invariant, and the node outputs are permuted. Hence FπG(πy)=πFG(y)F_π G(π y)=π F_G(y), as also stated in (16). The assumed bound is exactly the remaining contraction hypothesis. The two cited results now apply. ∎ Relation to the multi-attractor results. For the triangle family, the preceding proof established (△)⋆=△(G_ )_ =G_ . Consequently, Theorem 4 applies. If the initialization law has a density and full support on 3Y_3, every strict Voronoi basin has positive probability. Indeed, for each solution s, a sufficiently small open ball centered at s lies in Vs(G)V_s(G) because the other solutions are finitely many and distinct. The ideal operator therefore has almost-sure validity and positive capture probability for every solution at every fixed graph in △G_ . By contrast, every contractive operator above has one limit pF(G)p_F(G) for each fixed G. Its limiting law under any initialization distribution is a point mass. It cannot reach two distinct solutions with positive probability at the same graph. Moreover, Corollary 1 rules out exact almost-sure validity under the independent-weight graph law, even if coverage of multiple solutions is not required. Appendix D Single-source diffusion: protocol and additional results Datasets. Each instance is a complete undirected graph on three nodes. Its fixed input consists only of the three positive edge weights; static node features are identically zero and omitted from the implemented encoder. Draw w12,w23,w31w_12,w_23,w_31 independently from Unif[0.5,1.5]Unif[0.5,1.5]. Randomly permute node labels consistently with the adjacency matrix. No node indices, drawing coordinates, prescribed source, or exact diffusion fields are supplied as input features. We generate 500/100/100 train/validation/test instances. Model architecture. Across the three models, we use the same architecture: hi(0) h_i^(0) =enc(yi), =enc(y_i), mi(ℓ) m_i^( ) =∑j≠iMℓ(hi(ℓ),hj(ℓ),wij), = _j≠ iM_ (h_i^( ),h_j^( ),w_ij ), hi(ℓ+1) h_i^( +1) =Uℓ(hi(ℓ),mi(ℓ)),ℓ=0,1, =U_ (h_i^( ),m_i^( ) ), =0,1, (41) [Nθ(G,y)]i [N_θ(G,y)]_i =O(hi(2)),GNNθ(G,y)=P(ay+Nθ(G,y)). =O(h_i^(2)), _θ(G,y)=P(ay+N_θ(G,y)). Each component is a two-affine-layer MLP with a width-32 hidden layer and a SiLU activation between the affine layers. Messages use sum aggregation over the two neighbors. Layers within a block have separate parameters, but the entire block is shared across recurrent time t. The scalar skip a is also shared over nodes, coordinates, and time. The output projection is P([ub])=[uclip[0,1](b)]P([u\ b])=[u\ clip_[0,1](b)]; it does not impose the mass constraint or choose a source. Both columns of y0y_0 are i.i.d. (0,1)N(0,1) samples, independently across nodes and trajectories. Training. For supervised training, sample kGk_G uniformly from 1,2,3\1,2,3\ once per training graph, independently of the neural initialization, and retain sG=[HG−1ekGekG]s_G=[H_G^-1e_k_G\ \ e_k_G] for every subsequent visit. The target is not resampled between updates or epochs. We train the model by minimizing the loss function ‖yT−sG‖F2/6\|y_T-s_G\|_F^2/6 with T=10T=10. For energy based training, the training loss is EG(uT,bT)E_G(u_T,b_T), where energy function EGE_G is defined in (3) and T=10T=10. Our multi-attractor model is directly trained with this energy function without constraints on its weights, while contractive models are trained over the same energy function but applying the Lipschitz control technique described below. During training, a batch contains 32 different training graphs and eight starts per graph (256 trajectories); the final batch has 20 graphs (160 trajectories). With 500 graphs, there are 16 optimizer updates per epoch. We train 1,200 epochs (19,200 updates) with Adam, its default momentum parameters (0.9,0.999)(0.9,0.999), no weight decay, and no gradient clipping. The initial learning rate is 10−310^-3, multiplied by 0.20.2 after epoch 600 and by 0.050.05 relative to the initial rate after epoch 960. Lipschitz control. We enforce global contraction with respect to the recurrent state in the infinity norm. For a node-state array h, this means ‖h‖∞=maxi‖hi‖∞\|h\|_∞= _i\|h_i\|_∞. For a weight matrix, we use the corresponding induced norm ‖W‖∞=max∑sr|Wrs|\|W\|_∞= _r _s|W_rs|. Let σS(r)=r/(1+e−r) _ S(r)=r/(1+e^-r) denote SiLU. We use the fixed, conservative Lipschitz bound cS:=1+1/ec_ S:=1+1/e, since |σS′(r)|≤1+|r|e−|r|≤cS| _ S (r)|≤ 1+|r|e^-|r|≤ c_ S where r∈ℝr . Consequently, a two-affine-layer MLP g(z)=W2σS(W1z+β1)+β2g(z)=W_2 _ S(W_1z+ _1)+ _2 satisfies Lip∞(g)≤cS‖W2‖∞‖W1‖∞.Lip_∞(g)≤ c_ S\|W_2\|_∞\|W_1\|_∞. For a message MLP Mℓ(hi,hj,wij)M_ (h_i,h_j,w_ij), the graph is held fixed when measuring state sensitivity. Thus its first-layer norm is taken only over the combined columns multiplying hih_i and hjh_j; the edge-weight columns and all biases are excluded from this norm bound. In particular, LMℓL_M_ bounds the joint Lipschitz constant in [hi;hj][h_i;h_j], rather than a separate constant for each argument. Concatenation uses a maximum in the infinity norm: ‖[hi;hj]−[h~i;h~j]‖∞=max‖hi−h~i‖∞,‖hj−h~j‖∞.\|[h_i;h_j]-[ h_i; h_j]\|_∞= \\|h_i- h_i\|_∞,\|h_j- h_j\|_∞\. Therefore, summing the two incoming messages gives a bound 2LMℓ2L_M_ , with no additional factor for the two message arguments. The update MLP receives [hi;mi][h_i;m_i], giving the layer bound LUℓmax1,2LMℓL_U_ \1,2L_M_ \. Writing LencL_ enc and LOL_O for the encoder and output MLP bounds, the complete residual block satisfies Lipy(Nθ)≤LOLenc∏ℓ=01LUℓmax1,2LMℓ.Lip_y(N_θ)≤ L_OL_ enc _ =0^1L_U_ \1,2L_M_ \. We impose the following constraints once after parameter initialization and immediately after every Adam update, before the next forward pass. Set ρ=0.9ρ=0.9, c=cS−1/2c=c_ S^-1/2, and η=1−10−5η=1-10^-5. First clip a to [−ρ+0.05,ρ−0.05][-ρ+0.05,ρ-0.05], then use this clipped value to set the output-layer budget. The absolute row-sum caps and resulting MLP bounds are MLPFirst-layer capSecond-layer capLipschitz boundencηcηcLenc≤η2Mℓηc/2ηcLMℓ≤η2/2UℓηcηcLUℓ≤η2Oηcη(ρ−|a|)cLO≤η2(ρ−|a|) array[]c|c|c|cMLP&First-layer cap&Second-layer cap&Lipschitz bound\\ &η c&η c&L_ enc≤η^2\\ M_ &η c/2&η c&L_M_ ≤η^2/2\\ U_ &η c&η c&L_U_ ≤η^2\\ O&η c&η(ρ-|a|)c&L_O≤η^2(ρ-|a|) array For MℓM_ , the first-layer cap applies to the sum of absolute entries across both state-input blocks together, not to each block separately. All other caps apply to the full weight-matrix rows. For a constrained row or row segment w with cap B>0B>0, we rescale w←wmax1,‖w‖1/B.w← w \1,\|w\|_1/B\. Thus rows already satisfying the cap are unchanged; oversized rows are scaled radially to the boundary. This projection is an optimizer-side operation, performed without differentiating through it. Since 2LMℓ≤η2<12L_M_ ≤η^2<1, the two message-passing layers yield Lipy(Nθ)≤η8(ρ−|a|)Lip_y(N_θ)≤η^8(ρ-|a|). The coordinatewise clipping map P is nonexpansive in the infinity norm, so the full update GNNθ(G,y)=P(ay+Nθ(G,y))GNN_θ(G,y)=P(ay+N_θ(G,y)) satisfies Lipy(GNNθ)≤|a|+η8(ρ−|a|)≤ρ.Lip_y(GNN_θ)≤|a|+η^8(ρ-|a|)≤ρ. The constraints therefore certify contraction throughout training, rather than only after training or at sampled states. The multi-attractor and supervised models have no such projection. Single-target supervision itself does not guarantee a unique equilibrium. Additional visualizations. Figure 3 in the main text shows the evolution of the source coordinates bt∈ℝ3b_t ^3. Since the full recurrent state is yt=[utbt]y_t=[u_t\ b_t], Figure 4 complements this visualization by showing the diffusion coordinates ut∈ℝ3u_t ^3. The same pattern emerges: the multi-attractor trajectories form three distinct clusters, whereas the supervised and contractive controls each concentrate near a single point that remains visibly separated from all three exact diffusion solutions. Figure 4. u-counterpart of Figure 3: the same test graph, seed, 128 trajectories, methods, and time steps. Blue points show [(ut)1,(ut)2,(ut)3][(u_t)_1,(u_t)_2,(u_t)_3]; magenta circles mark the three true solutions. Appendix E Ising Experiments: Details and Additional Results This section provides details and additional results regarding the Ising model experiments. Datasets. We generate lattice graphs by networkx.triangular_lattice_graph with open boundaries, following [30]. Each undirected edge receives an independent uniform coupling from −0.5,−1,−1.5\-0.5,-1,-1.5\; both directed orientations are stored with equal weight. We generate graphs with controlled sizes to assess scalability. We consider four levels of sizes (see Table 4). Every size has 2,500 train-source and 500 test graphs; the final 250 train-source graphs are held out for validation. For each graph used in supervised training, Gurobi [25] provides one optimal spin configuration with a reported MIP gap of zero. Enumerating the full solution set is computationally infeasible. Labeling extra-large graphs can take more than 5,000 seconds per instance, so we omit supervised methods at that scale. Certified test references are still used to evaluate the extra-large graphs. Table 4. Graph sizes. Level Nodes Mean nodes Edges Mean edges Compact 6–55 26.972 9–135 61.068 Medium 136–253 189.174 360–693 510.918 Large 561–990 770.580 1584–2838 2195.904 Extra-large 2016–3916 2965.198 5859–11484 8667.174 Model Architecture. For the recurrent dynamics yt+1=GNNθ(G,yt)y_t+1=GNN_θ(G,y_t), we use a graph convolutional backbone augmented with physical node features. At each step, we compute for each node i (here N(i)N(i) represents the neighbors of i): (42) fi f_i =∑j∈N(i)Jijyj,ai=∑j∈N(i)|Jij|,qi=fi/max(ai,10−8), = _j∈ N(i)J_ijy_j, a_i= _j∈ N(i)|J_ij|, q_i=f_i/ (a_i,10^-8), (43) ci(β) c_i^(β) =tanh[∑j∈N(i)atanh(clip(tanh(βJij)yj,−1+10−6,1−10−6))max(ai,10−8)], = \! [ _j∈ N(i)atanh(clip( (β J_ij)y_j,-1+10^-6,1-10^-6)) (a_i,10^-8) ], with β∈0.5,1,2,4β∈\0.5,1,2,4\, giving the seven features ϕi=[fi,qi,ai,ci(0.5),ci(1),ci(2),ci(4)] _i=[f_i,q_i,a_i,c_i^(0.5),c_i^(1),c_i^(2),c_i^(4)]. Here, fif_i is the local interaction field, aia_i measures the total incident coupling strength, and qiq_i provides a normalized field. The four ci(β)c_i^(β) features provide nonlinear summaries of neighboring states at different coupling scales. Their functional form is inspired by the cavity equations for Ising belief propagation [49], with node states replacing directed cavity messages and an additional normalization by the incident coupling strength. Clipping keeps the argument of atanhatanh within (−1,1)(-1,1), including for Gaussian initial states. These quantities are supplied as input features; they do not define a separate belief-propagation solver or prescribe the update direction. Given these features, the recurrent cell y+=GNNθ(G,y)y^+=GNN_θ(G,y), with ϕ=ϕ(G,y)φ=φ(G,y), is defined as follows. The encoder is hi(0)=LayerNorm(We,2σ(We,1[yi,ϕi]+be,1)+be,2),h_i^(0)=LayerNorm\! (W_e,2σ(W_e,1[y_i, _i]+b_e,1)+b_e,2 ), with σ=SiLUσ=SiLU, dimensions 8→16→168→ 16→ 16, and affine LayerNorm with epsilon 10−510^-5. No zero-padding coordinates or random identifiers are supplied. The three weighted GraphConv layers are hi(ℓ+1)=σ(Wr,ℓhi(ℓ)+Wn,ℓ∑j∈N(i)Jijhj(ℓ)+bn,ℓ),ℓ=0,1,2.h_i^( +1)=σ\! (W_r, h_i^( )+W_n, _j∈ N(i)J_ijh_j^( )+b_n, ), =0,1,2. They use sum aggregation, a biased neighbor transform, and a bias-free root transform, without intermediate normalization or hidden residual connections. Independent 16→16→1616→ 16→ 16 SiLU MLPs produce ui=U(hi(3))u_i=U(h_i^(3)) and vi=R(hi(3))v_i=R(h_i^(3)). The readout combines local and graph-level information: (44) v¯=1n∑jvj,Δi=wo⊤[ui,v¯]+bo,yi+=tanh(yi+Δi). v= 1n _jv_j, _i=w_o [u_i, v]+b_o, y_i^+= (y_i+ _i). Training. For minibatch ℬB with Nℬ=∑G∈ℬ|VG|N_B= _G |V_G|, the proposed model uses the following loss in training: (1/Nℬ)∑ℬEG(yT)(1/N_B) _BE_G(y_T) with T=10T=10. Training uses such node-weighted energy within each minibatch, whereas test energies are averaged equally across graphs. In particular, we use Adam for 200 epochs, LR 10−310^-3, default betas (0.9,0.999)(0.9,0.999), epsilon 10−810^-8, and zero weight decay. There is no LR schedule, clipping, dropout, rollback, or early stopping. Graph batches are 32/32/4/1 for compact/medium/large/extra-large; each update uses one fresh trajectory per graph. The loader shuffles each epoch. The best validation checkpoint is saved after every completed epoch. Baselines. We compare single-target supervision, unique-equilibrium models, untied feedforward networks, and learned solution distributions. All methods are evaluated on the same test graphs. Table 5 summarizes model capacities and requested training budgets; resource limits and completed budgets are reported below. • Single-target supervision. We use the same architecture as our multi-attractor GNN, replacing the energy objective with MSE against one Gurobi-certified optimal spin configuration per training graph. Obtaining even one certified solution can exceed 5,000 seconds per graph at the extra-large scale, so we omit this supervised baseline at that scale. We allocate 2,000 training epochs and select the checkpoint with the lowest validation MSE. This extended budget, compared with 200 epochs for our energy-trained model, accommodates the slower optimization observed under single-target supervision: validation MSE continued to improve beyond 200 epochs. • Unique-equilibrium controls. We adopt the implicit graph neural network (IGNN) of Gu et al. [24], using its released implementation for the implicit layers and enforcing a contraction bound of ρ=0.9ρ=0.9 to ensure a unique equilibrium. To ensure fairness, we use the same physical-feature family and scalar spin-output interface as our model and train the adapters with the same energy objective. We evaluate hidden widths H∈16,32,64H∈\16,32,64\, including a capacity-comparable model and two larger variants, to assess whether increased capacity improves performance under the unique-equilibrium constraint. • Untied feedforward depth controls. We stack independently parameterized copies of our recurrent block to obtain yt+1=GNNθt(G,yt)y_t+1=GNN_ _t(G,y_t), with distinct parameters θt _t at each stage. We consider 1, 3, and 10 blocks, corresponding to 3, 9, and 30 message-passing layers. These models retain the block architecture, physical features, Gaussian node initialization, and terminal energy objective of our model. They test whether increasing feedforward depth and parameter count can match recurrent performance under the reported budgets. • Distributional and diffusion-based models. We include a mean-field GNN, an annealed Bernoulli GNN, DIFUSCO, and Variational Annealing on Graphs (VAG-CO; denoted VAG). The mean-field GNN uses our width-16, three-layer backbone to map a fresh Gaussian initialization y0∼(0,In)y_0 (0,I_n) to predicted spin means y1∈(−1,1)ny_1∈(-1,1)^n in one forward pass. It defines conditionally independent spins with Pr(si=+1∣G,y0)=(1+y1,i)/2 (s_i=+1 G,y_0)=(1+y_1,i)/2 and is trained by differentiating their analytic expected energy, [EG(s)∣G,y0]=EG(y1)E[E_G(s) G,y_0]=E_G(y_1), following the training of Karalias and Loukas [32]. The otherwise identical annealed Bernoulli GNN additionally rewards conditional Bernoulli entropy, with its coefficient decreasing linearly from 0.50.5 to zero over 200 training epochs, following the annealed variational training of Sun et al. [58]. DIFUSCO [59] is adapted to Ising using binary categorical diffusion and a denoiser conditioned on noisy spins, diffusion time, and graph features. It is trained with clean-label cross entropy against one Gurobi-certified spin configuration and generates candidates by reverse diffusion from uniform random bits. VAG [51] generates spins autoregressively in a randomized breadth-first order, jointly sampling five spins at each decoding step. It is trained with an energy-minus-entropy objective using score-function gradients and a leave-one-out baseline, with temperature annealed during training. For both DIFUSCO and VAG, the matched setting replaces their GNN backbone with our width-16, three-layer GraphConv backbone and local/global-mean readout, while retaining method-specific inputs and output heads. DIFUSCO additionally receives diffusion-time and structural features and predicts binary denoising probabilities. VAG additionally receives partial spin assignments, assignment masks, and decoding priorities, and predicts a distribution over the 252^5 assignments of each five-spin block. These additions yield 3,250 and 5,904 parameters, respectively, compared with 3,153 for our model. The large settings adopt the model-size configurations from their original GitHub repositories: width 256 with 12 layers for DIFUSCO and width 120 with three layers for VAG; see Table 5. Table 5. Learned configurations. H: hidden width; L: total learned message passing depth. Model H L Parameters LR Epochs Ours 16 3 3,153 10−310^-3 200 Single-target Supervision 16 3 3,153 10−310^-3 2000 IGNN adapter (H16) 16 3 2,851 10−310^-3 200 IGNN adapter (H32) 32 3 10,819 10−310^-3 200 IGNN adapter (H64) 64 3 42,115 10−310^-3 200 Feedforward (×1× 1) 16 3 3,153 10−310^-3 200 Feedforward (×3× 3) 16 9 9,459 10−310^-3 200 Feedforward (×10× 10) 16 30 31,530 10−310^-3 200 Mean-field GNN 16 3 3,153 10−310^-3 200 Annealed Bernoulli GNN 16 3 3,153 10−310^-3 200 DIFUSCO (matched) 16 3 3,250 10−310^-3 200 DIFUSCO (large) 256 12 1,909,762 2⋅10−42· 10^-4 100 VAG (matched) 16 3 5,904 10−310^-3 200 VAG (large) 120 3 281,192 5⋅10−45· 10^-4 30 Resource limits. Structural models, matched DIFUSCO, and energy diffusion use training batches 32/32/4/1 by size; large DIFUSCO uses 8/8/1/1, and VAG uses 32/1/1/1. Training samples for each graph are processed in parallel, but VAG tokens are sequential. Large VAG is OOM on large and extra-large graphs even alone on a 24-GB GPU at graph batch one. Larger-graph training is capped at 30 accumulated wall-clock hours, including validation and retries; capped rows use the best validated checkpoint. Matched VAG completes 68/17/3 of 200 epochs at medium/large/extra-large sizes. Metrics. We evaluate K=500K=500 test graphs, generating M=20M=20 candidate spin configurations per graph. Here, g∈1,…,Kg∈\1,…,K\ indexes test graphs and m∈1,…,Mm∈\1,…,M\ indexes candidates. Let sg,m∈−1,+1ngs_g,m∈\-1,+1\^n_g denote the mmth candidate for graph GgG_g, where ngn_g is its number of nodes. For models with continuous outputs, we set sg,m=sign+(yg,m,T)s_g,m=sign_+(y_g,m,T), where T is the evaluation horizon and sign+(0)=+1sign_+(0)=+1; discrete generators use their native spin outputs. We measure energy per node, eg,m=EGg(sg,m)/nge_g,m=E_G_g(s_g,m)/n_g. Let eg∗e_g^* be the Gurobi-certified optimal energy per node for GgG_g. The relative optimality gap is γg,m=max(0,eg,m−eg∗max(|eg∗|,ϵ)), _g,m= \! (0, e_g,m-e_g^* (|e_g^*|,ε) ), where ϵ>0ε>0 safeguards the denominator. We report three quality metrics: e¯=1KM∑g=1K∑m=1Meg,m,Best20=1K∑g=1Kmin1≤m≤Meg,m,Hit10=100KM∑g=1K∑m=1Mγg,m≤0.10. e= 1KM _g=1^K _m=1^Me_g,m,~Best_20= 1K _g=1^K _1≤ m≤ Me_g,m,~Hit10= 100KM _g=1^K _m=1^M1\ _g,m≤ 0.10\. Thus, e¯ e averages energy over all candidates, Best20Best_20 averages the best candidate energy within each graph, and Hit10Hit10 is the percentage of candidates within 10% of the certified optimum. All graphs receive equal weight. For iterative models, candidates are evaluated at the terminal step; we do not select intermediate states along a trajectory. For the experiment in Table 1, we measure diversity among optimal spin configurations, allowing a numerical tolerance of 10−610^-6: DM(g)=|sg,m:1≤m≤M,γg,m≤10−6|,GM=100K∑g=1KDM(g)≥2.D_M(g)= |\s_g,m:1≤ m≤ M,\ _g,m≤ 10^-6\ |, G_M= 100K _g=1^K1\D_M(g)≥ 2\. Here, DM(g)D_M(g) counts the distinct optimal configurations found among M candidates for graph GgG_g. We report its graph-average DM=K−1∑g=1KDM(g)D_M=K^-1 _g=1^KD_M(g). The percentage GMG_M measures how often the model finds at least two different optimal configurations for the same graph, indicating its ability to recover multiple solutions. For larger graphs, we relax the quality threshold and count distinct configurations within 10% of the certified optimum among 20 candidates: V2010%=1K∑g=1K|sg,m:1≤m≤20,γg,m≤0.10+10−6|.V_20^10\%= 1K _g=1^K |\s_g,m:1≤ m≤ 20,\ _g,m≤ 0.10+10^-6\ |. For example, V2010%=20V_20^10\%=20 means that every graph has 20 distinct candidates meeting this threshold. Repeated configurations count only once, while global spin flips are counted as distinct configurations. These metrics measure the diversity of sufficiently accurate outputs; they do not require the underlying trajectories to have converged. Recurrence is numerically converged at T when n−1‖yT−yT−1‖22<10−4 n^-1\|y_T-y_T-1\|_2^2<10^-4; CTC_T is the percentage meeting this criterion. Among nonconverged trajectories, periods p=2,…,10p=2,…,10 are tested using [1pn∑k=0p−1‖yT−k−yT−k−p‖22]1/2<10−4, [ 1pn _k=0^p-1\|y_T-k-y_T-k-p\|_2^2 ]^1/2<10^-4, assigning the smallest qualifying period. Remaining trajectories are unclassified at the tested horizon, not necessarily divergent or chaotic. Convergence is not applicable to feedforward or time-inhomogeneous generative models. Longer test-time trajectories. Table 1 reports a convergence rate slightly below 100% at t=200t=200. To determine whether the remaining trajectories exhibit persistent oscillations or simply require more iterations, we extend the same 30,000 (500×20×3500× 20× 3) trajectories across three training seeds to t=2000t=2000, without changing the learned weights. Only 11 trajectories fail the numerical fixed-point criterion at t=200t=200; all satisfy it by the measured horizon t=500t=500 and remain numerically stationary at subsequent reported horizons through t=2000t=2000 (Table 6). The extended checks also confirm sustained numerical stationarity at the final horizon, supporting the interpretation that these exceptions reflect longer transients. Although training uses only T=10T=10 iterations, the learned update remains stable and settles to numerical fixed points when applied for substantially more iterations than seen during training. These are numerical observations rather than asymptotic guarantees. The longer runs serve only to examine convergence: all reported energy, Hit10, D20D_20, and G20G_20 metrics use the terminal outputs at t=200t=200. Thus, the reported solution quality and diversity are already obtained by stopping at t=200t=200, without waiting for every trajectory to satisfy the convergence criterion. Table 6. Ising test-time dynamics. Percentages are means ± population standard deviations across three training seeds. t CtC_t (%) Period 2–10 (%) Unclassified (%) 5 3.3733±0.63703.3733± 0.6370 — — 10 32.5833±2.095632.5833± 2.0956 — — 20 77.2733±2.012977.2733± 2.0129 0.0000±0.00000.0000± 0.0000 22.7267±2.012922.7267± 2.0129 50 97.7400±0.261797.7400± 0.2617 0.0000±0.00000.0000± 0.0000 2.2600±0.26172.2600± 0.2617 100 99.7100±0.063899.7100± 0.0638 0.0000±0.00000.0000± 0.0000 0.2900±0.06380.2900± 0.0638 200 99.9633±0.017099.9633± 0.0170 0.0000±0.00000.0000± 0.0000 0.0367±0.01700.0367± 0.0170 500 100.0000±0.0000100.0000± 0.0000 0.0000±0.00000.0000± 0.0000 0.0000±0.00000.0000± 0.0000 1000 100.0000±0.0000100.0000± 0.0000 0.0000±0.00000.0000± 0.0000 0.0000±0.00000.0000± 0.0000 2000 100.0000±0.0000100.0000± 0.0000 0.0000±0.00000.0000± 0.0000 0.0000±0.00000.0000± 0.0000 Larger Graphs. Table 7 evaluates scalability on larger graphs while retaining the same width-16 architecture and 3,153 trainable parameters. Our model achieves the lowest mean and best-of-20 energies, the highest Hit10, and the largest number of distinct valid configurations among the evaluated methods at every scale. Its Hit10 reaches 96.88%96.88\%, 99.96%99.96\%, and 100%100\% on medium, large, and extra-large graphs, respectively, with an average of 19.37619.376, 19.99219.992, and 20.00020.000 distinct configurations within 10% of optimum among 20 candidates. Thus, a small shared recurrent model continues to produce both accurate and diverse solutions as graph size increases. This is consistent with the set-valued representation established by our theory: a single shared update can accommodate multiple solutions for the same input graph. Here, the empirical evidence concerns finite-horizon outputs, rather than convergence to distinct optimal attractors. Increasing baseline capacity also helps in some cases, but does not consistently recover this combination of quality and diversity. The 3-block feedforward model substantially improves over one-shot prediction, and large DIFUSCO is competitive on medium and large graphs, yet both yield worse mean and best-of-20 energies than our model. Further increasing feedforward depth to 10 untied blocks reduces Hit10 and produces the same rounded configuration across all 20 initializations for each tested graph. These results support recurrent weight sharing as an effective way to scale solution generation without increasing parameter count, while the time-capped and OOM entries identify practical resource limits of the evaluated baseline configurations. Table 7. Scaling to larger Ising graphs with fixed model capacity. The proposed model retains its width-16 architecture and 3,153 parameters across medium, large, and extra-large datasets, each containing 500 test graphs. Metrics use 20 candidates per graph: e¯ e and Best20 are mean and best-of-20 energies per node; Hit10 is the percentage of candidates within 10% of the certified optimum; and V2010%V_20^10\% is the mean number of distinct configurations meeting this threshold. † marks training capped at 30 hours; OOM denotes a memory failure. ‡ marks supervised methods omitted because of training-labeling costs. Method e¯↓ e Best20↓_20 Hit10 (%) ↑ V2010%↑V_20^10\% Medium Ours −1.2431-1.2431 −1.2862-1.2862 96.8896.88 19.37619.376 Single-target Supervision 1.52481.5248 1.29131.2913 0.000.00 0.0000.000 IGNN adapter (H16) 0.57130.5713 0.57130.5713 0.000.00 0.0000.000 IGNN adapter (H32) −0.7917-0.7917 −0.7917-0.7917 0.000.00 0.0000.000 IGNN adapter (H64) 0.00520.0052 0.00520.0052 0.000.00 0.0000.000 Feedforward (×1× 1) −1.0970-1.0970 −1.1792-1.1792 1.601.60 0.3200.320 Feedforward (×3× 3) −1.1990-1.1990 −1.2590-1.2590 58.7858.78 11.75611.756 Feedforward (×10× 10) −1.1966-1.1966 −1.1966-1.1966 56.6056.60 0.5660.566 Mean-field GNN −1.0963-1.0963 −1.1781-1.1781 1.611.61 0.3220.322 Annealed Bernoulli GNN −1.1002-1.1002 −1.1811-1.1811 1.731.73 0.3460.346 DIFUSCO (matched) −1.0864-1.0864 −1.2224-1.2224 10.4910.49 2.0982.098 DIFUSCO (large) −1.1934-1.1934 −1.2653-1.2653 55.2555.25 11.05011.050 VAG (matched)† −0.9949-0.9949 −1.0899-1.0899 0.010.01 0.0020.002 VAG (large) 0.03060.0306 −0.2003-0.2003 0.000.00 0.0000.000 Large Ours −1.2789-1.2789 −1.3023-1.3023 99.9699.96 19.99219.992 Single-target Supervision† 0.40710.4071 −0.2762-0.2762 0.000.00 0.0000.000 IGNN adapter (H16) 0.08640.0864 0.08640.0864 0.000.00 0.0000.000 IGNN adapter (H32) −0.9817-0.9817 −0.9817-0.9817 0.000.00 0.0000.000 IGNN adapter (H64) 0.12160.1216 0.12160.1216 0.000.00 0.0000.000 Feedforward (×1× 1) −1.1296-1.1296 −1.1722-1.1722 0.000.00 0.0000.000 Feedforward (×3× 3) −1.2444-1.2444 −1.2734-1.2734 75.8575.85 15.17015.170 Feedforward (×10× 10) −1.2289-1.2289 −1.2289-1.2289 41.0041.00 0.4100.410 Mean-field GNN −1.1287-1.1287 −1.1720-1.1720 0.000.00 0.0000.000 Annealed Bernoulli GNN −1.1352-1.1352 −1.1778-1.1778 0.000.00 0.0000.000 DIFUSCO (matched) −1.1285-1.1285 −1.2057-1.2057 0.590.59 0.1180.118 DIFUSCO (large) −1.2514-1.2514 −1.2876-1.2876 80.2680.26 16.05216.052 VAG (matched)† −0.9778-0.9778 −1.0318-1.0318 0.000.00 0.0000.000 VAG (large) OOM OOM OOM OOM Extra-large Ours −1.3111-1.3111 −1.3227-1.3227 100.00100.00 20.00020.000 Single-target Supervision‡ – – – – IGNN adapter (H16) 0.29370.2937 0.29370.2937 0.000.00 0.0000.000 IGNN adapter (H32) −0.8143-0.8143 −0.8143-0.8143 0.000.00 0.0000.000 IGNN adapter (H64) 2.91782.9178 2.91782.9178 0.000.00 0.0000.000 Feedforward (×1× 1) −1.1443-1.1443 −1.1665-1.1665 0.000.00 0.0000.000 Feedforward (×3× 3) −1.2753-1.2753 −1.2912-1.2912 98.8598.85 19.77019.770 Feedforward (×10× 10) −1.2390-1.2390 −1.2390-1.2390 5.005.00 0.0500.050 Mean-field GNN −1.1435-1.1435 −1.1658-1.1658 0.000.00 0.0000.000 Annealed Bernoulli GNN −1.1549-1.1549 −1.1767-1.1767 0.000.00 0.0000.000 DIFUSCO (matched)‡ – – – – DIFUSCO (large)‡ – – – – VAG (matched)† 0.09800.0980 0.03550.0355 0.000.00 0.0000.000 VAG (large) OOM OOM OOM OOM Appendix F Structural module detection in protein graphs: details and additional results This appendix gives the experiment details for Section 3.2. Dataset. We download the PROTEINS dataset through PyTorch Geometric’s TUDataset. For each graph we remove self-loops, collapse duplicate directed entries to one undirected edge, and use both directions only for message passing. The original graph classification label is retained for audit purposes but is never supplied to the models. The task graph contains 1,113 graphs, 43,471 nodes, and 81,044 undirected objective edges. We compute exact graph-isomorphism groups, randomly order the groups with split seed 2027, and allocate complete groups to 80/10/10 splits. This gives 891 training, 111 validation, and 111 test graphs. No exact topology is shared across validation and test or with another split; duplicate topologies inside training remain grouped. Feature-normalization statistics are fit on the 891 training graphs only; standardized added topology features are clipped to [−8,8][-8,8]. Table 8. Graph statistics. Entries for nodes and edges are min/mean/median/max. Split Graphs Nodes Undirected edges Train 891 4 / 36.53 / 24 / 620 5 / 67.98 / 46 / 1049 Validation 111 5 / 52.56 / 31 / 504 6 / 98.33 / 61 / 894 Test 111 8 / 45.84 / 32 / 285 16 / 86.15 / 57 / 424 For exact references, binary xikx_ik assign node i to group k, ∑kxik=1 _kx_ik=1. Auxiliary wijkw_ijk linearize xikxjkx_ikx_jk on each observed edge, and vk=∑idixikv_k= _id_ix_ik. Gurobi maximizes (45) 4m∑i,j∈E∑kwijk−∑kvk2,4m _\i,j\∈ E _kw_ijk- _kv_k^2, whose value divided by 4m24m^2 is QGQ_G. All 1,113 instances terminate with Gurobi’s optimal status. One certified partition per graph is retained as the supervised target. Node features. All neural models in this experiment receive the same node and edge features. We augment the original node attributes with standard structural descriptors drawn from local degree profiles [12], random-walk encodings [18], and graph statistics documented in NetworkX [26]. The input si∈ℝ19s_i ^19 comprises five original-attribute/degree coordinates and fourteen additional structural coordinates: si=[a~i,onehot3(ℓi),d^i,node(ψi)].s_i=[\, a_i,onehot_3( _i), d_i,Z_ node( _i)\,]. Here aia_i is the original continuous node attribute, ℓi _i is its original three-category node label, and d^i=di/max(1,maxjdj) d_i=d_i/ (1, _jd_j). The standardized attribute a~i a_i uses the training-node mean and sample standard deviation (denominator floored at 10−810^-8). The original node labels are attributes, not target module assignments. Let N(i)N(i) be the neighbor set, n=|V|n=|V|, Wij=Aij/max(di,1)W_ij=A_ij/ (d_i,1), τi _i the triangle count at i, κi _i its core number, and πi _i its PageRank. The fourteen raw structural coordinates are ψi=[ _i= [ log(1+di),(W2)ii,(W3)ii,(W4)ii,(W8)ii,ci,log(1+τi),κimax(1,maxjκj),nπi, (1+d_i),(W^2)_i,(W^3)_i,(W^4)_i,(W^8)_i,c_i, (1+ _i), _i (1, _j _j),n _i, (46) log(1+μi),log(1+σi),|N2(i)|max(n−1,1),|(i)|n,bimax(di,1)]. (1+ _i), (1+ _i), |N_2(i)| (n-1,1), |C(i)|n, b_i (d_i,1) ]. The four diagonal powers are random-walk return probabilities at lengths 2,3,4,82,3,4,8. ci=2τi/[di(di−1)]c_i=2 _i/[d_i(d_i-1)] is the local clustering coefficient, set to zero when di<2d_i<2. The neighbor-degree mean and population standard deviation are μi _i and σi _i (both zero for an empty neighborhood). N2(i)N_2(i) contains nodes reachable in two hops after excluding i and its one-hop neighbors; (i)C(i) is the connected component containing i; bib_i counts incident bridges. Core numbers use the undirected simple graph. PageRank is approximated by 80 iterations from a uniform initialization, with damping .85.85, uniform teleportation, and uniform redistribution of dangling-node mass. The operator nodeZ_ node standardizes each added coordinate using all training nodes, floors population standard deviations at 10−610^-6, and clips the result to [−8,8][-8,8]. Edge features. We combine endpoint attributes with standard neighborhood-similarity and connectivity descriptors. For a directed message edge j→ij→ i, set Cji=N(j)∩N(i)C_ji=N(j)∩ N(i) and Uji=N(j)∪N(i)U_ji=N(j)∪ N(i). The thirteen raw coordinates are χji=[ _ji= [ log(1+dj),log(1+di),d^j,d^i,[ℓj=ℓi],|a~j−a~i|,log(1+|Cji|),|Cji|max(|Uji|,1), (1+d_j), (1+d_i), d_j, d_i,1[ _j= _i],| a_j- a_i|, (1+|C_ji|), |C_ji| (|U_ji|,1), ∑w∈Cji1log(max(dw,2)+10−12),∑w∈Cji1max(dw,1),[j,i is a bridge], _w∈ C_ji 1 ( (d_w,2)+10^-12), _w∈ C_ji 1 (d_w,1),1[\j,i\ is a bridge], OPEN1max(djdi,1),|Cji|max(min(dj,di),1)]. 1 (d_jd_i,1), |C_ji| ( (d_j,d_i),1) ]. The neighborhood-similarity terms include common-neighbor counts and Jaccard similarity [37], the Adamic–Adar index [2], and the resource-allocation and hub-promoted indices [70]. The last coordinate is the hub-promoted index. Standard implementations of Jaccard, Adamic–Adar, and resource-allocation scores are also available in NetworkX [26]. We use eji=edge(χji)e_ji=Z_ edge( _ji), with training-edge population moments, standard-deviation floor 10−610^-6, and clipping to [−8,8][-8,8]. Both edge orientations are included in message passing. All added features are deterministic functions of graph topology and original attributes. They contain no node identifiers, target partitions, or reference optimizer outputs, and do not depend on the current assignment state. State and initialization. The recurrent variable is the soft partition Pt∈ℝn×4P_t ^n× 4. For each trajectory we independently draw node-wise and graph-wise Gaussian logits and initialize (47) P0,i=softmax(ϵi+bG2),ϵi,bG∼i.i.d.(0,I4).P_0,i=softmax\! ( _i+b_G 2 ), _i,b_G i.i.d. N(0,I_4). The same bGb_G is shared by the nodes of a graph within a trajectory, but is redrawn for each trajectory. Message-passing cell. We use a width-32, three-layer PNA-style cell [16]. In particular, let σ=ReLUσ=ReLU, at each time step the hidden embeddings are recomputed as hi(0)=LayerNorm(We,2σ(We,1[si,Pt,i]+be,1)+be,2),h_i^(0)=LayerNorm (W_e,2σ(W_e,1[s_i,P_t,i]+b_e,1)+b_e,2 ), with encoder dimensions 23→32→3223→ 32→ 32. The state-pair edge vector is (48) rji(Pt)=[Pt,j,Pt,i,Pt,j⊙Pt,i,|Pt,j−Pt,i|]∈ℝ16.r_ji(P_t)=[P_t,j,P_t,i,P_t,j P_t,i,|P_t,j-P_t,i|] ^16. For ℓ=0,1,2 =0,1,2, compute mji(ℓ) m_ji^( ) =Mℓ([hj(ℓ),hi(ℓ),eji,rji(Pt)]), =M_ ([h_j^( ),h_i^( ),e_ji,r_ji(P_t)]), ai(ℓ) a_i^( ) =[meanj∈N(i)mji(ℓ),stdj∈N(i)mji(ℓ),maxj∈N(i)mji(ℓ),sumj∈N(i)mji(ℓ)], =[\,mean_j∈ N(i)m_ji^( ),std_j∈ N(i)m_ji^( ),max_j∈ N(i)m_ji^( ),sum_j∈ N(i)m_ji^( )\,], (49) hi(ℓ+1) h_i^( +1) =LayerNorm(hi(ℓ)+Uℓ([hi(ℓ),ai(ℓ),Pt,i])). =LayerNorm (h_i^( )+U_ ([h_i^( ),a_i^( ),P_t,i]) ). MℓM_ is a 93→32→3293→ 32→ 32 MLP with ReLU after both affine maps; UℓU_ is a 164→32→32164→ 32→ 32 MLP with ReLU only between affine maps. Aggregations are elementwise. Global readout and update. The graph bottleneck and correction head are (50) gG=1n∑iR(hi(3))∈ℝ4,Δθ,i=O([hi(3),gG]),g_G= 1n _iR(h_i^(3)) ^4, _θ,i=O([h_i^(3),g_G]), where R is a 32→32→432→ 32→ 4 MLP and O is a 36→32→436→ 32→ 4 MLP, each with ReLU between its two affine maps. We update (51) Pt+1=softmax[log(Pt∨10−6)+Δθ(G,Pt)],P_t+1=softmax\! [ (P_t 10^-6)+ _θ(G,P_t) ], with row-wise softmax and coordinate-wise clipping before the logarithm. The three MP layers have distinct parameters within a cell; the entire cell is tied across time. There are 35,784 trainable parameters: 1,888 in the encoder, 3×10,4643× 10,464 in message passing, 1,188 in the global map, and 1,316 in the output head. The full model sizes are provided in Table 9. Training and objectives. For an undirected graph G=(V,E)G=(V,E) with m=|E|m=|E| edges and node degrees did_i, let zi∈1,…,Kz_i∈\1,…,K\ denote the module of node i. At resolution γ=1γ=1, we minimize the negative modularity: EG(z)=−QG(z)=∑k=1K(∑idi[zi=k]2m)2−1m∑i,j∈E[zi=zj].E_G(z)=-Q_G(z)= _k=1^K ( _id_i1[z_i=k]2m )^2- 1m _\i,j\∈ E1[z_i=z_j]. The edge term rewards placing connected nodes in the same module; alone, it is minimized by assigning all nodes to one module. The degree-volume term accounts for the within-module connectivity expected from node degrees under the modularity null model. Combining the two terms therefore favors modules with more internal connectivity than this baseline predicts. For a soft partition P, where each PiP_i is a probability vector over the K modules, we define the differentiable energy (52) eG(P)=−1m∑i,j∈E⟨Pi,Pj⟩+14m2[‖∑idiPi‖22−∑idi2∥Pi∥22].e_G(P)=- 1m _\i,j\∈ E P_i,P_j + 14m^2 [ \| _id_iP_i \|_2^2- _id_i^2\|P_i\|_2^2 ]. Specifically, for independent assignments Zi∼Categorical(Pi)Z_i (P_i), eG(P)=[EG(Z)]−14m2∑idi2.e_G(P)=E[E_G(Z)]- 14m^2 _id_i^2. Thus, the soft objective equals the expected discrete energy up to a graph-dependent constant. During training, each update samples ten training graphs with replacement and four independent initial soft partitions P0P_0 per graph. We unroll T=5T=5 tied cells and minimize [eG(PT)]E[e_G(P_T)]. We use Adam with learning rate 3×10−43× 10^-4, β1=.9 _1=.9, β2=.999 _2=.999, ϵ=10−8ε=10^-8, zero weight decay, global gradient-norm clipping at 1, 3,500 optimizer updates, batch size ten base graphs, and no learning-rate decay. The final update-3,500 checkpoint is used; there is no best-test or best-validation checkpoint substitution. Our model and all baselines in Table 2 use training seeds 0, 1, and 2. The unsupervised model never receives the Gurobi labels. Baselines. Table 9 summarizes the learned configurations used in the reported comparisons. In this PROTEINS experiment, we use the same baseline groups as in the Ising experiment: • Single-target supervision. We use the same architecture as our multi-attractor GNN, replacing the energy objective eG(P)e_G(P) with cross entropy against one Gurobi-optimal partition per graph. • Unique-equilibrium models. As in the Ising experiments, we adopt the IGNN with hidden widths H∈16,32,64H∈\16,32,64\ and enforce a contraction bound of ρ=0.9ρ=0.9. To train IGNN, we use the same unsupervised loss function [eG(PT)]E[e_G(P_T)]. Again, the input features are the same as in the multi-attractor setting. • Untied feedforward depth controls. The feedforward controls replace the same three-layer rich-PNA cell by B∈1,2,3,4,5B∈\1,2,3,4,5\ independent copies. They receive the same random P0P_0, losses, four training trajectories, optimizer, and 3,500-update budget. • Distributional and diffusion-based models. We include the same distributional models as in the Ising experiments: a mean-field GNN, an annealed Bernoulli GNN, DIFUSCO, and Variational Annealing on Graphs (VAG-CO; denoted VAG). They are still applicable here. “Matched” setting uses the same model backbone, learning rate, and training updates. “Large” setting uses the upstream original repository’s configurations. Details are in Table 9. Table 9. Learned configurations on PROTEINS. H: hidden width; L: total untied depth. Model H L Parameters LR Updates Ours 32 3 35,784 3⋅10−43· 10^-4 3,500 Single-target Supervision 32 3 35,784 3⋅10−43· 10^-4 3,500 IGNN (H16) 16 3 11,340 3⋅10−43· 10^-4 3,500 IGNN (H32) 32 3 39,564 3⋅10−43· 10^-4 3,500 IGNN (H64) 64 3 146,700 3⋅10−43· 10^-4 3,500 Feedforward (×1× 1) 32 3 35,784 3⋅10−43· 10^-4 3,500 Feedforward (×2× 2) 32 6 71,568 3⋅10−43· 10^-4 3,500 Feedforward (×3× 3) 32 9 107,352 3⋅10−43· 10^-4 3,500 Feedforward (×4× 4) 32 12 143,136 3⋅10−43· 10^-4 3,500 Feedforward (×5× 5) 32 15 178,920 3⋅10−43· 10^-4 3,500 Mean-field GNN 32 3 35,784 3⋅10−43· 10^-4 3,500 Annealed Bernoulli GNN 32 3 35,784 3⋅10−43· 10^-4 3,500 DIFUSCO (matched) 32 3 35,880 3⋅10−43· 10^-4 3,500 DIFUSCO (large) 256 12 5,555,466 2⋅10−42· 10^-4 2,800 VAG (matched) 32 3 40,270 3⋅10−43· 10^-4 3,500 VAG (large) 40 3 101,790 5⋅10−45· 10^-4 10,000 Metrics. The caption of Table 2 defines every displayed metric. Longer test-time trajectories. Table 2 reports a one-step convergence rate of C200=88.05±2.23%C_200=88.05± 2.23\% at t=200t=200. To determine whether the remaining trajectories exhibit persistent oscillations or simply require more iterations, we extend the same 6,660 (111 × 20 × 3) trajectories across three training seeds to t=10000t=10000, without changing the learned weights. Here t denotes test-time recurrence steps, and results are reported in Table 10. The convergence rate increases to C800=98.75±0.62%C_800=98.75± 0.62\% and C3000=99.49±0.29%C_3000=99.49± 0.29\%. By t=5000t=5000, every trajectory is either convergent or classified as periodic at the reported horizons, with 6,632 trajectories satisfying the numerical convergence criterion. The extended checks confirm that the vast majority of non-converged trajectories at t=200t=200 simply reflect longer transients, while a small remainder settles into stable periodicity rather than divergent behavior. Although training uses only T=5T=5 unrolled iterations, the learned update remains stable and largely settles to numerical fixed points when applied for substantially more iterations than seen during training. The longer runs serve only to examine convergence: all reported modularity, Hit10, and discovery metrics use the terminal outputs at t=200t=200. Thus, the reported solution quality and diversity are already obtained by stopping at t=200t=200, without waiting for every trajectory to strictly satisfy the convergence criterion. Table 10. PROTEINS test-time dynamics. Metrics are the same as in Appendix E, Table 6. t CtC_t (%) Period 2-10 (%) Unclassified (%) 10 3.83±1.783.83± 1.78 – – 20 16.65±3.4516.65± 3.45 – – 50 47.18±1.9247.18± 1.92 – – 100 68.72±0.8268.72± 0.82 – – 200 88.05±2.2388.05± 2.23 0.02±0.020.02± 0.02 11.94±2.2511.94± 2.25 500 97.51±0.4597.51± 0.45 0.39±0.280.39± 0.28 2.10±0.342.10± 0.34 800 98.75±0.6298.75± 0.62 0.42±0.280.42± 0.28 0.83±0.480.83± 0.48 1000 99.07±0.3099.07± 0.30 0.42±0.280.42± 0.28 0.51±0.240.51± 0.24 2000 99.49±0.3199.49± 0.31 0.42±0.280.42± 0.28 0.09±0.070.09± 0.07 3000 99.49±0.2999.49± 0.29 0.42±0.280.42± 0.28 0.09±0.040.09± 0.04 5000 99.58±0.2899.58± 0.28 0.42±0.280.42± 0.28 0.00±0.000.00± 0.00 10000 99.58±0.2899.58± 0.28 0.42±0.280.42± 0.28 0.00±0.000.00± 0.00 Appendix G Chemical reaction networks: details and additional results This appendix provides the experiment details and additional results for Section 3.3. Dataset. We construct kinetic instances on the published reaction-network structures from [68]. We select networks labeled Multi from the Atom and Joshi construction families, preserving their original train/validation/test assignments. These labels indicate structural capacity for multistationarity but do not provide kinetic rates or stationary concentrations. For each selected structure, we add one inflow ∅→Xi → X_i and one outflow Xi→∅X_i→ per species. We then randomly propose two positive concentration vectors c(1),c(2)c^(1),c^(2) in log space and solve a linear feasibility problem for shared positive rate constants κ satisfying fG(c(1))=fG(c(2))=0f_G(c^(1))=f_G(c^(2))=0; the equations are linear in κ when the concentrations are fixed. We retain instances whose two witnesses have normalized residuals at most 10−710^-7, concentration-Jacobian condition numbers at most 10910^9, and log-RMS separation at least 0.250.25. After feasibility filtering, we balance the two construction families within each split, obtaining 1,242 training, 266 validation, and 236 test instances, with 2–5 species and 6–19 reactions. Each instance contains the network structure, constructed rate constants, and two numerically verified stationary witnesses, which need not exhaust its solution set. Physics-only training uses only (α,β,κ)(α,β,κ); single-target supervision and supervised diffusion use one witness selected once and held fixed as the label for each training instance. Table 11. Dataset construction. Balance is enforced within each split after feasibility filtering. Every retained instance has two constructed witnesses. Split Atom attempts Joshi attempts Atom feasible Joshi feasible Retained Reactions Train 790 2560 771 621 1242 6–19 Validation 180 570 172 133 266 6–18 Test 300 600 293 118 236 6–18 Nodes and edges. There is one node per species and one per reaction, including inflow/outflow reactions. For every pair (i,r)(i,r) with αri>0 _ri>0 or βri>0 _ri>0, we include both directed incidence edges. There is no extra node for ∅ , no stoichiometric-copy expansion, and no direct species–species edge. The six static node coordinates are xi x_i =(1,0,0,0,0,0), =(1,0,0,0,0,0), species, , (53) xr x_r =(0,1,logκr,∑iαri,∑iβri,∑i|βri−αri|), = (0,1, _r, _i _ri, _i _ri, _i| _ri- _ri| ), reaction. . The first two coordinates distinguish node types; the remaining coordinates encode the log rate, total reactant/product orders, and net stoichiometric magnitude. Apart from the log-rate transform, these graph inputs are used directly, without dataset-fitted standardization. The four edge coordinates are ei→r=(αri,βri,1,0),er→i=(αri,βri,0,1).e_i→ r=( _ri, _ri,1,0), e_r→ i=( _ri, _ri,0,1). The last two entries encode message direction, not reactant versus product. A species can be both a reactant and a product, with its two coefficients retained on the same incidence pair. Reaction roles are thus edge attributes, not a single global species label. All compared neural models receive these same base graph features. State and initialization. The recurrent state is only the species vector zt=logct∈[−20,20]sz_t= c_t∈[-20,20]^s. With width w=32w=32, compute the static encoding bv=tanh(LayerNorm(Wxxv+bx)).b_v= \! (LayerNorm(W_xx_v+b_x) ). Each candidate draws independent node vectors ξv∼(0,I8) _v (0,I_8). A bias-free 8→328→ 32 map and a scalar affine head initialize z0,i=clip[−20,20](winit(bi+Wξξi)+binit).z_0,i=clip_[-20,20] (w_ init T(b_i+W_ξ _i)+b_ init ). Subsequent calls receive no identifier forcing. At each call, reconstruct a temporary node matrix H(0)H^(0) from b, replace coordinate zero of each species row with its current ziz_i, and leave reaction rows equal to brb_r. All temporary hidden coordinates are discarded after the call except the new species log concentrations. There is no auxiliary carried hidden state. Message-passing backbone. One update contains three distinct learned GATv2 [11] blocks; the entire stack is tied across recurrent iterations. For ℓ=1,2,3 =1,2,3, compute U(ℓ)=LayerNormℓ(GATv2ℓ(H(ℓ−1),E,e)),H(ℓ)=tanhU(ℓ),ℓ<3,U(ℓ),ℓ=3.U^( )=LayerNorm_ (GATv2_ (H^( -1),E,e) ), H^( )= cases U^( ),& <3,\\ U^( ),& =3. cases Each block is PyG GATv2Conv(32,32,heads=4,concat=False,edge_dim=4). Attention heads are averaged, not concatenated. Self-loops are enabled, with mean incident edge features for their attributes; these internal loops are not chemical reactions. Other settings are zero attention dropout, LeakyReLU slope 0.20.2, separate source/target transforms, and biases enabled. LayerNorm is affine with epsilon 10−510^-5. There is no spectral constraint on these attention blocks. Dynamic physical features. For c=expzc= z, compute ℓr=logκr+∑iαrizi _r= _r+ _i _riz_i and protected fluxes v~r=exp(clip[−30,30](ℓr)) v_r= (clip_[-30,30]( _r)). Set a~=|N|v~ a=|N| v and q~i=(Nv~)i/(a~i+10−12) q_i=(N v)_i/( a_i+10^-12). The node-wise dynamic features are ψi=(zi,q~i,log(a~i+10−12)),ψr=(ℓr,0,0). _i=(z_i, q_i, ( a_i+10^-12)), _r=( _r,0,0). A shared affine 3→323→ 32 map followed by tanh encodes ψv _v. The reaction log-flux feature itself is not clipped; clipping protects the exponentials used to construct the internal residual. The independent evaluation residual uses original coefficients and unclipped mass-action fluxes. Gated readout and update. A bias-free conditioning map WbW_b and the physical-feature encoder give H^=GRUCell(tanh[H(3)+Wbb+tanh(Wψψ+bψ)],H(0)). H=GRUCell\! ( \! [H^(3)+W_bb+ (W_ψ+b_ψ) ],H^(0) ). This width-32 GRU [15] is a within-call gating computation, not an additional inter-iteration memory. Three scalar affine heads applied to each species row of H H produce di=12σ(wdH^i+bd),gi=2tanh(wgH^i+bg),aiNN=tanh(waH^i+ba).d_i= 12σ(w_d T H_i+b_d), g_i=2 (w_g T H_i+b_g), a_i N= (w_a T H_i+b_a). The update is (54) zi+=clip[−20,20][zi+di(giq~i+tanh(R~)aiNN)],R~=s−1∑iq~i2+10−24.z_i^+=clip_[-20,20] [z_i+d_i (g_i q_i+ ( R)a_i N ) ], R= s^-1 _i q_i^2+10^-24. The signed-gain head starts with zero weights and bias atanh(1/2)atanh(1/2), so its initial gain is one. The model has 35,716 trainable parameters. It is a physics-assisted learned solver, not a generic GNN with only a physics loss, and is not guaranteed contractive. All aggregation/readout operations are permutation-equivariant. Initialization is equivariant when ξ is permuted with the nodes, and equivariant in distribution under independent identifier sampling. After initialization, the shared map depends only on G and z. Training and residual objective. Using the original mass-action fluxes, define (55) qi(c)=(Nv(c))i(|N|v(c))i+10−12,R(c)=1s∑iqi(c)2,Eres(c)=1s∑iϕ0.1(qi(c)),q_i(c)= (Nv(c))_i(|N|v(c))_i+10^-12,~~R(c)= 1s _iq_i(c)^2,~~E_ res(c)= 1s _i _0.1(q_i(c)), where ϕ0.1(u)=u2/0.2 _0.1(u)=u^2/0.2 when |u|<0.1|u|<0.1, and |u|−0.05|u|-0.05 otherwise. We train the model by minimizing [Eres(cT)]E[E_ res(c_T)]. In particular, we use AdamW with LR 0.0020.002, weight decay 10−510^-5, default (β1,β2)=(0.9,0.999)( _1, _2)=(0.9,0.999), epsilon 10−810^-8, gradient-norm clipping at 5, at most 50 epochs, and early-stopping patience 12. Accumulate gradients over four graphs, with 20 candidate initializations per graph. Graph’s recurrent unroll length T is sampled uniformly from 4,8,16,24,32\4,8,16,24,32\. Baselines. We retain the four baseline groups used in the Ising and PROTEINS experiments: single-target supervision, unique-equilibrium models, untied feedforward depth controls, and distributional models. Table 12 summarizes their configurations. All models receive the same base graph inputs xi,xr,ei→r,er→ix_i,x_r,e_i→ r,e_r→ i and are evaluated on the same 236 test graphs. • Single-target supervision. We use the same architecture as our multi-attractor model, with the same random initialization, randomized unroll lengths, optimizer, and stopping rule, but replace the residual objective by the concentration-space MSE s−1‖cT−ctarget‖22s^-1\|c_T-c target\|_2^2. For each training instance, ctargetc target is drawn uniformly from the two witnesses once and then held fixed throughout training. The supervised diffusion models below use the same label. • Unique-equilibrium controls. As in the other two experiments, we adopt IGNN [24] and enforce a contraction bound ρ=0.9ρ=0.9, which guarantees a unique equilibrium. It is trained with the same residual objective [Eres]E[E_ res] as our model. We evaluate hidden widths H∈16,32,64H∈\16,32,64\, which include a capacity-comparable model (H=32H=32) and a smaller and a larger variant. Because its equilibrium is uniquely determined by G, IGNN produces a single candidate per graph; replicated copies are not counted toward diversity. • Untied feedforward depth controls. We replace the tied recurrence by B∈1,2,3,4,5B∈\1,2,3,4,5\ independently parameterized copies of our update, zt+1=GNNθt(G,zt)z_t+1=GNN_ _t(G,z_t), with distinct parameters θt _t at each stage. Each stage has its own input encoder, three GATv2 blocks, GRU, physical-feature encoder, and scalar heads. Only the identifier projection and initialization head, which produce z0z_0, are shared. The model thus has 3B3B message-passing blocks and 289+35,427B289+35,427B parameters, and it executes exactly B stages at both training and test time. It uses the same random initialization, 20 candidates per graph, residual objective, and optimizer as our model. These controls compare increasing untied depth and parameter count with shared recurrence; inference computation is not matched. None is iterated beyond its trained depth, so no settling rate is reported. • Distributional models. Since concentrations are continuous, we replace the discrete generators used for Ising and PROTEINS by DDPM [28], DDIM [57], and the time-reversed diffusion sampler (DIS) [8]. DDPM is trained with the standard noise-prediction loss against the same fixed witness used for single-target supervision, after mapping log concentrations to the unbounded latent u=20atanh(logc/20)u=20\,atanh( c/20). It uses 1,000 linear-beta noise levels and 100 reverse steps. DDIM shares the trained DDPM checkpoint and samples deterministically (η=0η=0) with 100 steps. DIS uses no witness. Instead, it learns to sample from the residual-based density ∝exp[−Eres(c)/τ] [-E_ res(c)/τ] with τ=10−4τ=10^-4, using 100 steps. The matched setting keeps our width-32, three-block GATv2 backbone, GRU, physical-feature encoder, and heads. The only change is that the input encoder also receives the noisy species state and the normalized diffusion time, giving 35,780 parameters compared with 35,716 for our model. It uses the same learning rate, 50-epoch budget, and early-stopping patience as our model. The large setting follows the upstream repositories. DDPM/DDIM use DIFUSCO’s gated residual graph backbone with width 256 and 12 layers (5,205,249 parameters). DIS uses the original FourierMLP controller, conditioned on a width-32, three-block GAT encoder (59,330 parameters), and is trained with its native log-variance loss. Every generator returns 20 candidates per graph after 100 sampling steps. Table 12. Learned configurations for chemistry. H: hidden width; L: total untied depth. Model H L Parameters LR Epochs Ours 32 3 35,716 2⋅10−32· 10^-3 50 Single-target Supervision 32 3 35,716 2⋅10−32· 10^-3 50 IGNN (H16) 16 3 10,273 2⋅10−32· 10^-3 50 IGNN (H32) 32 3 37,441 2⋅10−32· 10^-3 50 IGNN (H64) 64 3 142,465 2⋅10−32· 10^-3 50 Feedforward (×1× 1) 32 3 35,716 2⋅10−32· 10^-3 50 Feedforward (×2× 2) 32 6 71,143 2⋅10−32· 10^-3 50 Feedforward (×3× 3) 32 9 106,570 2⋅10−32· 10^-3 50 Feedforward (×4× 4) 32 12 141,997 2⋅10−32· 10^-3 50 Feedforward (×5× 5) 32 15 177,424 2⋅10−32· 10^-3 50 DDPM / matched 32 3 35,780 2⋅10−32· 10^-3 50 DDIM / matched 32 3 35,780 2⋅10−32· 10^-3 50 DIS / matched 32 3 35,780 2⋅10−32· 10^-3 50 DDPM / large 256 12 5,205,249 2⋅10−42· 10^-4 200 DDIM / large 256 12 5,205,249 2⋅10−42· 10^-4 200 DIS / large 32/64 3/4 59,330 5⋅10−35· 10^-3 100 Metrics. All metrics use the n=236n=236 test graphs with equal weight per graph. Graph g has MgM_g candidates cgmc_gm, with Mg=20M_g=20. A candidate is accepted, Igm=1I_gm=1, if it is finite, all its coordinates exceed 10−1210^-12, and R(cgm)≤10−3R(c_gm)≤ 10^-3. Here R and EresE_ res are defined as in the training objective, but are recomputed in float64 from the original α,β,κα,β,κ rather than from the clipped internal features of the cell. We report E¯res=1n∑g1Mg∑mEres(cgm),Hit10−3=100n∑g1Mg∑mIgm. E_ res= 1n _g 1M_g _mE_ res(c_gm), _10^-3= 100n _g 1M_g _mI_gm. Nonfinite residuals are kept in these means, and acceptance involves no polishing or matching to reference roots. To count distinct predictions, we process the accepted candidates of each graph in saved order. Each one joins the first existing representative within distance ‖logc−logc′‖2/sg≤10−3\| c- c \|_2/ s_g≤ 10^-3, or otherwise becomes a new representative. With KgK_g the resulting number of representatives (Kg=0K_g=0 if no candidate is accepted), we report (56) D20raw=1n∑gKg,G20raw=100n∑gKg≥2.D_20 raw= 1n _gK_g, G_20 raw= 100n _g 1\K_g≥ 2\. Numerically separated accepted predictions can approximate the same root; these counts therefore do not certify distinct stationary solutions. For recurrent models, we measure numerical settling at horizon T by (57) CT=100n∑g1Mg∑m‖zT+1(g,m)−zT(g,m)‖2sg<10−3,C_T= 100n _g 1M_g _m 1\! \ \|z_T+1^(g,m)-z_T^(g,m)\|_2 s_g<10^-3 \, and the column C in Table 3 reports C200C_200. This one-step test on the log-state is independent of the physical residual and does not certify asymptotic convergence. For IGNN, C is instead the percentage of graphs for which all three implicit layers satisfy ‖Fℓ,g(Hℓ,200)−Hℓ,200‖∞,p<10−3\|F_ ,g(H_ ,200)-H_ ,200\|_∞,p<10^-3 after 200 updates per layer. This measures hidden-state settling, not chemical validity. Feedforward and generative models have no autonomous state update, so C is not reported for them. Longer test-time dynamics. To investigate the trajectories that fail the T=200T=200 one-step criterion, we extend the same three frozen multi-attractor models and initialization seeds to t=20,000t=20,000. No model weights, initialization distributions, or physical acceptance thresholds change. The main table’s quality and diversity remain evaluated at T=200T=200, not at a retrospectively chosen longer horizon. Table 13. Chemistry test-time dynamics, mean ± population SD across three training seeds. CtC_t is the fraction with next-step species log-state RMS change below 10−310^-3. Oscillatory denotes detected approximate periods 22–6464 among trajectories failing that criterion; unresolved is the remainder. All percentages use the full trajectory population, and the three categories partition that population for t≥200t≥ 200, up to rounding. Dashes denote unperformed early-horizon oscillation tests. t CtC_t (%) Oscillatory (%) Unresolved (%) 10 2.26 ± 0.65 — — 20 13.77 ± 3.72 — — 50 63.11 ± 2.05 — — 100 83.99 ± 0.79 — — 200 91.99 ± 0.35 1.75 ± 1.50 6.26 ± 1.27 500 95.20 ± 1.72 2.23 ± 1.70 2.56 ± 0.54 1000 95.64 ± 1.67 2.49 ± 2.04 1.88 ± 0.46 2000 95.95 ± 1.83 2.77 ± 2.09 1.28 ± 0.26 5000 96.15 ± 2.01 2.77 ± 2.09 1.08 ± 0.27 10000 96.13 ± 2.01 2.77 ± 2.09 1.10 ± 0.27 20000 96.12 ± 2.01 2.77 ± 2.09 1.11 ± 0.31 Oscillation test: For t≥200t≥ 200, use the 257257 saved log states wj=zt+jw_j=z_t+j, j=0,…,256j=0,…,256. Let A=maxj‖wj−w¯‖2/sA= _j\|w_j- w\|_2/ s. Among trajectories failing the one-step settling criterion, select the smallest p∈2,…,64p∈\2,…,64\ satisfying max0≤j≤256−p‖wj+p−wj‖2s≤min(10−3,0.02A). _0≤ j≤ 256-p \|w_j+p-w_j\|_2 s≤ (10^-3,0.02A). Additionally require A>0A>0, half-window amplitude ratio in [0.8,1.25][0.8,1.25], and RMS distance at most 10−310^-3 between the first/last cycle centers. Half-window amplitudes are maximum RMS deviations from their respective means over w0,…,w127w_0,…,w_127 and w128,…,w256w_128,…,w_256; cycle centers average the first and last p states. The window contains at least four cycles for every searched period. Findings: Settling increases from 91.99±0.35%91.99± 0.35\% at t=200t=200 to 96.15±2.01%96.15± 2.01\% at t=5,000t=5,000 and then plateaus. At t=20,000t=20,000, 96.12±2.01%96.12± 2.01\% are step-settled, 2.77±2.09%2.77± 2.09\% are oscillatory, and 1.11±0.31%1.11± 0.31\% are unresolved. The unresolved group contains 157 trajectories across the three seeds. The oscillatory group includes 40 small-amplitude trajectories excluded by the earlier amplitude-gated diagnostic. No nonfinite states occurred. Settling need not increase monotonically because each horizon is tested independently. The log-state box precludes finite-state norm divergence but not bounded nonconvergence. Unresolved does not establish slow convergence, chaos, or divergence: the category can include drift, irregular fluctuations, periods beyond 64, or transients. These observations concern the learned solver, not physical chemical kinetics, and do not establish asymptotic attraction. Appendix H Related Work GNNs and their theory. Message-passing GNNs are at most as expressive as the 1-WL test [66, 43], and a line of work characterizes which invariant and equivariant maps GNNs can approximate [41, 14, 33, 4, 21]. Unique node identifiers or i.i.d. random node features restore universality with high probability [40, 52, 1]. Recent work extends this to measurable targets with explicit approximation rates [23], and to recurrent GNNs, which with random initialization can uniformly express any polynomial-time computable function on connected graphs [50]. All of these results concern single-valued targets G↦f(G)G f(G), with randomness serving only to tell nodes apart. In our setting, the same random initial state must do double duty: it identifies nodes and selects the solution branch. Only the output-dimensional state yty_t is carried between steps, so node distinguishability must be preserved by dynamics that also converge to a solution (Lemma 6, Theorem 2). Recurrent/equilibrium GNNs. The original GNN was defined as the fixed point of a contraction [53]. Later weight-tied models unroll a shared update [35, 17], while several implicit GNNs impose sufficient conditions for equilibrium existence and uniqueness, including contraction and strong monotonicity [24, 39, 47, 6]. Yang et al. [67] also analyze unfolded GNN regimes with stationarity guarantees that do not require uniqueness. Deep equilibrium models compute fixed points [5], with monotone variants guaranteeing uniqueness [63]; path independence, meaning convergence to the same state across initializations, has also been advocated [3]. Weight-tied message passing is also used for algorithmic reasoning and size extrapolation [62, 60, 7] and for SAT and constraint satisfaction [55, 61]. RUN-CSP in particular starts from random states and keeps the best of several parallel runs [61], but without a representational account of the solutions reached. On the theory side, Liu et al. [38] show that regular implicit operators (contractive in the state, Lipschitz in the input) represent exactly the locally Lipschitz maps, so that iterations rather than parameters supply expressivity. Jore and Liu [30] represent set-valued maps through attractor landscapes and demonstrate multi-solution discovery with recurrent GNNs on Ising problems, but do not establish permutation-equivariant dynamics or their message-passing realization. We develop this theory on graphs, where symmetry changes the picture in three ways: (i) the solution set is equivariant although individual branches need not be (Section 2.1); (i) our message-passing construction realizes the dynamics through symmetry breaking carried by the state itself, without auxiliary inputs (Theorem 2, Lemma 6); and (i) continuous, equivariant, globally contractive updates can miss every valid solution with positive probability (Theorems 3 and 6). As in Liu et al. [38], complexity comes from iteration: each step of the ideal update is globally Lipschitz, while its limit selector jumps across basin boundaries (Theorem 1). Multiple attractors have long served as a resource in Hopfield networks [29, 48], where they store fixed patterns rather than solutions of an input-dependent problem. Distributional models. Another way to handle multiple solutions is to learn a distribution over them. Unsupervised neural combinatorial optimization trains GNNs from an energy, either as mean-field relaxations [32, 54] or as annealed, autoregressive, diffusion, or GFlowNet samplers [58, 51, 69]; variational autoregressive networks play the same role for statistical-mechanics models [64]. Supervised diffusion solvers and diffusion samplers generate candidates through time-inhomogeneous reverse processes [59, 28, 57, 8]. Other work keeps supervision but relaxes the single target: multiple output heads with a hindsight loss [36], learned selection among one-of-many solutions [45], label realignment under formulation symmetry [13], and equivariant flow matching with symmetry-aware target matching for bifurcation problems [27]. Our approach learns a single autonomous update whose attracting limits are intended to approximate valid solutions. Numerical stationarity and task validity are evaluated separately; iterations can be extended at test time without retraining, and training uses energies or residuals rather than solution labels. Classical multi-solution methods such as deflation and homotopy continuation [20, 10] operate instance by instance, whereas our model is trained once and amortized over a graph family. Symmetry breaking. Deterministic equivariant networks cannot map a symmetric input to a less symmetric output [56, 31]. Existing remedies break symmetry through augmented inputs: random or positional node features [52, 1, 18], symmetry-breaking sets and objects [56, 65, 22], and probabilistic symmetry breaking through equivariant conditional distributions and sampled canonicalizations [34], the notion closest to our set-valued equivariance. We instead break symmetry through the model’s own evolving state, with no auxiliary input: randomness enters once through y0y_0, and Lemma 6 shows that the ideal dynamics preserve node distinguishability at every finite step. The two mechanisms are compatible. Augmented features can be appended to X or used to seed y0y_0 (our PROTEINS model uses random-walk encodings, and our chemistry model seeds y0y_0 from random node vectors), and the construction in Lemma 5 only needs distinct node codes, whatever their source. Augmentation alone, however, does not settle which solution the network should output: deterministic encodings alone do not provide a mechanism for sampling multiple valid outputs, and continuous canonicalization is impossible in general [19]. Theorems 3 and 6 concern updates without auxiliary inputs; combining multi-attractor dynamics with augmented features is a natural direction for future work.