Paper deep dive
Dynamical Systems Theory Behind a Hierarchical Reasoning Model
Vasiliy A. Es'kin, Mikhail E. Smorkalov
Intelligence
Status: succeeded | Model: google/gemini-3.1-flash-lite-preview | Prompt: intel-v1 | Confidence: 96%
Last extracted: 3/26/2026, 1:40:53 AM
Summary
The paper introduces the Contraction Mapping Model (CMM), a novel architecture for algorithmic reasoning that reformulates discrete recursive reasoning into continuous Neural Ordinary and Stochastic Differential Equations (NODEs/NSDEs). By enforcing convergence to stable equilibrium states and utilizing a hyperspherical repulsion loss to prevent feature collapse, the CMM achieves state-of-the-art performance on benchmarks like Sudoku-Extreme and Maze with extreme parameter efficiency (as low as 0.26M parameters), outperforming larger Hierarchical Reasoning Models (HRM) and Tiny Recursive Models (TRM).
Entities (5)
Relation Signals (3)
Contraction Mapping Model → evaluatedon → Sudoku-Extreme
confidence 95% · On the Sudoku-Extreme benchmark, a 5M-parameter CMM achieves a state-of-the-art accuracy of 93.7 %
Contraction Mapping Model → outperforms → Hierarchical Reasoning Model
confidence 95% · a 5M-parameter CMM achieves a state-of-the-art accuracy of 93.7 %, outperforming the 27M-parameter HRM (55.0 %)
Contraction Mapping Model → uses → Neural Ordinary Differential Equations
confidence 95% · reformulates discrete recursive reasoning into continuous Neural Ordinary and Stochastic Differential Equations (NODEs/NSDEs)
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Current large language models (LLMs) primarily rely on linear sequence generation and massive parameter counts, yet they severely struggle with complex algorithmic reasoning. While recent reasoning architectures, such as the Hierarchical Reasoning Model (HRM) and Tiny Recursive Model (TRM), demonstrate that compact recursive networks can tackle these tasks, their training dynamics often lack rigorous mathematical guarantees, leading to instability and representational collapse. We propose the Contraction Mapping Model (CMM), a novel architecture that reformulates discrete recursive reasoning into continuous Neural Ordinary and Stochastic Differential Equations (NODEs/NSDEs). By explicitly enforcing the convergence of the latent phase point to a stable equilibrium state and mitigating feature collapse with a hyperspherical repulsion loss, the CMM provides a mathematically grounded and highly stable reasoning engine. On the Sudoku-Extreme benchmark, a 5M-parameter CMM achieves a state-of-the-art accuracy of 93.7 %, outperforming the 27M-parameter HRM (55.0 %) and 5M-parameter TRM (87.4 %). Remarkably, even when aggressively compressed to an ultra-tiny footprint of just 0.26M parameters, the CMM retains robust predictive power, achieving 85.4 % on Sudoku-Extreme and 82.2 % on the Maze benchmark. These results establish a new frontier for extreme parameter efficiency, proving that mathematically rigorous latent dynamics can effectively replace brute-force scaling in artificial reasoning.
Tags
Links
- Source: https://arxiv.org/abs/2603.22871v1
- Canonical: https://arxiv.org/abs/2603.22871v1
Trouble viewing inline? Open PDF directly →
Full Text
82,660 characters extracted from source content.
Expand or collapse full text
Dynamical Systems Theory Behind a Hierarchical Reasoning Model Vasiliy A. Es’kin1,2, Mikhail E. Smorkalov2,3 1Department of Radiophysics, University of Nizhny Novgorod, Nizhny Novgorod 603950, Russia 2Huawei Nizhny Novgorod Research Center, Nizhny Novgorod 603006, Russia 3Skolkovo Institute of Science and Technology, Moscow 121205, Russia vasiliy.eskin@gmail.com, smorkalovme@gmail.com Corresponding author: Vasiliy Alekseevich Es’kin (vasiliy.eskin@gmail.com) Abstract Current large language models (LLMs) primarily rely on linear sequence generation and massive parameter counts, yet they severely struggle with complex algorithmic reasoning. While recent reasoning architectures, such as the Hierarchical Reasoning Model (HRM) and Tiny Recursive Model (TRM), demonstrate that compact recursive networks can tackle these tasks, their training dynamics often lack rigorous mathematical guarantees, leading to instability and representational collapse. We propose the Contraction Mapping Model (CMM), a novel architecture that reformulates discrete recursive reasoning into continuous Neural Ordinary and Stochastic Differential Equations (NODEs/NSDEs). By explicitly enforcing the convergence of the latent phase point to a stable equilibrium state and mitigating feature collapse with a hyperspherical repulsion loss, the CMM provides a mathematically grounded and highly stable reasoning engine. On the Sudoku-Extreme benchmark, a 5M-parameter CMM achieves a state-of-the-art accuracy of 93.7%, outperforming the 27M-parameter HRM (55.0%) and 5M-parameter TRM (87.4%). Remarkably, even when aggressively compressed to an ultra-tiny footprint of just 0.26M parameters, the CMM retains robust predictive power, achieving 85.4% on Sudoku-Extreme and 82.2% on the Maze benchmark. These results establish a new frontier for extreme parameter efficiency, proving that mathematically rigorous latent dynamics can effectively replace brute-force scaling in artificial reasoning. It was through hardship that the ape became human Folk wisdom Figure 1: Stylized visualizations of neural network models represented as volumetric Aizawa attractor trajectories. The number of model parameters is proportional to the volume of the sphere circumscribed around the corresponding attractor. Keywords Dynamical Systems Theory ⋅· Hierarchical Reasoning ⋅· Tiny Recursive Models ⋅· Neural ODEs ⋅· Stochastic Differential Equations ⋅· Algorithmic Reasoning ⋅· Parameter Efficiency ⋅· Contraction Mapping 1 Introduction The evolution of deep learning has been historically characterized by the systematic stacking of successive layers to enhance representational power, yet modern large language models (LLMs) based on the Transformer architecture remain paradoxically shallow in their fundamental computational structure [28]. Despite their success, these models operate primarily through a “System 1” modality, where the fixed-depth constraint of the Transformer block confines them to restricted computational complexity classes such as AC0AC^0 or TC0TC^0 [5, 20]. Specifically, models in TC0TC^0 are limited to solving problems that can be decomposed into parallelizable, constant-depth circuits, effectively preventing standard Transformers from executing algorithms that require a number of sequential steps proportional to the input size [9]. To circumvent these boundaries, the prevailing research paradigm has relied on Chain-of-Thought (CoT) prompting [31], which externalizes reasoning into linear, token-level sequences. However, CoT is increasingly viewed as a functional crutch rather than an architectural solution; it suffers from high inference latency and extreme data requirements [28]. Inspired by the hierarchical organization of the human brain, recent advancements have introduced Hierarchical Reasoning Models (HRM) that separate high-level strategic planning from low-level procedural execution [15]. A pivotal development in this space is the Tiny Recursive Model (TRM), which challenges traditional scaling laws by demonstrating that recursive weight-sharing can enable compact models to outperform much larger counterparts. Unlike standard Transformers, TRMs decouple computational depth from parameter count by iteratively passing latent representations through shared parameters for an arbitrary number of steps. This allows the model to expand its “thinking time” in the latent space, theoretically enabling it to transcend the TC0TC^0 barrier within a compact, recursive architecture [17]. Despite the promising performance of hierarchical and recursive architectures on discrete reasoning benchmarks like ARC-AGI, several critical bottlenecks continue to impede the development of universal reasoning systems. A primary challenge lies in the efficiency of credit assignment within these multi-level structures. Current reinforcement learning (RL) methods often fail to distinguish between critical high-level planning decisions and routine execution steps, thereby diluting the optimization signal and slowing the discovery of effective latent strategies [29]. Furthermore, there is an ongoing debate regarding the necessity of explicit hierarchical modules; recent empirical evidence suggests that deep supervision and recursive refinement might be the primary drivers of performance, raising questions about whether current hierarchical separation truly captures the underlying structure of reasoning tasks [15]. Training these models also faces systemic hurdles, including the high memory overhead of Backpropagation Through Time (BPTT) and the tendency of recurrent states to collapse into “shortcut” solutions that lack robustness [18, 17]. Perhaps most significantly, many current systems remain transductive, encoding specific task patterns into their weights rather than inducing generalizable algorithmic programs, which limits their adaptability to novel, out-of-distribution problems that require genuine “program synthesis” [2, 15]. Furthermore, the emergence of TRMs raises a fundamental question regarding the minimal parameter set required to sustain complex reasoning without sacrificing the accuracy typically associated with much larger models. Consequently, a key area of interest lies in maximizing the reasoning precision of models within the compact parameter constraints established by state-of-the-art TRMs. These issues, ranging from the optimization of recursive depth to the pursuit of extreme parameter efficiency, constitute the frontier of current research in machine reasoning. Figure 2: Architecture of hierarchical reasoning model. This work is devoted to enhancing the stability, generalization, and parameter efficiency of recursive reasoning models by analyzing them through the lens of continuous dynamical systems. We briefly list the main contributions of this paper: 1. We reformulate the discrete recursive steps of hierarchical and tiny recursive architectures as continuous Neural Ordinary Differential Equations (NODEs) and Neural Stochastic Differential Equations (NSDEs), providing a rigorous mathematical framework for latent reasoning dynamics. 2. We propose the Contraction Mapping Model (CMM), which enforces the convergence of the phase point of the system to a stable equilibrium state. This is achieved by introducing auxiliary loss terms derived directly from the equilibrium point conditions and Routh–Hurwitz stability criterion. 3. We introduce a hyperspherical repulsion loss function designed to prevent representational collapse, ensuring a uniform, orthogonal, and highly expressive distribution of task features within the latent space. 4. We develop computationally efficient polynomial approximations of the StableMax function (StableMax3 and StableMax5) and implement an adaptive loss balancing technique via Algebraic Gradient Normalization (AlgGradNorm) to stabilize the complex, multi-objective training process. 5. We empirically demonstrate that CMMs achieve highly competitive accuracy on complex algorithmic reasoning benchmarks (Sudoku-Extreme and Maze) under extreme parameter constraints, successfully compressing the reasoning engine to just 0.26M parameters while matching or exceeding the performance of significantly larger baselines. The remainder of the paper is structured as follows. Section 2 reviews the problem statement and the current state of recursive reasoning architectures, specifically detailing the Hierarchical Reasoning Model (HRM) and Tiny Recursion Model (TRM) alongside their deep supervision training paradigms. In Section 3, we present our mathematical and architectural modifications, transitioning from discrete iterative equations to continuous NODEs and NSDEs, and define the proposed CMM with its auxiliary loss functions. Section 4 showcases the results of our extensive numerical experiments and ablation studies. Finally, concluding remarks are provided in Section 5. 2 Current State and Statement of the Problem Consider a set of puzzles consisting of an input set ii=1Np \ x_i \^N_p_i=1 (questions) and an output set ii=1Np \ y_i \^N_p_i=1 (answers). Here, each input vector i x_i corresponds to the output vector i y_i, and NpN_p denotes the number of puzzles in the set. We will consider puzzles such as Sudoku-Extreme, Maze-Hard, ARC-AGI-1 and ARC-AGI-2 (see in [28]). It is necessary to find a mapping that translates the input data i x_i into the output i y_i. The process of solving the given problem can be considered as an operator :→G:X between spaces X and Y. In the case of the puzzle problem, the space of questions X and the space of answers Y are finite spaces. We assume that there exists a solution operator †G for the given problem: †:→. :X . (1) Our goal is to approximate †G by constructing a parametric map :→,∈Θ, _ θ:X , θ∈ , (2) for some finite-dimensional parameter space Θ by choosing †∈Θ θ ∈ so that †≈†G_ θ . Such an operator can be represented by a large language model, or relatively small models such as hierarchical reasoning [28] and tiny recursive [17] models. 2.1 Hierarchical Reasoning Model Consider the Hierarchical Reasoning Model (HRM) as an approximator of the problem solving operator. The architecture of the HRM is shown on Figure 2. The HRM model consists of the following learnable components: an input embedding ^input(⋅;in) E_input(·; θ_in), a low-level recurrent module ^L(⋅;L) F_ L(·; θ_ L) (L-module), a high-level recurrent module ^H(⋅;H) F_ H(·; θ_ H) (H-module), and an output network output(⋅;out) E_output(·; θ_out). Dynamics of the model unfold over NHN_ H high-level cycles of NLN_ L low-level timesteps each. The discrete timesteps of a single forward pass are indexed by i=1,…,NL×NHi=1,…,N_ L× N_ H. The modules ^L F_ L and ^H F_ H each produce hidden state tensors ^L(i) z^(i)_ L for ^L F_ L and ^H(i) z^(i)_ H for ^H F_ H, which are initialized with the tensors ^L(0) z^(0)_ L and ^H(0) z^(0)_ H, respectively. The HRM maps an input vector x to an output prediction vector y as follows. First, the input x is projected into a hidden tensor of working representation x by the input network: ^=^input(;in). x= E_input( x; θ_in). (3) At each timestep i, the L-module updates its state conditioned on its own previous state ^L(i−1) z^(i-1)_ L, current state of the H-module ^H(i−1) z^(i-1)_ H (which remains fixed throughout the cycle), and the input representation x. The H-module updates ^H(i) z^(i)_ H only once per cycle (i.e., every NLN_ L timesteps) using final state of the L-module at the end of that cycle: ^L(i)=^L(^H(i−1)+^L(i−1)+^;L), z^(i)_ L= F_ L ( z^(i-1)_ H+ z^(i-1)_ L+ x; θ_ L ), (4) ^H(i)=^H(^H(i−1)+^L(i−1);H)ifi≡0(modNL),^H(i−1)otherwise. z^(i)_ H= cases F_ H( z^(i-1)_ H+ z^(i-1)_ L; θ_ H)&if i≡ 0 (mod\,\,\,N_ L),\\[6.0pt] z^(i-1)_ H&otherwise. cases (5) Finally, after NHN_ H full cycles, a prediction y is extracted from the hidden state of the H-module: =output(^;out), y= E_output( y; θ_out), (6) where ^=^H(NLNH) y= z^(N_ LN_ H)_ H. This entire (NLNH)(N_ LN_ H)–timestep process represents a single forward pass of the HRM. 2.2 Tiny Recursion Model Tiny Recursive Models (TRMs), proposed in [17], represent a simple recursive reasoning approach using a single tiny network recursing on its latent reasoning feature and progressively improving its final answer. The L-module and H-module are represented by the same neural network (i.e., L≡H θ_ L≡ θ_ H and ^L≡^H F_ L≡ F_ H). This single module consists of two transformer layers (see Figure 2) with a self-attention layer or an MLP-mixer (MLP) [27], which is applied on the sequence length. It significantly reduces the number of parameters by halving the number of layers and replacing the two networks with a single tiny network. Another feature of the TRM that mitigates overfitting and improves the generalizing abilities of the neural network is the use of an Exponential Moving Average (EMA). As a result of applying these approaches, the TRM significantly outperformed HRM in terms of accuracy. 2.3 Training Figure 3: Diagram of training of HRM (a) and TRM (b). Listing 1: Pseudocode of HRM ⬇ def hrm(z, x, N_L=2, N_H=2): # hierarchical reasoning zH, zL = z with torch.no_grad(): for i in range(N_L*N_H - 2): zL = L_net(zL, zH, x) if (i + 1) % N_H == 0: zH = H_net(zH, zL) # 1-step grad zL = L_net(zL, zH, x) zH = H_net(zH, zL) return (zH, zL), output_head(zH), Q_head(zH) def ACT_halt(q, y_hat, y_true): target_halt = (y_hat == y_true) loss = 0.5*binary_cross_entropy(q[0], target_halt) return loss def ACT_continue(q, last_step): if last_step: target_continue = sigmoid(q[0]) else: target_continue = sigmoid(max(q[0], q[1])) loss = 0.5*binary_cross_entropy(q[1], target_continue) return loss # Deep Supervision for x_input, y_true in train_dataloader: z = z_init for step in range(N_supervision): # deep supervision x = input_embedding(x_input) z, y_pred, q = hrm(z, x) loss = softmax_cross_entropy(y_pred, y_true) # Adaptive computational time (ACT) using Q-learning loss += ACT_halt(q, y_pred, y_true) _, _, q_next = hrm(z, x) # extra forward pass loss += ACT_continue(q_next, step == N_sup - 1) z = z.detach() loss.backward() opt.step() opt.zero_grad() if q[0] > q[1]: # early-stopping break Listing 2: Pseudocode of TRM ⬇ def latent_recursion(x, y, z, N_L=6): for i in range(N_L): # latent reasoning z = net(x, y, z) y = net(y, z) # refine output answer return y, z def deep_recursion(x, y, z, N_L=6, N_H=3): # recursing N_H-1 times to improve y and z (no gradients needed) with torch.no_grad(): for j in range(N_H-1): y, z = latent_recursion(x, y, z, N_L) # recursing once to improve y and z y, z = latent_recursion(x, y, z, N_L) return (y.detach(), z.detach()), output_head(y), Q_head(y) # Deep Supervision for x_input, y_true in train_dataloader: y, z = y_init, z_init for step in range(N_supervision): x = input_embedding(x_input) (y, z), y_hat, q_hat = deep_recursion(x, y, z) loss = softmax_cross_entropy(y_hat, y_true) loss += binary_cross_entropy(q_hat, (y_hat == y_true)) loss.backward() opt.step() opt.zero_grad() if q_hat > 0: # early-stopping break Figure 4: Pseudocode of training of the HRM (left) and TRM (right). The authors of HRM proposed a one-step approximation of the HRM gradient, using the gradient of the last state of each module and treating other states as constants (see Figure 3a). The authors of the TRM modified this training step by using the last state of the H-module for gradient calculation, as well as the last states of the L-modules after the previous update of the state of the H-module (see Figure 3b). 2.3.1 Deep supervision for the HRM For a given data sample ( x, y), multiple forward passes of the HRM model are run, each of which is referred to as a segment. For each segment m∈1,…,Nsm∈ \1,…,N_s \, let ^(m)=(^H(mNLNH),^L(mNLNL)) z^(m)= ( z^(mN_ LN_ H)_ H, z^(mN_ LN_ L)_ L ) represent the hidden state at the end of segment m. At each segment m, a deep supervision step of training is applied as follows: (^(m),(m)) ( z^(m), y^(m) ) ←HRM(^∘(m−1),;), ( z^(m-1), x; θ ), ℒLM(m)() ^(m)_ LM( θ) ←Loss((m),), ←Loss ( y^(m), y ), θ ←OptimizationStep(,∇ℒLM(m)()). ←OptimizationStep ( θ, _ θL^(m)_ LM( θ) ). (7) Here, HRM is the HRM model, ^(m−1) z^(m-1) is the hidden state on the previous segment m−1m-1, (m) y^(m) is the output vector associated with the current hidden state (at segment m) ^(m) z^(m), θ is the set of all trainable parameters of HRM (=in,H,L,out θ= \ θ_in, θ_H, θ_L, θ_out \). Here and below Ra symbol over variables indicates that the computed values (^∘(m−1) z^(m-1)) are detached from the computational graph (this operation corresponds to detach in PyTorch, and lax.stop_gradient in JAX). The critical aspect of this process is that the hidden state, ^∘(m) z^(m), is “detached” from the computation graph before being used as the input state for the next segment. Consequently, gradients from segment m+1m+1 do not propagate back through segment m, effectively creating a one-step approximation of the gradient of the recursive deep supervision process [3, 24]. This approach provides more frequent feedback to the H-module and acts as a regularization method, demonstrating better empirical performance and improved stability in deep equilibrium models compared to more complicated, Jacobian-based regularization methods [3, 4]. Figure 4 shows the pseudocode for training the HRM using deep supervision. 2.3.2 Adaptive computational time (ACT) The authors of the HRM incorporated an adaptive halting strategy into the HRM. This integration leverages deep supervision and uses the Q-learning algorithm to adaptively determine the number of segments. A Q-head uses the final state of the H-module to predict the Q-values Q^(m)=(Q^halt(m),Q^continue(m)) Q^(m)= ( Q^(m)_ halt, Q^(m)_ continue ) of the “halt” and “continue” actions: Q^(m)=σ(QT^H(mNLNH)), Q^(m)=σ ( θ T_Q z^(mN_ LN_ H)_ H ), (8) where σ is a sigmoid function, QT θ T_Q are trainable parameters. Let MmaxM_ max and MminM_ min are the maximum (a fixed hyperparameter) and minimum (a random variable) numbers of segments, respectively. The value of MminM_ min is determined randomly: with a probability ε , it is chosen uniformly from the set 2,…,Mmin\2,…,M_ min\, and with a probability 1−ε1- , it is set to 11. The halt action is taken under two conditions: when the number of segments exceeds the maximum threshold MmaxM_ max, or when the estimated value of the halt Q^halt(m) Q^(m)_ halt is greater than the estimated value for continuing Q^continue(m) Q^(m)_ continue, and the number of segments has reached at least the minimum threshold MminM_ min. The Q-head is updated through a Q-learning algorithm, which is defined by Markov Decision Process (MDP). The state of the MDP at segment m is ^(m) z^(m), and the action space is halt,continue\halt,continue\. Choosing the action “halt” terminates the episode and returns a binary reward indicating prediction correctness, i.e., (m)= 1\,\ y^(m)= y\. Choosing “continue” yields a reward of 0 and the state transitions to ^(m+1) z^(m+1). Thus, the Q-learning targets for the two actions G^(m)=(G^halt(m),G^continue(m)) G^(m)= ( G^(m)_ halt, G^(m)_ continue ) are given by G^halt(m) G^(m)_ halt =(m)=, = 1\,\ y^(m)= y\, G^continue(m) G^(m)_ continue =Q^halt(m+1),ifm≥Mmaxmax(Q^halt(m+1),Q^continue(m+1)),otherwise. = cases Q^(m+1)_ halt,&if m≥ M_ max\\ ( Q^(m+1)_ halt, Q^(m+1)_ continue ),&otherwise. cases (9) The overall loss for each supervision segment combines both the Q-head loss and the sequence-to-sequence loss: ℒACT(m)():=ℒLM(m)()+ℒBCE(m)(), ^(m)_ ACT( θ):=L^(m)_ LM( θ)+L^(m)_ BCE( θ), (10) where ℒBCE(m)():=BinaryCrossEntropy(Q^(m),G^(m)). ^(m)_ BCE( θ):=BinaryCrossEntropy ( Q^(m), G^(m) ). (11) Minimizing this loss allows for both accurate predictions and near-optimal stopping decisions. The “halt” action ends the training loop. In practice, sequences are processed in batches, which can be easily handled by replacing any halted sample in a batch with a new sample from the data loader. The optimization problem can be defined as follows: ∗=argminℒACT(m)(), θ^*= \, θ \,L^(m)_ ACT( θ), (12) where ∗ θ^* are optimal parameters of the neural networks which minimize the discrepancy between the exact solution y and the approximate one (m) y^(m). 2.3.3 Training of TRM The training of the TRM is similar to the training of the HRM described above. There are two main differences. The first difference is the reduction of the need for the expensive second forward pass. The ACT in HRM through Q-learning requires two forward passes, which slows down training. The authors of the TRM proposed a simple solution, which is to eliminate the continue loss (from the Q-learning) and only learn a halting probability through a Binary-Cross-Entropy loss of having reached the correct solution. By removing the continue loss, the authors eliminated the need for the expensive second forward pass, while still being able to determine when to halt with relatively good accuracy. The second difference is the use of an Exponential Moving Average (EMA). To reduce the problem of overfitting in the HRM due to small datasets and to improve stability, the authors integrated Exponential Moving Average (EMA) of the weights, a common technique in GANs and diffusion models to improve stability. It is found that it prevents sharp collapse and leads to higher generalization. 3 Modifications of the HRM 3.1 From Discrete Equations to Neural Ordinary Differential Equations We need the continuous counterparts of equations (4) and (5) for further analysis. Next, we find the differential equations that generate these difference equations. Let the discrete variables be related to continuous time as ^L(i)=^L(ti),^H(i)=^H(ti),ti=iΔt. z^(i)_ L= z_ L(t_i), z^(i)_ H= z_ H(t_i), t_i=i t. (13) The time derivative is approximated by the forward difference: dzdt|t=ti≈z(ti)−z(ti−1)Δt=z(i)−z(i−1)Δt. . dz dt |_t=t_i≈ z(t_i)-z(t_i-1) t= z^(i)-z^(i-1) t. (14) From (4) we obtain the corresponding differential equation ∂^L(t)∂t=1Δt[^L(^H(t)+^L(t)+^;L)−^L(t)]. ∂ z_ L(t)∂ t= 1 t [ F_ L ( z_ H(t)+ z_ L(t)+ x; θ_ L )- z_ L(t) ]. (15) The equation (5) corresponds to an impulse differential equation ∂^H(t)∂t=1Δt[^H(^H(t)+^L(t);H)−^H(t)]∑j∈ℤδ(t−jNLΔt), ∂ z_ H(t)∂ t= 1 t [ F_ H ( z_ H(t)+ z_ L(t); θ_ H )- z_ H(t) ] _j δ (t-jN_ L t ), (16) where δ is the Dirac delta function. Averaging over the period NLΔtN_ L t leads to the following equation: ∂^H(t)∂t=1NLΔt[^H(^H(t)+^L(t);H)−^H(t)]. ∂ z_ H(t)∂ t= 1N_ L t [ F_ H ( z_ H(t)+ z_ L(t); θ_ H )- z_ H(t) ]. (17) Assuming that Δt=1 t=1, we get the following system of differential equations: ∂^L(t)∂t=^L(^H(t)+^L(t)+^;L)−^L(t), ∂ z_ L(t)∂ t= F_ L ( z_ H(t)+ z_ L(t)+ x; θ_ L )- z_ L(t), (18) ∂^H(t)∂t=1NL[^H(^H(t)+^L(t);H)−^H(t)]. ∂ z_ H(t)∂ t= 1N_ L [ F_ H( z_ H(t)+ z_ L(t); θ_ H)- z_ H(t) ]. (19) Note that the system of equations (4) and (5) represents a discrete Euler integration scheme with step size Δt=1 t=1 applied to the continuous-time system of equations (18)–(19). Let us rewrite this continuous system in the following more compact form ∂^(t)∂t=^(^(t),^;)−^(t), ∂ z(t)∂ t= F ( z(t), x; θ )- z(t), (20) where z denotes either ^L z_ L or ^H z_ H, and F stands for ^L F_ L or ^H F_ H, respectively. For concreteness, we consider t∈[0,T]t∈[0,T], where T is the end time point. According to the problem formulation, the final value of ^H z_ H at t=Tt=T must be equal to y y, while the initial values of ^L z_ L and ^H z_ H are not specified and can therefore be chosen arbitrarily. We assume the initial condition ^H(0)= z_ H(0)= x. Thus, the initial and final values of ^H z_ H can be expressed as follows: ^H(0)=^,^H(T)=^. z_ H(0)= x, z_ H(T)= y. (21) The system (20) belongs to the class of Neural Ordinary Differential Equations (NODEs) [6], which, subject to conditions (21), can be solved using standard ODE solvers. The objective is to optimize the neural networks within (20) such that the resulting trajectory satisfies both the differential equations and the specified initial and final conditions. It is important to highlight a key technical detail: the term x (the input tensor of the puzzle) in the system of equations (20) acts as a system parameter. This implies that for each specific puzzle, the model solves a distinct system of equations parameterized by the input tensor x. Consequently, the problem can be reformulated with a more relaxed initial condition compared to (21): the objective is to find the neural network parameters ∗ θ^* of F such that, for any initial conditions of ^L z_ L and ^H z_ H, the solution to the system of equations (20) parameterized by the input tensor x yields the puzzle solution y at time T. 3.2 Modifications of the training and model Note that a neural network can achieve the desired result (solution y) at the end of one of the first segments of adaptive computational time (Section 2.3.2). Consequently, adaptive computational time in the following segments must not cause the state to diverge from this solution. That is, such a state of the system should be stable. In the theory of dynamical systems, such a state corresponds to stable equilibrium state or an attractive set (attractor) [1, 26]. Hereafter, we will call it a stable equilibrium point, encompassing the notion of attractive sets. In this work, we show that the stability and convergence of recursive reasoning can be significantly improved by framing models like HRM and TRM as contraction mapping systems. To fully leverage this mathematical framework, we propose the Contraction Mapping Model (CMM). Unlike previous empirical methods, CMM enforces the contraction mapping property through specific architectural modifications and introduces auxiliary loss terms associated with equilibrium points. This provides a mathematically rigorous foundation for convergence of the model to a unique fixed point. Hereafter, for brevity, we will denote the CMM trained to solve the NODEs (20) subject to the conditions (21) as CMMDNODECMM^NODE_D, where D is the dimensionality of the hidden state. 3.2.1 Equilibrium points The equilibrium points of dynamic system (20) are determined from the following equation ^∗(t)=^(^∗(t),^;). z^*(t)= F ( z^*(t), x; θ ). (22) Hereafter, by equilibrium points, we mean sets that can be attractive (stable equilibrium points) or repulsive (unstable equilibrium points). Let the function F be expanded into a Taylor series in some neighborhood of the equilibrium state ^∗ z^*. Then the system (20) in this neighborhood can be rewritten as ∂zi(t)∂t=[∑k=1K∂Fi(^(t),^;)∂zk|^=^∗−1](zi(t)−zi∗). ∂z_i(t)∂ t= [Σ _k=1^K . ∂F_i ( z(t), x; θ )∂ z_k |_ z= z^*-1 ] (z_i(t)-z^*_i ). (23) Here, ziz_i and FiF_i are the iith components of the tensor z and the tensor function F (i=1,…,Ki=1,...,K). In equation (23), we neglected the nonlinear terms above the first order of smallness. The characteristic equation of the linearized system (23) is as follows D(λ)=|J11−λJ12⋯J1KJ21J22−λ⋯J2K⋮⋱⋮JK1JK2⋯JKK−λ|=0, D(λ)= vmatrixJ_11-λ&J_12&·s&J_1K\\ J_21&J_22-λ&·s&J_2K\\ & & & \\ J_K1&J_K2&·s&J_K-λ vmatrix=0, (24) where Jik=∂^i(^(t),^;)∂zk−δikJ_ik= ∂ F_i ( z(t), x; θ )∂ z_k- _ik, δik _ik is the Kronecker symbol. Expanding the determinant in (24) yields a characteristic equation whose left side is a K-th degree polynomial with respect to (λ)(λ): a0λK+a1λK−1+⋯+aK−1λ+aK=0, a_0λ^K+a_1λ^K-1+…+a_K-1λ+a_K=0, (25) where a0,a1,…,aKa_0,a_1,…,a_K are real coefficients depending on the system parameters (a0=(−1)Ka_0=(-1)^K). Criterion of Routh-Hurwitz. Let a0>0a_0>0. We form the determinants: Δ1=a1;Δ2=|a1a0a3a2|;Δ3=|a1a00a3a2a1a5a4a3|;… _1=a_1; _2= vmatrixa_1&a_0\\ a_3&a_2 vmatrix; _3= vmatrixa_1&a_0&0\\ a_3&a_2&a_1\\ a_5&a_4&a_3 vmatrix; … (26) ΔK=|a1a000⋯0a3a2a1a0⋯0a5a4a3a2⋯0⋮⋱⋮a2K−1a2K−2a2K−3a2K−4⋯aK|=aKΔK−1, _K= vmatrixa_1&a_0&0&0&·s&0\\ a_3&a_2&a_1&a_0&·s&0\\ a_5&a_4&a_3&a_2&·s&0\\ & & & & & \\ a_2K-1&a_2K-2&a_2K-3&a_2K-4&·s&a_K vmatrix=a_K _K-1, (27) where aj=0a_j=0, if j>Kj>K. For all roots of equation (25) (with real values aja_j, j=0,1,…,Kj=0,1,…,K, and given a0>0a_0>0) to have negative real parts, it is necessary and sufficient that the following inequalities hold: Δ1>0,Δ2>0,…,ΔK−1>0;aK>0. _1>0, _2>0, …, _K-1>0; a_K>0. (28) As is known, the equilibrium state of a linear time-invariant system is considered asymptotically stable (i.e., the system returns to this state after small perturbations) if all roots (λi _i) of its characteristic equation lie in the open left half-plane of the complex plane. This means that the real part (Re) of every root must be strictly negative: Re(λi)<0Re( _i)<0 for all i. Thus, if the conditions of the Routh-Hurwitz criterion are met, it guarantees that all roots of the characteristic equation lie in the left half-plane, and, consequently, the equilibrium state of the dynamic system is stable. We should note that while satisfying all K conditions of the Routh-Hurwitz criterion is strictly necessary and sufficient for asymptotic stability, computing all K determinants is computationally intractable for high-dimensional latent spaces (e.g., K∈128,512K∈\128,512\) within the inner loop of deep learning optimization. However, if the first condition Δ1>0 _1>0 is not met, the system cannot be stable. Therefore, we adopt a computationally efficient relaxation by explicitly enforcing this first necessary condition. Geometrically, enforcing a positive Δ1 _1 (which corresponds to a negative trace of the linearized system matrix) ensures that the local phase volume contracts around the equilibrium point. For the given problem, we have the following equation for Δ1 _1: Δ1=(−1)K−1a0(∑i=1KJii−K)=K−∑i=1KJii. _1= (-1)^K-1a_0 (Σ _i=1^KJ_i-K )=K-Σ _i=1^KJ_i. (29) Additional Loss Terms. Based on equations (22) and (29), we formulate the following additional loss terms for the values ^H(m) z^(m)_ H of m-th segment. Equation (22) yields the term: ℒequil(m)() ^(m)_ equil( θ) :=[^∘H(n)−^H(^H(n)+^L(n);H)]2, := [ z^(n)_ H- F_ H( z^(n)_ H+ z^(n)_ L; θ_ H) ]^2, (30) where n=mNLNHn=mN_ LN_ H. Equation (29) yields the following auxiliary loss term. Rather than enforcing the computationally intractable full criterion, this term acts as a feasible regularizer corresponding to the first necessary condition (the trace penalty) derived from the Routh-Hurwitz criterion. This encourages the formation of a stable equilibrium point at ^H∗= z^*_ H= y: ℒRHstable(m)() ^(m)_ RH\,stable( θ) :=[ReLU(1K∑i=1KJii−1)]2. := [ ReLU ( 1KΣ _i=1^KJ_i-1 ) ]^2. (31) Similarly, we introduce a corresponding loss term that encourages an unstable equilibrium point (a repeller) at ^H∗= z^*_ H= x: ℒRHunstable(m)() ^(m)_ RH\,unstable( θ) :=[ReLU(1−1K∑i=1KJii)]2. := [ ReLU (1- 1KΣ _i=1^KJ_i ) ]^2. (32) Hereafter, for brevity, we refer to the loss terms derived from this necessary condition simply as “RH stable” and “RH unstable”. Discussion on changing integral trajectories during the training process. During the training of the neural network F, and taking into account the aforementioned auxiliary loss terms, the phase portrait of the dynamical system evolves such that a stable equilibrium point ^H∗= z^*_ H= y emerges, toward which all integral trajectories converge. Figure 5 presents a sketch of the evolution of the phase portrait for the dynamical system (20)–(21). Figure 5: Evolution of the phase portrait of the dynamical system (20), (21) during neural network training. It should be noted that as the system of differential equations is integrated and the phase point approaches the stable equilibrium point ^H∗= z^*_ H= y, the time-dependent behaviour of ^L z_ L (across integration steps) should become oscillatory. This is due to the fact that upon reaching the equilibrium, if ^L z_ L deviates from its equilibrium value at step i, it must take a value at step i+1i+1 that is close to its value at step i−1i-1. If this condition is violated, the phase point leaves the neighborhood of the equilibrium ^H∗= z^*_ H= y, and the solution of NODEs fails to satisfy the second condition (21). This likely explains why an even number of steps for the L-module yields the best results for the TRM [17]. An illustration of such behaviour with an even number of steps for the L-module is presented in Fig. 6. Figure 6: Integration of the system of differential equations (20)–(21) under the trained neural network, which determines the landscape of phase space. Green and red arrows correspond to the integration stages using the L- and H-modules, respectively. 3.2.2 Repulsion Loss Term To prevent representation collapse, ensure a uniform distribution of puzzle features within the latent space, and enhance the expressive power of the neural network, a repulsion-based loss function defined on a unit hypersphere is employed. Let the hidden representation (denoted generically as ^H z_H, ^L z_L, or x) be a tensor of dimensions B×S×DB× S× D, where B is the batch size, S is the sequence length (the dimension of the input message), and D is the hidden dimension. Initially, the feature tensor of the i-th sample in the batch is flattened into a single one-dimensional vector i∈ℝMz_i ^M, where M=S×DM=S× D. These flattened vectors are subsequently projected onto the surface of a unit hypersphere via L2L_2-normalization: i=i‖i‖2,i∈1,…,B.u_i= z_i\|z_i\|_2, i∈\1,…,B\. (33) Following normalization, the Gram matrix (cosine similarity matrix) is constructed for all pairs of normalized vectors. To minimize the correlation between distinct samples, a penalty is imposed on the squared off-diagonal elements of this matrix. The final repulsion loss, ℒrepL_rep, is defined as the average of the squared inner products for all pairs i≠ji≠ j: ℒrep:=1B(B−1)∑i=1B∑j≠i⟨i,j⟩2.L_rep:= 1B(B-1) _i=1^B _j≠ i _i,u_j ^2. (34) Squaring the inner product ensures that both parallel and anti-parallel representations are equally penalized, thereby driving the latent vectors toward mutual orthogonality. The diagonal elements are explicitly excluded from the summation to avoid penalizing the inherent self-similarity of the features. An illustration of the data transformation stages and their behaviour under the influence of this loss term is presented in Fig. 7. Here, the cow’s horns, eyes, and nostrils serve as the endpoints of the feature vector (“spot cow” is taken from the repository associated with [10]). Figure 7: Illustration to repulsion loss term. (a) is initial hidden representation, (b) is projection to the surface of a unit hypersphere, (c) is projection after using repulsion loss term, (d) is hidden representation after using repulsion loss term. 3.2.3 Modifications of StableMax The authors of the HRM and TRM trained their models with the StableMax loss instead of Softmax for improved stability [23]: StableMax(xi):=s(xi)∑js(xj), StableMax(x_i):= s(x_i)Σ _js(x_j), (35) s(x):=1+x,x≥0,11−x,x<0. s(x):= \ array[]l1+x, x≥ 0,\\[5.69054pt] 11-x, x<0. array . (38) To maintain the stability inherent in StableMax and increase the accuracy inherent in Softmax, the following modifications of StableMax were used in training: StableMax3 and StableMax5. StableMax3 is defined as follows StableMax3(xi):=s3(xi)∑js3(xj), StableMax3(x_i):= s_3(x_i)Σ _js_3(x_j), (39) s3(x):=1+x(1+0.5x(1+x/3)),x≥0,11−x(1−0.5x(1−x/3)),x<0. s_3(x):= \ array[]l1+x(1+0.5x(1+x/3)), x≥ 0,\\[5.69054pt] 11-x(1-0.5x(1-x/3)), x<0. array . (42) StableMax5 is defined as follows StableMax5(xi):=s5(xi)∑js5(xj), StableMax5(x_i):= s_5(x_i)Σ _js_5(x_j), (43) s5(x):=1+x(1+0.5x(1+x(1+x(1+x/5)/4)/3)),x≥0,11−x(1−0.5x(1−x(1−x(1−x/5)/4)/3)),x<0. s_5(x):= \ array[]l1+x(1+0.5x(1+x(1+x(1+x/5)/4)/3)), x≥ 0,\\[5.69054pt] 11-x(1-0.5x(1-x(1-x(1-x/5)/4)/3)), x<0. array . (46) In Eqns. (42) and (46), 1+x(1+0.5x(1+x/3))1+x(1+0.5x(1+x/3)) and 1+x(1+0.5x(1+x(1+x(1+x/5)/4)/3))1+x(1+0.5x(1+x(1+x(1+x/5)/4)/3)) are computationally efficient Taylor series expansions of exp(x) (x) up to and including the third and fifth powers, respectively, implemented via the Horner method for polynomial evaluation. This motivates the use of the numbers 3 and 5 in the names of the proposed methods, StableMax3 and StableMax5. 3.3 Neural Stochastic Differential Equations The training process descriptions for HRM and TRM (see Figure 4 and [28, 17]) state that the training procedure includes inner loops over N_supervision steps, which implement the ACT scheme. However, the official code repositories for these models [28, 17] omit this explicit inner loop. Instead, they employ only an outer loop over batches via a train_dataloader, while the current ACT segment is tracked internally within the model. Consequently, the batch iterations and ACT segment iterations are fused. This setup implies that the initial hidden states ^i(m) z^(m)_i for a given segment m and batch i are the final hidden states ^i−1(m−1) z^(m-1)_i-1 from the previous segment and batch. This training implementation causes the phase point (representing the system states ^L z_L and ^H z_H at a given moment) to jump from one integral trajectory to another, rather than moving smoothly along a single trajectory toward a stable equilibrium point. A demonstration of the difference in behaviour between the integration described in the papers and the integration implemented in the code is shown in Fig. 8. Fig. 8(a) illustrates the behaviour of the phase point under the integration procedure when an explicit loop over NsupervisionN_supervision is present. In the case where the NsupervisionN_supervision loop and the dataset batch loop are fused, the integration trajectory of the phase point changes direction in accordance with the shifting stable equilibrium point (see Fig. 8(b)–(d)). Figure 8: (a) represents the integration of the system of differential equations (20) during training according to the descriptions in [28, 17]; (b), (c), and (d) represent the integration of the system of differential equations (20) during training according to the official code for these papers. Panels (b), (c), and (d) correspond to batches i, i+1i+1, and i+2i+2, respectively. The behaviour described above is equivalent to a situation where, for a given puzzle, the initial values of z for the m-th segment take on random values every NL×NHN_ L×N_ H integration steps. Considering this circumstance and the remark at the end of section 3.1, we have modified the NODEs (20) into the following neural stochastic differential equations (NSDEs) [19, 25]: d^(t)=[^(^(t),^;)−^(t)]dt+G(^(t),t)d^(t), d z(t)= [ F ( z(t), x; θ )- z(t) ] dt+G ( z(t),t ) d W(t), (47) where W(t)W(t) is a standard Wiener process (Brownian motion) [21], which is a continuous-time stochastic process such that ^(t+s)−^(s) W(t+s)- W(s) follows a Gaussian distribution (0,s^)N(0,s I). Here, G(^(t),t)G ( z(t),t ) represents the noise injection (which can be either additive G(t)G (t ) or multiplicative G(^(t),t)G ( z(t),t )). The discrete version of equation 47 for the additive noise can be written as ^(i)=^(^(i−1),^;)+σ^(i), z^(i)= F ( z^(i-1), x; θ )+σ ζ^(i), (48) where ^(i)∼(0,^) ζ^(i) (0, I ) is standard Gaussian noise, and σ is the magnitude of the stochastic component. The dimensions of ^(i) ζ^(i) and ^(i) z^(i) coincide. The discrete version of equation 47 for the multiplicative noise can be written as ^(i)=^(^(i−1),^;)⊙(^+σ^(i)). z^(i)= F ( z^(i-1), x; θ ) ( I+σ ζ^(i) ). (49) Hereafter, for brevity, we will denote the CMM trained to solve the NSDEs (47) as CMMDNSDECMM^NSDE_D. It is evident that the injection of noise improves both the accuracy and the generalization capabilities of the model [19, 25]. This stochasticity creates a manifold of latent states from which the phase point is forced to converge toward the stable equilibrium. As specified in equations (48) and (49), the magnitude σ acts as a regularization hyperparameter that directly determines the spatial extent of this latent state distribution. Without the explicit noise component, the model would likely be trained to reach the target state only along narrow, specific trajectories, leading to rapid overfitting. In this context, the inherent “jumps” between batches in the fused iteration cycle, combined with the explicitly added noise σ^(i)σ ζ^(i), act synergistically to prevent representational collapse. This combination promotes the formation of a more robust and smooth manifold for the latent dynamics, ultimately enhancing the performance of the solver on unseen data. 3.4 Weights of Loss Terms Taking into account the additional loss components in the optimization problem (12) (which was formulated for the HRM and TRM), for the CMM, it is necessary to replace ℒACT(m)()L^(m)_ ACT( θ) with the sum of all loss terms. Naturally, each component of the loss function must be included with a certain weight. Thus, the final loss is expressed as follows: ℒTotal(m)():=∑n∈termsλnℒn(m)(), ^(m)_ Total( θ):=Σ _ subarraycn subarray _nL^(m)_n( θ), (50) where λn _n is the weight of the n-th loss term, “terms” is the set LM,BCE,“rep x”,“rep z”,“equil x”,“equil z”,“RH stable z”,\LM,BCE,``rep x′,``rep z′,``equil x′,``equil z′,``RH stable z′, “RH unstable x”``RH unstable x′\, “x” and “z” in the names of the terms denote the loss terms for x and ^H z_ H, respectively (for example, “rep x” and “rep z” correspond to repulsion loss for x and ^H z_ H). The choice of weights λn _n can be performed empirically to balance the contribution of each loss term. The task-specific terms (LM and BCE) require larger magnitudes compared to the auxiliary stability terms (repulsion and equilibrium) to ensure that the model prioritizes the correct puzzle solution. The final values were kept constant across a portion of our experiments to maintain consistency. 3.4.1 Adaptive Loss Balancing via Algebraic Gradient Normalization (AlgGradNorm) Experience in balancing loss terms for training physics-informed neural networks (PINNs) [12, 13, 11] suggests that the most effective approach is the use of automatic balancing techniques, such as the Gradient Normalization method (GradNorm) [8]. This method dynamically adjusts the weights λn _n during training to ensure that the gradients from different loss components have similar magnitudes, preventing any single term from dominating the optimization process. To train the CMM effectively with the multi-component loss defined in (50), we employ an algebraic modification of the GradNorm algorithm [8]. For each loss component ℒn(m)L^(m)_n, the L2L_2 norm of its gradient with respect to the parameters of the last shared layer w is computed: Gn(t)=λn(t)⋅‖∇wLn(t)‖2, G_n(t)= _n(t)· \| _wL_n(t) \|_2, (51) where λn(t) _n(t) is the current weight of the n-th task. The average gradient norm across all active tasks is defined as G¯(t)=n[Gn(t)] G(t)=E_n[G_n(t)]. To synchronize the convergence speeds, the algorithm monitors the training progress of each task relative to its initial state ℒn(m)(0)L^(m)_n(0): ℒ~n(m)(t)=ℒn(m)(t)ℒn(m)(0),rn(t)=ℒ~n(m)(t)j[ℒ~j(m)(t)], L^(m)_n(t)= L^(m)_n(t)L^(m)_n(0), r_n(t)= L^(m)_n(t)E_j[ L^(m)_j(t)], (52) where rn(t)r_n(t) represents the relative inverse training rate. A value of rn(t)>1r_n(t)>1 indicates that the n-th task is converging slower than the average, necessitating an increase in its weight λn _n. The target gradient norm for each task is modulated by the asymmetry parameter α: Gntarget(t)=G¯(t)[rn(t)]α. G_n^target(t)= G(t)[r_n(t)]^α. (53) The weights are updated by calculating the ratio between the target and current norms. To ensure numerical stability, the update factor is clamped within the range [0.1,10.0][0.1,10.0]: λntemp=λn(t)clamp(Gntarget(t)Gn(t)+ϵ,0.1,10.0). _n^temp= _n(t)clamp ( G_n^target(t)G_n(t)+ε,0.1,10.0 ). (54) To prevent scale drift, the weights are renormalized such that their sum equals the number of tasks N. Finally, an exponential moving average with a smoothing factor ρ is applied to stabilize the trajectories of λn _n across iterations: λ^n(t+1)=λntempN/∑j=1Nλjtemp, λ_n(t+1)= _n^tempN /Σ _j=1^N _j^temp ., (55) λn(t+1)=ρλn(t)+(1−ρ)λ^n(t+1). _n(t+1)=ρ _n(t)+(1-ρ) λ_n(t+1). (56) To adapt to different stages of the optimization process, the reference losses ℒn(m)(0)L^(m)_n(0) in (52) are softly reset every TresetT_reset steps. This allows the system to recalibrate its expectations as the phase point enters the neighborhood of the stable equilibrium y, as described by the second condition (21). 3.5 Training The training procedure for the TRM [17] serves as the primary baseline for our experimental setup. We maintained the internal tracking of supervision steps NsuperN_ super as implemented in the original repository [17]. Furthermore, we introduced an inner accumulation loop NaccumN_ accum within each batch iteration, where Naccum≤NsuperN_ accum≤ N_ super. Notably, the case where Naccum=NsuperN_ accum=N_ super corresponds to the training methodology described in [28, 17]. Given the GPU memory constraints, which often prevent fitting the required number of samples into a single batch, we employed gradient accumulation. Under this approach, the model weights were updated only after every NgradN_ grad segments. For the same memory-efficiency reasons, gradient checkpointing [7, 16] was utilized in several experimental runs to further reduce the peak VRAM footprint. 3.6 Architecture The L-module and H-module of the TRM are represented by the same neural network (L≡H θ_ L≡ θ_ H and ^L≡^H F_ L≡ F_ H) [17]. This single module consists of two transformer layers (see Figure 2). This architecture can be further simplified by enforcing weight sharing between these two transformer layers, effectively utilizing the same layer sequentially. We mark this modification as “identical transformer layers”. To mitigate overfitting and improve generalization, the training of TRM employs an exponential moving average of the weights. Another approach to prevent overfitting and simultaneously accelerate training and inference is the aggressive reduction of the total number of model parameters. Furthermore, for the CMM, which operates as a contractive mapping system, bounded nonlinear activation functions such as tanh may be preferable to unbounded functions like SiLU to ensure stability and convergence. 4 Numerical Experiments The performance of the proposed methods was evaluated on the Sudoku-Extreme dataset [28] as the main dataset. The neural networks and their training were implemented in the PyTorch framework [22] (version 2.6 under CUDA 12.4). Recognizing that many independent researchers and startups have limited access to high-performance computing clusters, efforts were focused on achieving results for the CMM that are superior or comparable to those of the TRM while operating within a restricted memory budget. Specifically, the training can be performed on hardware with 16 GB of VRAM (such as V100 or T4 GPUs) without requiring support for advanced precision formats, like bfloat16, and latest versions of flash attention. To accelerate the computations, automatic mixed precision (AMP) was employed. All models were trained using the Adam-Atan2 optimizer [14, 30]. The designations in the text are as follows: B is the batch size, D is the hidden state size. The weights of the loss terms are as follows: λLM=1.0 _LM=1.0, λBCE=0.5 _BCE=0.5, λrep x=103 _rep x=10^3, λrep z=103 _rep z=10^3, λequil x=1 _equil x=1, λequil z=1 _equil z=1, λRH stable z=104 _RH stable z=10^4, λRH unstable x=10 _RH unstable x=10. Both training and inference were executed using compiled code via torch.compile (as in [28, 17]). It is noteworthy that models trained without such compilation demonstrate a significant degradation in accuracy. This phenomenon is attributed to the operator fusion and the utilization of Fused Multiply-Add (FMA) instructions during the compilation process, which minimize intermediate rounding errors. By performing multiple operations within GPU registers before writing back to global memory, the compiled mode preserves higher numerical precision, which is particularly critical for the stability of the iterative dynamics in our model. 4.1 Baseline tests The baselines for the present study are established by comparing the proposed models against existing architectures. Tables 1 and 2 summarize the performance of DeepSeek R1, the original HRM, and the original TRM architectures on the Sudoku-Extreme and Maze benchmarks, respectively. The baseline data for these reference models are taken directly from the literature [28, 17]. To validate the experimental setup under the constrained 16 GB VRAM budget, the HRM and TRM baselines were reproduced using Automatic Mixed Precision (AMP) and gradient accumulation. While the reproduced TRM results for the Maze benchmark closely match the originally reported metrics, a noticeable performance drop is observed in the Sudoku-Extreme reproduction for both HRM and TRM. This discrepancy indicates the sensitivity of these baseline architectures to hardware-imposed training constraints, such as reduced numerical precision and limited batch configurations. Nevertheless, establishing these hardware-constrained baselines provides a fair and realistic ground for evaluating the subsequent architectural modifications. Method Accuracy % Number of Features Parameters DeepSeek R1 [28] 0.0 671B – HRM [28] 55.0 27M – TRM [17] 87.4 5M MLP-mixer TRM [17] 74.7 7M Attention HRM (our run) 45.9 27M AMP TRM (our run) 79.1 5M AMP, MLP-mixer Table 1: Test accuracy on Sudoku-Extreme benchmark. Hidden size of HRM and TRM is 512. Method Accuracy % Number of Features Parameters DeepSeek R1 [28] 0.0 671B – HRM [28] 74.5 27M – TRM [17] 85.3 7M Attention TRM (our run) 84.2 7M AMP (gradient accumulation on 10 segments) Table 2: Test accuracy on Maze benchmark. Hidden size of HRM and TRM is 512. 4.2 Experiments with CMM512NODECMM^NODE_512 Tables 3-6 detail the ablation studies and optimization of the CMM512NODECMM^NODE_512 architecture. An early finding (Table 3) indicates that replacing the unbounded SiLU activation function with the bounded tanh function can improve the stability of the latent trajectories, which naturally aligns with the requirements of a contraction mapping system. Further tuning of the high-level (NHN_ H) and low-level (NLN_ L) steps (Table 4) demonstrates that deeper recursive configurations yield the highest performance. Integrating the modified StableMax3 function alongside the Routh-Hurwitz stability loss terms (stable for ^H z_ H and unstable for x) systematically pushes the metric upwards. Table 6 illustrates the compound effect of these modifications. Note that calculating the Routh-Hurwitz stability loss terms is a computationally expensive operation and, as our experiments show, gives an accuracy close to that of a model trained using equilibrium losses. The introduction of the repulsion loss on the hypersphere, combined with the equilibrium loss, significantly improves the latent space representation. As a result, the CMM512NODECMM^NODE_512 reaches a substantially higher peak performance when utilizing an exponential decay learning rate scheduler. Furthermore, it is observed that StableMax5 provides performance on par with StableMax3 and slightly better than the original StableMax. Features Accuracy % B=256B=256 69.6 SiLU to tanh , B=256B=256 81.8 SiLU to tanh , B=64B=64, Ngrad=4N_ grad=4, using the all states of the H-module for gradient calculation 56.8 Table 3: Test accuracy of TRM on Sudoku-Extreme benchmark. General features: D=512D=512, AMP, Nsuper=16N_ super=16, NH=3N_ H=3, NL=6N_ L=6 MLP-mixer, number of parameters is 5M, Naccum=NsuperN_ accum=N_ super. NHN_ H NLN_ L Accuracy % Features 1 2 58.9 tanh is activation function 1 2 63.8 1 4 68.9 1 6 68.7 2 1 64.6 2 2 70.5 2 4 75.3 2 6 76.4 2 4 78.3 StableMax3 2 4 77.1 StableMax3, tanh is activation function 1 2 64.2 StableMax3, RH unstable x and RH stable ^H z_ H 2 4 77.9 StableMax3, RH unstable x and RH stable ^H z_ H 2 4 80.3 StableMax3, Naccum=2N_ accum=2 3 6 84.1 StableMax3, Naccum=2N_ accum=2 2 4 80.5 StableMax3, Naccum=2N_ accum=2, equilibrium in x and ^H z_ H 3 6 84.9 StableMax3, Naccum=2N_ accum=2, equilibrium in x and ^H z_ H Table 4: Dependence of accuracy of CMM512NODECMM^NODE_512 on NHN_ H and NLN_ L for Sudoku-Extreme benchmark. General features: D=512D=512, AMP, Nsuper=16N_ super=16, Ngrad=1N_ grad=1, B=256B=256, MLP-mixer, number of parameters is 5M, SiLU is activation function for main part of experiments, StableMax is loss for main part of experiments. Features Accuracy % Naccum=3N_ accum=3, Ngrad=1N_ grad=1 79.2 Naccum=2N_ accum=2, Ngrad=1N_ grad=1, using the all states of the H-module for gradient calculation 74.6 Naccum=2N_ accum=2, Ngrad=1N_ grad=1 78.8 Naccum=2N_ accum=2, Ngrad=3N_ grad=3 81.8 Naccum=2N_ accum=2, Ngrad=4N_ grad=4 81.3 Naccum=2N_ accum=2, Ngrad=3N_ grad=3, repulsion for ^H z_ H, 84.1 Table 5: Test accuracy of CMM512NODECMM^NODE_512 on Sudoku-Extreme benchmark. General features: D=512D=512, B=256B=256, AMP, Nsuper=16N_ super=16, NrmH=2N_rmH=2, NrmH=4N_rmH=4, StableMax3, equilibrium in x and ^H z_ H, MLP-mixer, number of parameters is 5M. Features Accuracy % B=256B=256, Ngrad=3N_ grad=3, without repulsion for x 88.0 B=128B=128, Ngrad=6N_ grad=6 89.9 B=256B=256, Ngrad=3N_ grad=3, tanh is activation function 79.4 B=256B=256, Ngrad=3N_ grad=3 90.1 B=384B=384, Ngrad=3N_ grad=3 89.9 B=250B=250, Ngrad=4N_ grad=4 90.5 B=500B=500, Ngrad=2N_ grad=2 90.2 B=250B=250, Ngrad=4N_ grad=4, StableMax5 90.5 B=250B=250, Ngrad=4N_ grad=4, StableMax 89.9 B=250B=250, Ngrad=4N_ grad=4, identical transformer layers (number of parameters is 2.5M) 89.9 B=250B=250, Ngrad=4N_ grad=4, D=768D=768, identical transformer layers 88.1 B=250B=250, Ngrad=4N_ grad=4, learning rate is 3×10−43× 10^-4, without equilibrium in x 75.2 B=250B=250, Ngrad=4N_ grad=4, learning rate is 6×10−56× 10^-5, without equilibrium in x 91.5 B=250B=250, Ngrad=4N_ grad=4, without equilibrium in x scheduler: exponential decay, from 10−410^-4 to 5⋅10−55· 10^-5 at 30000 epochs, 3⋅10−53· 10^-5 from 30000 to 50000 epochs 91.6 Table 6: Test accuracy of CMM512NODECMM^NODE_512 on Sudoku-Extreme benchmark. General features: D=512D=512, AMP, Nsuper=16N_ super=16, Naccum=2N_ accum=2, NH=3N_ H=3, NL=6N_ L=6, StableMax3, equilibrium in x and ^H z_ H, repulsion for ^H z_ H and x, MLP-mixer, number of parameters is 5M, learning rate is 10−410^-4. In Table 6, it is noted that the accuracy using StableMax5 is similar to the accuracy using StableMax3. The maximum accuracy of the model when learning with StableMax5 is achieved later than when training with StableMax3, and earlier than when training with StableMax. 4.3 Experiments with CMM512NSDECMM^NSDE_512 Regarding the stochastic formulation, Tables 7 and 8 present the results for CMM512NSDECMM^NSDE_512, which utilizes additive Gaussian noise (σ=0.01σ=0.01). In deterministic models, rapid overfitting is a common issue due to the phase point collapsing into narrow trajectories. The injection of noise creates a robust manifold of initial states, forcing the system to learn a broader basin of attraction towards the stable equilibrium. By scaling up the gradient accumulation steps (Ngrad=64N_ grad=64) and increasing the accumulation sequence in one batch (Naccum=16N_ accum=16), a highly stable optimization process is achieved. This configuration yields a significant performance improvement, setting a new peak score on the Sudoku-Extreme dataset for models of this scale. Furthermore, applying the AlgGradNorm technique and fine-tuning the epsilon parameter of Adam optimizer (to 10−1410^-14) with weight decay (Table 8) provides a highly stable training dynamic. This allows the model to consistently reach superior results without suffering from gradient domination by any single auxiliary loss term. Features Accuracy % Nsuper=8N_ super=8, learning rate 3×10−43× 10^-4 87.4 Nsuper=16N_ super=16, learning rate 2×10−42× 10^-4 91.7 Nsuper=16N_ super=16, learning rate 10−410^-4 91.1 Nsuper=16N_ super=16, Naccum=16N_ accum=16, Ngrad=64N_ grad=64 learning rate 2×10−42× 10^-4 93.7 B=1000B=1000, Nsuper=16N_ super=16, Naccum=16N_ accum=16, Ngrad=16N_ grad=16 learning rate 2×10−42× 10^-4, gradient checkpointing 92.7 Table 7: Test accuracy of CMM512NSDECMM^NSDE_512 on Sudoku-Extreme benchmark. General features: D=512D=512, B=250B=250, Ngrad=4N_ grad=4, AMP, Naccum=2N_ accum=2, NH=3N_ H=3, NL=6N_ L=6, StableMax3, equilibrium in ^H z_ H, repulsion for ^H z_ H and x, MLP-mixer, σ=0.01σ=0.01. Training the model using AlgGradNorm is significantly faster than the baseline training, which causes the D=512D=512 network to overfit quickly. In contrast, a network with D=128D=128 avoids overfitting despite the accelerated training pace, owing to its substantially reduced parameter count. Features Accuracy % Nsuper=16N_ super=16, Ngrad=64N_ grad=64, λLM=2 _ LM=2 80.2 Nsuper=2N_ super=2, Ngrad=8N_ grad=8, the optimizer Adam: epsilon is 10−1410^-14, weight decay is 2 89.8 Nsuper=3N_ super=3, Ngrad=12N_ grad=12, the optimizer Adam: epsilon is 10−1410^-14, weight decay is 2 91.7 Table 8: Test accuracy of CMM512NSDECMM^NSDE_512 on Sudoku-Extreme benchmark with gradient normalization (AlgGradNorm). General features: D=512D=512, B=250B=250, AMP, Naccum=16N_ accum=16, NH=3N_ H=3, NL=6N_ L=6, StableMax3, equilibrium in ^H z_ H, repulsion for ^H z_ H and x, MLP-mixer, σ=0.01σ=0.01, weight decay is 1, learning rate 2×10−42× 10^-4, . 4.4 Experiments with CMM128NODECMM^NODE_128 and CMM128NSDECMM^NSDE_128 Tables 9 and 10 explore the extreme limits of parameter efficiency. By aggressively reducing the hidden dimension to D=128D=128, the baseline TRM and standard CMMNODECMM^NODE models struggle severely, experiencing a drastic drop in performance due to the lack of expressivity in the tiny 0.52M parameter space. However, the application of the NSDE framework enables a dramatic recovery. By increasing the batch size to B=1000B=1000 to stabilize the noise estimation, employing AlgGradNorm to balance the complex multi-component loss, and freezing the input embeddings after 2500 epochs, the predictive capability is largely restored. Remarkably, when weight-sharing is enforced across the recursive steps (identical transformer layers), the parameter count is halved to an ultra-tiny 0.26M parameters. Adjusting the Adam epsilon parameter to an extremely small value (10−1410^-14) to accommodate the micro-adjustments in the gradients allows this miniature CMM128NSDECMM^NSDE_128 model to reach a highly competitive performance level (Table 10). This demonstrates that enforcing a mathematically rigorous contraction mapping allows a model with merely a quarter of a million parameters to rival the performance of baseline architectures that are substantially larger. NHN_ H NLN_ L Features Accuracy % Number of Parameters 2 4 CMM128NODECMM^NODE_128, B=256B=256 20.3 0.52M 3 6 TRM, B=250B=250, Ngrad=4N_ grad=4, Nsuper=8N_ super=8, Attention 23.8 0.52M 3 6 CMM128NODECMM^NODE_128, B=250B=250, Ngrad=4N_ grad=4, learning rate is 3×10−43×10^-4, Attention 31.4 0.52M 2 4 CMM128NSDECMM^NSDE_128, B=250B=250 Ngrad=4N_ grad=4, Nsuper=8N_ super=8 63.4 0.52M 3 6 CMM128NSDECMM^NSDE_128, B=250B=250 Ngrad=4N_ grad=4, Nsuper=8N_ super=8, learning rate is 3×10−43×10^-4 71.4 0.52M 3 6 CMM128NSDECMM^NSDE_128, B=1000B=1000 Ngrad=1N_ grad=1, Nsuper=8N_ super=8, learning rate is 3×10−43×10^-4 74.4 0.52M 3 6 CMM128NSDECMM^NSDE_128, B=250B=250, Naccum=8N_ accum=8, Ngrad=32N_ grad=32, Nsuper=8N_ super=8, learning rate is 10−410^-4 34.8 0.52M 3 6 CMM128NSDECMM^NSDE_128, B=250B=250, Naccum=16N_ accum=16, Ngrad=64N_ grad=64, Nsuper=16N_ super=16, learning rate is 2×10−32× 10^-3 70.9 0.52M 3 6 CMM128NSDECMM^NSDE_128, B=250B=250, Naccum=4N_ accum=4, Ngrad=16N_ grad=16, Nsuper=4N_ super=4, learning rate is 2×10−32× 10^-3 61.0 0.52M 3 6 CMM128NSDECMM^NSDE_128, B=250B=250, Naccum=4N_ accum=4, Ngrad=16N_ grad=16, Nsuper=4N_ super=4, learning rate is 2×10−32× 10^-3 freeze embedding after 2500 epoch 61.5 0.52M 3 6 CMM128NSDECMM^NSDE_128, B=1000B=1000, Naccum=16N_ accum=16, Ngrad=16N_ grad=16, Nsuper=16N_ super=16, learning rate is 2×10−32× 10^-3 freeze embedding after 2500 epoch, AlgGradNorm, 85.0 0.52M 3 6 CMM128NSDECMM^NSDE_128, B=1000B=1000, Naccum=16N_ accum=16, Ngrad=16N_ grad=16, Nsuper=16N_ super=16, learning rate is 2×10−32× 10^-3 freeze embedding after 2500 epoch, AlgGradNorm, identical transformer layers 84.7 0.26M Table 9: Dependence of accuracy of CMM on NHN_ H and NLN_ L for Sudoku-Extreme benchmark. General features: D=128D=128, AMP, Nsuper=16N_ super=16, Naccum=2N_ accum=2, Ngrad=3N_ grad=3, MLP-mixer, number of parameters is 0.5M, StableMax3 is the loss function for the majority of the experiments, equilibrium in x and ^H z_ H, repulsion for ^H z_ H and x, learning rate is 10−410^-4, epsilon of the optimizer Adam is 10−1210^-12, ^H z_ H and ^L z_ L are initialized by x and 0, respectively, at the beginning of the 1st segment of training. Features Accuracy % Number of Parameters initialization of ^H z_ H and ^L z_ L as in baseline, epsAdam =10−810^-8 77.8 0.52M epsAdam = 10−810^-8 81.8 0.52M epsAdam = 10−1210^-12 85.0 0.52M epsAdam = 10−1410^-14 86.1 0.52M epsAdam = 10−1410^-14, identical transformer layers 85.4 0.26M Table 10: Dependence of accuracy of CMM on different features for Sudoku-Extreme benchmark. General features: D=128D=128, NH=3N_ H=3, NL=6N_ L=6, AMP, Nsuper=16N_ super=16, Naccum=16N_ accum=16, Ngrad=16N_ grad=16, B=1000B=1000, MLP-mixer, number of parameters is 0.5M, StableMax3 is the loss function for the majority of the experiments, equilibrium in x and ^H z_ H, repulsion for ^H z_ H and x, learning rate is 2×10−32× 10^-3, freeze embedding after 2500 epoch, AlgGradNorm, epsilon of the optimizer Adam is epsAdam = 10−1210^-12, ^H z_ H and ^L z_ L are initialized by x and 0, respectively, at the beginning of the 1st segment of training. To verify that the proposed architectural enhancements and loss formulations are not restricted to the Sudoku dataset, the ultra-tiny CMM architecture is also evaluated against the baselines established in Table 2 on the Maze benchmark (Table 11). Under the highly constrained setting (D=128D=128, identical transformer layers, 0.26M parameters, and Adam epsilon =10−14=10^-14), the model successfully converges and maintains strong predictive performance. This result confirms that the NSDE framework, and the equilibrium and repulsion loss penalties provide robust, task-agnostic regularization that significantly enhances the algorithmic reasoning capabilities of extreme small-scale neural networks. Features Accuracy % Number of Parameters λLM=10 _ LM=10, epsAdam = 10−1410^-14, identical transformer layers 82.2 0.26M Table 11: Dependence of accuracy of CMM on different features for Maze benchmark. General features: D=128D=128, NH=3N_ H=3, NL=4N_ L=4, AMP, Nsuper=16N_ super=16, Naccum=4N_ accum=4, Ngrad=4N_ grad=4, B=100B=100, MLP-mixer, StableMax3 is the loss function for the majority of the experiments, equilibrium in x and ^H z_ H, repulsion for ^H z_ H and x, learning rate is 2×10−42× 10^-4, freeze embedding after 5000 epoch, AlgGradNorm, epsilon of the optimizer Adam is epsAdam = 10−1210^-12, ^H z_ H and ^L z_ L are initialized by x and 0, respectively, at the beginning of the 1st segment of training. 5 Conclusions In this paper, we introduced the Contraction Mapping Model (CMM), a novel architecture that bridges the gap between compact recursive reasoning models and the rigorous mathematical framework of continuous dynamical systems. By reformulating the discrete iterative steps of baseline hierarchical models into Neural Ordinary and Stochastic Differential Equations, we provided a principled approach to stabilizing latent reasoning dynamics. We demonstrated that enforcing contraction mapping properties, by using equilibrium constraints and the Routh-Hurwitz stability criterion, effectively guides the phase point of the system toward a stable, correct solution. Furthermore, the integration of a hyperspherical repulsion loss successfully mitigates representational collapse, while the use of adaptive loss balancing via AlgGradNorm and polynomial StableMax approximations ensures robust multi-objective optimization. Our empirical results on complex algorithmic benchmarks, specifically Sudoku-Extreme and Maze, underscore the profound advantages of this approach. The injection of noise within the NSDE framework not only prevented rapid overfitting but also allowed the model to construct a resilient, broad basin of attraction in the latent space. In particular, we achieved highly accurate predictions while aggressively compressing our model to just 0.26 million parameters, a significant reduction compared to existing approaches, without sacrificing performance. These findings suggest that massive parameter scaling is not a strict prerequisite for complex algorithmic reasoning if the underlying latent dynamics are mathematically constrained and properly optimized. Future work will explore the application of CMMs to broader, open-ended reasoning benchmarks such as ARC-AGI, as well as the potential integration of these ultra-efficient, mathematically rigorous reasoning engines into resource-constrained edge devices and larger foundational models. References [1] A. A. Andronov, E. A. Leontovich, I. I. Gordon, and A. G. Maier (1973) Qualitative theory of second-order dynamic systems. J. Wiley. Cited by: §3.2. [2] ARC-Foundation (2025) The hidden drivers of hrm’s performance on arc-agi. External Links: Link Cited by: §1. [3] S. Bai, Z. Geng, Y. Savani, and J. Z. Kolter (2022) Deep equilibrium optical flow estimation. External Links: 2204.08442, Link Cited by: §2.3.1. [4] S. Bai, V. Koltun, and J. Z. Kolter (2021) Stabilizing equilibrium models by jacobian regularization. External Links: 2106.14342, Link Cited by: §2.3.1. [5] T. Bylander (1991) Complexity results for planning. In Proceedings of the 12th International Joint Conference on Artificial Intelligence - Volume 1, IJCAI’91, San Francisco, CA, USA, p. 274–279. External Links: ISBN 1558601600 Cited by: §1. [6] R. T. Q. Chen, Y. Rubanova, J. Bettencourt, and D. Duvenaud (2019) Neural ordinary differential equations. External Links: 1806.07366, Link Cited by: §3.1. [7] T. Chen, B. Xu, C. Zhang, and C. Guestrin (2016) Training deep nets with sublinear memory cost. External Links: 1604.06174, Link Cited by: §3.5. [8] Z. Chen, V. Badrinarayanan, C. Lee, and A. Rabinovich (2018) GradNorm: gradient normalization for adaptive loss balancing in deep multitask networks. External Links: 1711.02257, Link Cited by: §3.4.1, §3.4.1. [9] D. Chiang (2025) Transformers in uniform tc0. External Links: 2409.13629, Link Cited by: §1. [10] K. Crane, U. Pinkall, and P. Schröder (2013) Robust fairing via conformal curvature flow. ACM Transactions on Graphics (TOG) 32 (4), p. 1–10. Cited by: §3.2.2. [11] V. A. Es’ kin, D. Davydov, J. Gur’eva, A. O. Malkhanov, and M. E. Smorkalov (2025) Separable physics-informed neural networks for solving elasticity problems. Computational Mathematics and Mathematical Physics 65 (9), p. 2260–2275. Cited by: §3.4.1. [12] V. A. Es’kin, D. V. Davydov, E. D. Egorova, A. O. Malkhanov, M. A. Akhukov, and M. E. Smorkalov (2024-12) About Modifications of the Loss Function for the Causal Training of Physics-Informed Neural Networks. Doklady Mathematics 110 (S1), p. S172–S192. External Links: Document, ISSN 1064-5624, Link Cited by: §3.4.1. [13] V. A. Es’kin, A. O. Malkhanov, and M. E. Smorkalov (2025) Are two hidden layers still enough for the physics-informed neural networks?. Journal of Computational Physics 537, p. 114085. External Links: ISSN 0021-9991, Document, Link Cited by: §3.4.1. [14] K. Everett, L. Xiao, M. Wortsman, A. A. Alemi, R. Novak, P. J. Liu, I. Gur, J. Sohl-Dickstein, L. P. Kaelbling, J. Lee, and J. Pennington (2024) Scaling exponents across parameterizations and optimizers. External Links: 2407.05872, Link Cited by: §4. [15] R. Ge, Q. Liao, and T. Poggio (2025) Hierarchical reasoning models: perspectives and misconceptions. External Links: 2510.00355, Link Cited by: §1, §1. [16] A. Gruslys, R. Munos, I. Danihelka, M. Lanctot, and A. Graves (2016) Memory-efficient backpropagation through time. External Links: 1606.03401, Link Cited by: §3.5. [17] A. Jolicoeur-Martineau (2025) Less is more: recursive reasoning with tiny networks. External Links: 2510.04871, Link Cited by: §1, §1, §2.2, §2, Figure 8, §3.2.1, §3.3, §3.5, §3.6, §4.1, Table 1, Table 1, Table 2, §4. [18] T. P. Lillicrap and A. Santoro (2019) Backpropagation through time and the brain. Current Opinion in Neurobiology 55, p. 82–89. Note: Machine Learning, Big Data, and Neuroscience External Links: ISSN 0959-4388, Document, Link Cited by: §1. [19] X. Liu, T. Xiao, S. Si, Q. Cao, S. Kumar, and C. Hsieh (2019) Neural sde: stabilizing neural ode networks with stochastic noise. External Links: 1906.02355, Link Cited by: §3.3, §3.3. [20] W. Merrill and A. Sabharwal (2025) A logic for expressing log-precision transformers. External Links: 2210.02671, Link Cited by: §1. [21] B. Øksendal (2003) Stochastic differential equations. In Stochastic Differential Equations: An Introduction with Applications, p. 65–84. External Links: ISBN 978-3-642-14394-6, Document, Link Cited by: §3.3. [22] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, A. Desmaison, A. Köpf, E. Yang, Z. DeVito, M. Raison, A. Tejani, S. Chilamkurthy, B. Steiner, L. Fang, J. Bai, and S. Chintala (2019) PyTorch: An Imperative Style, High-Performance Deep Learning Library. In Proceedings of the 33rd International Conference on Neural Information Processing Systems, Cited by: §4. [23] L. Prieto, M. Barsbey, P. A. M. Mediano, and T. Birdal (2025) Grokking at the edge of numerical stability. External Links: 2501.04697, Link Cited by: §3.2.3. [24] Z. Ramzi, F. Mannel, S. Bai, J. Starck, P. Ciuciu, and T. Moreau (2023) SHINE: sharing the inverse estimate from the forward pass for bi-level optimization and implicit models. External Links: 2106.00553, Link Cited by: §2.3.1. [25] M. Shen and C. Cheng (2025) Neural sdes as a unified approach to continuous-domain sequence modeling. External Links: 2501.18871, Link Cited by: §3.3, §3.3. [26] L. P. Shilnikov, A. L. Shilnikov, D. V. Turaev, and L. O. Chua (2001) Methods of qualitative theory in nonlinear dynamics (part i). Vol. 5, World Scientific. Cited by: §3.2. [27] I. Tolstikhin, N. Houlsby, A. Kolesnikov, L. Beyer, X. Zhai, T. Unterthiner, J. Yung, A. Steiner, D. Keysers, J. Uszkoreit, M. Lucic, and A. Dosovitskiy (2021) MLP-mixer: an all-mlp architecture for vision. External Links: 2105.01601, Link Cited by: §2.2. [28] G. Wang, J. Li, Y. Sun, X. Chen, C. Liu, Y. Wu, M. Lu, S. Song, and Y. A. Yadkori (2025) Hierarchical reasoning model. External Links: 2506.21734, Link Cited by: §1, §2, §2, Figure 8, §3.3, §3.5, §4.1, Table 1, Table 1, Table 2, Table 2, §4, §4. [29] H. Wang, Q. Xu, C. Liu, J. Wu, F. Lin, and W. Chen (2025) Emergent hierarchical reasoning in llms through reinforcement learning. External Links: 2509.03646, Link Cited by: §1. [30] P. Wang (2024) Adam-atan2: a pytorch implementation. External Links: Link Cited by: §4. [31] J. Wei, X. Wang, D. Schuurmans, M. Bosma, B. Ichter, F. Xia, E. Chi, Q. Le, and D. Zhou (2023) Chain-of-thought prompting elicits reasoning in large language models. External Links: 2201.11903, Link Cited by: §1.