Scaling Up Thermodynamic AI Models

# Scaling Up Thermodynamic AI Models
Source: [https://arxiv.org/html/2607.00170](https://arxiv.org/html/2607.00170)
\(June 22, 2026\)

###### Abstract

Thermodynamic computing devices based on the Ising model show great promise for low\-power AI inference and edge computing, but scalable methods for training large models for such hardware remain limited\. Prior theory shows that the time\-averaged behavior of high\-temperature Gibbs\-sampled Ising systems can implement feed\-forward neural inference\. We turn this theoretical correspondence into a scalable and purely backpropagation\-based algorithm for training deep convolutional networks for thermodynamic inference on Ising machine hardware\. Our image classification models achieve accuracies of 94\.9% on CIFAR\-10 and 76\.0% on CIFAR\-100 under binary Gibbs sampling\. We then develop and experimentally validate a mathematical theory relating inference cost to accuracy and controlling autocorrelation times\. Subsequently, we calculate asymptotic results showing that inference cost is bounded by a well\-controlled tradeoff with performance and exhibit algorithms for computing optimal inference schedules\. Finally, we discuss implications for hardware development and the future of high\-temperature thermodynamic AI models\.

## 1Introduction

AI models aggregate many small computations into statistical inferences\. However, running an AI model on a GPU, even a quantized one, entails exact error\-free computation of every activation\. This exactitude at a microscopic level, however, is not necessary for macroscopic performance: small errors in activations generally do not cause catastrophic performance loss\. This is the basic principle behind weight quantization methods and dropout layers\. On the other hand, keeping those computations error\-free comes with an energy cost at the hardware level\. In other words, keeping each computation exact during inference is wasteful\. Thermodynamic computing hardware, including both stochastic silicon and other probabilistic substrates, promises to performapproximatelocal computations at a tiny fraction of the energy expenditure\. At the macro level, therefore, running a model on such hardware could perform computations far more efficiently without significant accuracy loss\. Given the massive and ever\-growing energy demands of AI inference, the need to find a more efficient solution is likely to become more pressing over time\.

Thermodynamic computing systems exploit thermal fluctuations, stochastic dynamics, Boltzmann distributions, or non\-equilibrium processes rather than suppressing noise as an error source\[CCC2019ThermodynamicComputingWorkshop\]\. A canonical example is the Ising computer: any type of hardware which implements a system of binary spins whose probability distribution is governed by the Ising model\. Ising computation has a long history with AI, going back to the 80s\. The Hopfield network, a canonical example of the connection between statistical mechanics and machine learning, was explicitly based on the low\-temperature Ising model\. The Boltzmann machine, a positive\-temperature cousin, is still relevant today\. However, the crown jewels of modern AI inference, such as image classification and natural language processing, have largely remained outside their reach\.

Ising\-based systems have always struggled with a lack of truly scalable and efficient training algorithms for general workloads\. This has mostly confined Ising computation in practice to those situations in which a good training method does exist\. In applications where a Boltzmann machine is suitable, contrastive divergence or equilibrium propagation provides a solution, but these methods do not scale efficiently to large benchmark problems\[Nguyen\_2017\]\. Ising computers have also found wide application in combinatorial optimization, where weights can easily be derived analytically from the structure of the problem\[Lucas\_2014\]\. However, the lack of a purely backpropagation\-based approach has hampered the ability of the technology to tackle large\-scale AI workloads: in short, the best methods we know for training models have not yet been applicable\. Thermodynamic AI should not require abandoning the modern deep\-learning software stack\. This paper aims to rectify that\. We take an unorthodox approach to Ising computation, develop new backpropagation\-based algorithms for training thermodynamic models, demonstrate our results on benchmark tasks, and provide theory to explain our results\.

### 1\.1Main Contributions

The direct precursor to this work is\[moore2025correspondenceisingmachinesneural\], which establishes a correspondence theorem connecting shallow tanh neural networks to Ising systems by a global scaling parameterδ\\delta\. In theory, it promises that neural networks can be run directly on Ising machines\. More specifically, by uniformly scaling down the interactions of the second layer interactions, the average behavior of the Ising system approximates the activations of a neural network with the same topology, weights, and tanh activation\. However, that result only shows that a trained shallow perceptron runs on an Ising machine, and only in a low\-coupling and infinite\-time limit\. Furthermore, it does not address any concerns about mixing times or dynamics\. In other words, the result remains a mathematical possibility, not a practical recipe\. To make thermodynamic AI models a reality at practical workload scales, we need to be able to train deep models, control inference times, and ensure that the correspondence holds in a meaningful way at practical coupling scales\.

This paper is an algorithmic continuation and application of that work\. We develop algorithms to train deep neural networks as a series ofthermodynamic blocks, where each block meets the criteria of the correspondence theorem \([Section2](https://arxiv.org/html/2607.00170#S2)\)\. Thus, each block in the network can be implemented as a stochastic Ising machine, average\-sampled at high temperature, with binary data passed between blocks\. Furthermore, we introduce a method for ensuring that our models function at practical coupling scales, not only in a low\-coupling limit \([Section3](https://arxiv.org/html/2607.00170#S3)\)\. Our training algorithm is based on pure backpropagation through regularization terms, ensuring that it scales with problem size\. Finally, we develop a convergence theory and error rate estimates \([Section5](https://arxiv.org/html/2607.00170#S5)\), demonstrating that inference cost and mixing times remain under control, a necessary ingredient in using a stochastic system for practical applications\. We develop asymptotic estimates and closed\-form approximate formulae for the error rate as a function of inference time, then show, within a linearized red\-black model, that the autocorrelation time goes to one in the low\-coupling limit\. This theory allows us to easily derive optimal inference\-time schedules and be assured that the cost remains reasonable\.

On the implementation side, we demonstrate the power of our methods by training thermodynamic models on image classification benchmarks \([Section4](https://arxiv.org/html/2607.00170#S4)\)\. When running our models with thermodynamic inference, we achieve accuracies of 98\.1% on MNIST, 93\.5% on FashionMNIST, 94\.9% on CIFAR\-10, and 76% on CIFAR\-100\. We also demonstrate with selected results that there is a mild tradeoff between inference\-time cost and performance: more extreme parameters and larger sweep counts can achieve better accuracy, but less extreme inference costs still produce good results\. Analyzing the models, we show that, for our larger models, over99\.99%99\.99\\%of the FLOPs have been off\-loaded to thermodynamic inference\. In[Section6](https://arxiv.org/html/2607.00170#S6), we continue the experimental investigation by validating our theoretical models, visually characterizing model performance, conducting ablation experiments to show the impact of our methods, and explaining the procedure for deriving optimal inference schedules\.

Finally, in[Section7](https://arxiv.org/html/2607.00170#S7), we discuss how our results relate to the current thermodynamic hardware landscape\. Any particular hardware system imposes unique constraints, especially on sparsity and quantization; therefore, we discuss how our algorithms can fit into that picture with an experiment\. Furthermore, while the new approach of high\-temperature sampling may be able to inject a new ingredient into the long\-standing relationship between Ising computation and combinatorial optimization, ultimately thermodynamic hardware must be developed which is well\-suited to the AI workloads of the future\.

Our methods focus primarily on convolutional networks\. One reason is that they fit the theory—an unrolled conv with tanh activation is essentially a sparse MLP\. Another reason is that hardware requires sparsity and locality, and convolutional architectures build that in naturally from the beginning\. But also, we believe that convolutional networks are especially well suited to the workloads we need thermodynamic AI for\. Low\-power image processing would be extremely useful for small surveillance drones or security cameras\. Applied to one\-dimensional convolutions, the same approach could benefit implantable medical sensors, biometric wearables, and audio signal\-processing devices\.

To achieve our results, we use an unorthodox model of the Ising system as a computation platform\. Instead of annealing to zero temperature or sampling a single state from a Boltzmann distribution, as is traditional in Ising computation, we rely on empirical time\-averages of the Boltzmann distribution, taken by averaging over many samples at a fixed temperature\. By building the theory of high\-temperature computation into a full\-scale application, we show that there is merit to this new approach: high\-temperature average sampling is a natural fit for neural inference, and first\-class support for this sampling technique should be a priority in hardware development\.

### 1\.2Relation to Previous Work

This work is also related to neural training methods based on equilibrium or bidirectional dynamics, including equilibrium\-propagation\-style approaches such as\[E\_P\_2025\]\. Those methods demonstrate that deep networks can be trained to tolerate iterative bidirectional relaxation\. Our setting is different in two important respects\. First, the communicating activations are stochastic binary spins rather than continuous states\. Second, we deliberately avoid making the main training procedure depend on long sampling\-based gradient estimators\. The Gibbs chain supplies targets for regularization, but the gradient path remains short, smooth, and compatible with standard backpropagation\.

The closest benchmark comparison is\[Niazi\_2024\], which demonstrates image classification using Gibbs\-sampled Ising systems and provides an important baseline for the field\. Our approach differs by using a more conventional deep\-learning training pipeline, by targeting deeper convolutional architectures, and by emphasizing the sweep\-count tradeoff needed for larger models\. The comparison is not merely one of accuracy: the goal is to show that thermodynamic inference can be made to look like a scalable neural software stack rather than a collection of separately engineered small Ising classifiers\.

On the hardware side, the architecture proposed in\[jelinčič2025efficientprobabilistichardwarearchitecture\]is especially relevant\. That work articulates a system\-level picture in which probabilistic sampling units communicate through binary states and are supplemented by a small amount of conventional computation at the input and output\. Our model is closely aligned with that picture: it uses a classical encoder and decoder around a sequence of Ising\-compatible thermodynamic blocks\. The contribution here is the training and inference recipe needed to make such blocks perform deep neural computation at nontrivial scale\.

Finally, this work should be distinguished from the broader Ising\-machine literature on combinatorial optimization\. Quantum annealers, CMOS annealers, digital annealers, coherent Ising machines, simulated bifurcation machines, and probabilistic\-bit systems all provide important evidence that Ising\-like substrates can be engineered and scaled\. However, most of that literature targets optimization: construct an energy function and search for a low\-energy state\. The present paper uses Ising dynamics differently\. We do not anneal to a ground state\. We use finite\-temperature Gibbs dynamics to estimate Boltzmann averages, and we train a neural architecture so that the signs of those averages implement feed\-forward inference\.

## 2Architecture

![Refer to caption](https://arxiv.org/html/2607.00170v1/figures/arch.drawio.png)Figure 1:Our CIFAR\-10 classifier\. The encoder is a single classical convolution with sign activation\. This is followed by six thermodynamic blocks\. The classical decoder is a global average pool followed by a fully\-connected linear map to the output layer\. Each thermodynamic block contains three convolutional layers, as depicted in[Figure2](https://arxiv.org/html/2607.00170#S2.F2)\.Our model is loosely based on the Wide Residual Network \(WRN\)\[zagoruyko2017wideresidualnetworks\], a classic and well\-tested convolutional architecture\. It consists of an encoder, a sequence of thermodynamic blocks, and a decoder\. The thermodynamic blocks account for between99%99\\%and over99\.99%99\.99\\%of the FLOPs\.111To speak precisely, theyreplacewhat would have been FLOPs in classical inference\.Each thermodynamic block is composed of two tanh\-activated convolutional layers with a skip\-convolution, followed by sign\-activation to ensure that data flowing between blocks is binary\. The model can be run as a normal feed\-forward model on a GPU\. On the other hand, each thermodynamic block can also be loaded onto an Ising machine and run with Gibbs sampling: the neurons are implemented as spins, and the weights and biases are simply reinterpreted as the coupling strengths and local magnetic fields, respectively\. We envision the hardware implementation of the model as follows: classical GPU computation or dedicated digital ASICs perform the lightweight encoder and decoder\. Each thermodynamic block is implemented by an Ising machine chip\. Binary information is extracted from a block and passed to the next block as input\. This setup is closely aligned with the overall hardware architecture proposed in\[jelinčič2025efficientprobabilistichardwarearchitecture\], although we use a completely different software model and training system\. A visual representation of our model architecture is shown in[Figure1](https://arxiv.org/html/2607.00170#S2.F1)\.

### 2\.1The Thermodynamic Block in Feed\-Forward Mode

When run on a GPU in feed\-forward mode, a thermodynamic block is a simple composition of convolutions, affine transforms, and hyperbolic tangent activations, with a final output passed through sign activation\. There is nothing intrinsically probabilistic or energy\-based about the block when run in this way\. However, its architecture is carefully chosen so that it is capable of mapping directly onto probabilistic and energy\-based Ising hardware\. Each block is parameterized by three convolution operatorsK1,K2,K3K\_\{1\},K\_\{2\},K\_\{3\}, one input affine map

Ain​\(x\)=ain⊙x\+din,A\_\{\\mathrm\{in\}\}\(x\)=a\_\{\\mathrm\{in\}\}\\odot x\+d\_\{\\mathrm\{in\}\},one internal affine map

Amid​\(u\)=amid⊙u\+dmid,A\_\{\\mathrm\{mid\}\}\(u\)=a\_\{\\mathrm\{mid\}\}\\odot u\+d\_\{\\mathrm\{mid\}\},and a hidden biasbb\. In our models, the convolutionsK1K\_\{1\}andK2K\_\{2\}use a3×33\\times 3kernel, while the skip\-convolutionK3K\_\{3\}uses a1×11\\times 1kernel\. Given a binary input tensorxx, the sign\-forward block computation of output binary tensoryyis

x′\\displaystyle x^\{\\prime\}=Ain​\(x\),\\displaystyle=A\_\{\\mathrm\{in\}\}\(x\),\(1\)z1\\displaystyle z\_\{1\}=K1​x′\+b,\\displaystyle=K\_\{1\}x^\{\\prime\}\+b,\(2\)z2\\displaystyle z\_\{2\}=K2​Amid​tanh⁡\(z1\)\+K3​x′,\\displaystyle=K\_\{2\}A\_\{\\mathrm\{mid\}\}\\tanh\(z\_\{1\}\)\+K\_\{3\}x^\{\\prime\},\(3\)y\\displaystyle y=sign⁡\(z2\)\.\\displaystyle=\\operatorname\{sign\}\(z\_\{2\}\)\.\(4\)Thus the internal state of the block is continuous, but the quantity exported to the next block is binary\. The continuous internal states correspond approximately to thermal mean\-field averages of the block when run on an Ising machine\. The exported output must be binary to satisfy the requirements of inter\-chip communication with Ising hardware, but this binarization also serves another purpose: clamping outputs to binary rounds away much of the error introduced by mapping the block into an Ising machine\.

### 2\.2The Thermodynamic Block as an Ising Machine

thermodynamic sampler𝒙\(i\)∈\{±1\}Cin×H×W\\bm\{x\}^\{\(i\)\}\\in\\\{\{\\pm\}1\\\}^\{C\_\{\\text\{in\}\}\\times H\\times W\}𝒔0\\bm\{s\}\_\{0\}pinned spinsK1∘Ain\\text\{K\}\_\{1\}\\ \\circ A\_\{\\text\{in\}\}𝒔1\\bm\{s\}\_\{1\}bbδ⋅K2∘Amid\\delta\\cdot\\text\{K\}\_\{2\}\\circ A\_\{\\text\{mid\}\}δ⋅K3∘Ain\\delta\\cdot\\text\{K\}\_\{3\}\\circ A\_\{\\text\{in\}\}𝒔2\\bm\{s\}\_\{2\}𝒂=sign⁡⟨𝒔2⟩\\bm\{a\}=\\operatorname\{sign\}\\langle\\bm\{s\}\_\{2\}\\rangle𝒚\(i\)∈\{±1\}Cout×H′×W′\\bm\{y\}^\{\(i\)\}\\in\\\{\{\\pm\}1\\\}^\{C\_\{\\text\{out\}\}\\times H^\{\\prime\}\\times W^\{\\prime\}\}Figure 2:Architecture of a single thermodynamic block, implemented as an Ising system\. It contains three layers of spins: input spinss0s\_\{0\}, hidden spinss1s\_\{1\}, and output spinss2s\_\{2\}\. The input layers0s\_\{0\}is pinned in place\. They are connected by bidirectional linear couplings, depicted in purple boxes\. Couplings in an Ising system are bidirectional, but the labels indicate the forward coupling, i\.e\. moving downward in the diagram\. The reverse couplings are the linear adjoint; for the exact relationship, see[Equation5](https://arxiv.org/html/2607.00170#S2.E5)\. The hidden layer is also affected by the local magnetic fieldsbb\. Finally, an average\-and\-sign circuitaaaccumulates time\-averages ofs2s\_\{2\}and outputs their signs\.In[Section2\.1](https://arxiv.org/html/2607.00170#S2.SS1), we describe the feed\-forward operation of the thermodynamic block\. In this section, we will describe how the parameters can be mapped directly to an Ising system\. A critical extra ingredient is added to perform this mapping: a real\-valued output coupling scaleδ∈\(0,1\]\\delta\\in\(0,1\], which we call thetemperature gradient parameter\. The parameterδ\\deltamultiplies all couplings and biases applied to the output layer\. The name ‘temperature gradient’ comes from an analogy: modifying the temperature effectively scales all the couplings at once through changes inβ\\beta\. Therefore, scaling down all the couplings on one side of the system is analogous to increasing the temperature on that side\. Thus, a smaller value ofδ\\deltacorresponds to a stronger ‘temperature gradient\.’ Smaller values ofδ\\deltaallow the interaction between the input layer and the hidden layer to be dominant over the interaction between the hidden layer and the output layer, encouraging forward information flow at the cost of increased noise\.

To perform the mapping, we begin by observing that the affine and convolutional operations compose to form a network of sparse linear couplings, where the activation sites are identified as Ising spins\. The three layers of activation sites are interpreted as three layers of Ising spins: the input spinss0s\_\{0\}, the hidden spinss1s\_\{1\}, and the output spinss2s\_\{2\}\. Written in inner\-product notation, theHamiltonian\(energy function\) of the unit is

H​\(s0,s1,s2\)=−⟨s1,K1​ain​s0⟩−δ​⟨s2,K3​ain​s0⟩−δ​⟨s2,K2​amid​s1⟩\\displaystyle H\(s\_\{0\},s\_\{1\},s\_\{2\}\)=\-\\langle s\_\{1\},K\_\{1\}a\_\{\\text\{in\}\}s\_\{0\}\\rangle\-\\delta\\langle s\_\{2\},K\_\{3\}a\_\{\\text\{in\}\}s\_\{0\}\\rangle\-\\delta\\langle s\_\{2\},K\_\{2\}a\_\{\\text\{mid\}\}s\_\{1\}\\rangle\(5\)−⟨s1,b\+K1​din⟩−δ​⟨s2,K3​din\+K2​dmid⟩,\\displaystyle\-\\langle s\_\{1\},b\+K\_\{1\}d\_\{\\text\{in\}\}\\rangle\-\\delta\\langle s\_\{2\},K\_\{3\}d\_\{\\text\{in\}\}\+K\_\{2\}d\_\{\\text\{mid\}\}\\rangle,\(6\)where the couplings are on the first line and the biases on the second line\. We assume that we are able to fix the input spinss0s\_\{0\}, so that they do not evolve over time\. This allows us to feed stable input data into the block\. On the other hand, the values of the hidden spinss1s\_\{1\}and the output spinss2s\_\{2\}are allowed to evolve\. We also fix the inverse temperatureβ\\betaat 1, the natural value given the coupling scale set by neural\-network training\. The statistical mechanics of Ising systems ensures that the values of these layers will be distributed according to a Boltzmann distribution on the energy functionHH\. LettingΣmid=\{±1\}C×H×W\\Sigma\_\{\\text\{mid\}\}=\\\{\\pm 1\\\}^\{C\\times H\\times W\}andΣout=\{±1\}C×H×W\\Sigma\_\{\\text\{out\}\}=\\\{\\pm 1\\\}^\{C\\times H\\times W\}be the configuration spaces of the hidden and output layers, whereCCis the number of channels andH×WH\\times Wis the image size, we can write the Boltzmann distribution as

ℙ​\[s1,s2\|s0\]=1Z​\[s0\]​exp⁡\(−H​\(s0,s1,s2\)\),\\displaystyle\\mathbb\{P\}\[s\_\{1\},s\_\{2\}\|s\_\{0\}\]=\\frac\{1\}\{Z\[s\_\{0\}\]\}\\exp\(\-H\(s\_\{0\},s\_\{1\},s\_\{2\}\)\),Z​\[s0\]=∑s1∈Σmid∑s2∈Σoutexp⁡\(−H​\(s0,s1,s2\)\)\.\\displaystyle Z\[s\_\{0\}\]=\\sum\_\{s\_\{1\}\\in\\Sigma\_\{\\text\{mid\}\}\}\\sum\_\{s\_\{2\}\\in\\Sigma\_\{\\text\{out\}\}\}\\exp\(\-H\(s\_\{0\},s\_\{1\},s\_\{2\}\)\)\.\(7\)
In order to operate the block, we set the input spinss0s\_\{0\}to our input binary tensorxx, then allow the system to fluctuate according to the Boltzmann distribution\. We observe the long\-term time average of the output spins, the conditional expectation⟨s2\|s0⟩\\langle s\_\{2\}\|s\_\{0\}\\rangle, according to this distribution\. Finally, we extract the output by applying the sign function to this average, i\.e\.y=sign​⟨s2\|s0⟩y=\\operatorname\{sign\}\\langle s\_\{2\}\|s\_\{0\}\\rangle\. The block architecture is diagrammed in[Figure2](https://arxiv.org/html/2607.00170#S2.F2)\. The reason that it is designed in this way is that previous theory\[moore2025correspondenceisingmachinesneural\]guarantees that asδ\\deltaapproaches0,⟨s1\|s0⟩\\langle s\_\{1\}\|s\_\{0\}\\rangleapproachestanh⁡\(z1\)\\tanh\(z\_\{1\}\)andsign​⟨s2\|s0⟩\\operatorname\{sign\}\\langle s\_\{2\}\|s\_\{0\}\\rangleapproachessign⁡\(z2\)\\operatorname\{sign\}\(z\_\{2\}\)\. In other words, for a strong enough temperature gradient parameter, the thermodynamic computing system closely approximates the result of the block when operated in feed\-forward mode\.

It should be noted that this approach to using an Ising system for computation differs significantly from those typical in the existing literature\. Typically, Ising systems are used in one of two ways\. Inannealing, the temperature is initially high, then decreased slowly as the system fluctuates, bringing the system toward a low\-energy state\. After the anneal, the state of the system is observed, extracting a candidate low\-energy state\. On the other hand, in Boltzmann\-style sampling, the temperature is fixed, but the object of interest is the probability distribution itself: the output is extracted by observing an output\-spin configuration at a sampling time, thus sampling directly from the Boltzmann distribution\. In our setup, we observeBoltzmann averages, which can be thought of as averaging over many samples from a Boltzmann machine\. This approach, while non\-traditional, does naturally fit into the existing hardware paradigm, and is necessary to make the correspondence with feed\-forward neural networks work\[moore2025correspondenceisingmachinesneural\]\.

### 2\.3Gibbs Sampling

So far, we have explained how the thermodynamic block works from the perspective of equilibrium thermodynamics\. However, we have not yet specified itsdynamics\. While any Markov chain whose stationary distribution is the Boltzmann distribution would work, for this paper we will focus on chromatic Gibbs sampling with two colors \(also called red\-black sampling\), corresponding to the hidden and output layers\. This is a natural approach because, with the input fixed, the dynamically changing component of the block is a sparse bipartite graph\. In order to run the block, we approximate the Boltzmann average ofs2s\_\{2\}using finitely many samples\. We will define an update of both layers to be a single ‘sweep’ of the sampler, and we will refer to the number of sweeps executed asGG\.

Thus, the Gibbs sampler provides the core inference algorithm by probabilistically evolving the spins in an Ising block\. We will now specify its operation precisely\. Letℬ​\(p\)\\mathcal\{B\}\(p\)be a variant of the Bernoulli distribution which takes11with probabilityppand−1\-1with probability1−p1\-p\. When applied to a vector, it is sampled element\-wise and independently\. Letσ\\sigmabe the sigmoid function, i\.e\.tanh\\tanhrescaled to have the range\(0,1\)\(0,1\)\. Then, the Gibbs sampling algorithm performs sweeps of the form

s1\\displaystyle s\_\{1\}←ℬ​\(σ​\(2​\(δ​amid⊙K2⊤​s2\+K1​Ain​x\+b\)\)\),\\displaystyle\\leftarrow\\mathcal\{B\}\(\\sigma\(2\(\\delta a\_\{\\text\{mid\}\}\\odot K\_\{2\}^\{\\top\}s\_\{2\}\+K\_\{1\}A\_\{\\text\{in\}\}x\+b\)\)\),\(8\)s2\\displaystyle s\_\{2\}←ℬ​\(σ​\(2​δ​\(K2​Amid​s1\+K3​Ain​x\)\)\)\.\\displaystyle\\leftarrow\\mathcal\{B\}\(\\sigma\(2\\delta\(K\_\{2\}A\_\{\\text\{mid\}\}s\_\{1\}\+K\_\{3\}A\_\{\\text\{in\}\}x\)\)\)\.\(9\)Outputs are extracted by time\-averaging the values ofs2s\_\{2\}, then applying the sign function\. There is also a burn\-in period ofGbG\_\{b\}sweeps in which the values are discarded\. Usually,Gb=⌊G/4⌋G\_\{b\}=\\lfloor G/4\\rfloor, though it is sometimes useful to useGb=0G\_\{b\}=0\. In short, the output, referred to as theGibbs consensus, is computed as

y=sign⁡\(1G−Gb​∑g=Gb\+1Gs2\(g\)\)\.y=\\operatorname\{sign\}\\left\(\\frac\{1\}\{G\-G\_\{b\}\}\\sum\_\{g=G\_\{b\}\+1\}^\{G\}s\_\{2\}^\{\(g\)\}\\right\)\.\(10\)wheres2\(g\)s\_\{2\}^\{\(g\)\}refers to the state ofs2s\_\{2\}after sweepgg\. There is no backpropagation through this procedure: it is used only for inference, and to generate teacher tensors for the Gibbs regularizer discussed in[Section3](https://arxiv.org/html/2607.00170#S3)\.

The values ofs1s\_\{1\}ands2s\_\{2\}are initialized differently depending on what we are using the Gibbs sampler for\. In a real device, and therefore in our real inference, we use random i\.i\.d\. Rademacher initialization\. However, when we use Gibbs sampling in the training regularizer, this extra ingredient of realism actually just wastes computation\. Instead, it is advantageous to start as close as possible to equilibrium in order to minimize the value ofGGnecessary to get a strong signal of the true equilibrium\. This can be done with feed\-forward initialization, which is described in[Section3\.2](https://arxiv.org/html/2607.00170#S3.SS2)\.

This setup introduces the core problem of inference hyperparameter selection: given a trained model, we must selectGGandδ\\deltato maximize accuracy under Gibbs sampling while minimizing sweep cost\. Generally speaking, higherGGgives better accuracy at any choice ofδ\\delta, while lowerδ\\deltagives better long\-term accuracy but requires higher values ofGGto achieve it\. This hyperparameter selection problem is critical to extracting good performance from our trained models; we will analyze this problem thoroughly and provide heuristics, modeling, and theory for hyperparameter selection in[Section5\.1](https://arxiv.org/html/2607.00170#S5.SS1)and[Section6](https://arxiv.org/html/2607.00170#S6)\. For now, the important qualitative point is this: while the theory in\[moore2025correspondenceisingmachinesneural\]promises that asδ→0\\delta\\rightarrow 0andG→∞G\\rightarrow\\inftythe block’s outputs under Gibbs sampling will approach the block’s outputs under feed\-forward evaluation, itdoes notmake any promises about the actual values ofδ\\deltaandGGrequired\. In fact, without specialized training, obtaining useful accuracy from the thermodynamic block would require extremely lowδ\\delta, and, as a consequence, impractically highGG\. Training the network to perform well at more moderateδ\\deltaand reasonably low sweep countGGis the purpose of our core training algorithm, to which we now turn\.

## 3Gibbs\-Regularization Training

The full training procedure has four phases\. The purpose of the first three phases, which we refer to as pretraining, is to turn a standard residual network into a network which is capable of being loaded onto an Ising machine—in other words, we create the blocks described in[Section2\.1](https://arxiv.org/html/2607.00170#S2.SS1)\. A detailed description of the pretraining pipeline can be found in[AppendixA](https://arxiv.org/html/2607.00170#A1)\. The output of the pretraining is a trained Straight\-Through Estimator \(STE\)\. The purpose of the fourth phase, which is our core algorithmic contribution, is to reshape that network’s weights so that its feed\-forward outputs remain stable under the noise and backward coupling of the Ising system\. The goal is to create a network which not only performs well in feed\-forward operation, but also produces good accuracy at reasonable cost under Gibbs sampling\. While the raw STE is technically a valid Ising system, as proven by the correspondence theorem, the purpose of this training is to turn a theoretical solution into a practically useful solution\. To achieve this, we continue to train the STE while introducing two new regularizers, thefixed\-point loss\([Section3\.1](https://arxiv.org/html/2607.00170#S3.SS1)\) and themagnitude regularizer\([Section3\.4](https://arxiv.org/html/2607.00170#S3.SS4)\)\. Taken together, we call this processGibbs regularization\.

The goal of Gibbs regularization is to reshape the weights so that the STE activations of each block become stable fixed points, or near\-fixed points, of the corresponding Gibbs dynamics\. To achieve this goal, we run a new cosine\-annealed STE training schedule while introducing our new regularizers according to a gentle curriculum\. The overall loss function is

ℒ=\(1−α\)​ℒCE\+α​ℒKD\+λFP​ℒFP\+λmag​ℒmag\\mathcal\{L\}=\(1\-\\alpha\)\\,\\mathcal\{L\}\_\{\\text\{CE\}\}\+\\alpha\\,\\mathcal\{L\}\_\{\\text\{KD\}\}\+\\lambda\_\{\\text\{FP\}\}\\,\\mathcal\{L\}\_\{\\text\{FP\}\}\+\\lambda\_\{\\text\{mag\}\}\\,\\mathcal\{L\}\_\{\\text\{mag\}\}\(11\)The cross\-entropy classification lossℒCE\\mathcal\{L\}\_\{\\text\{CE\}\}and the knowledge distillation lossℒKD\\mathcal\{L\}\_\{\\text\{KD\}\}are the same as in pure STE training\. The STE path carries the classification signal\.222As a variation, some of the classification pathway can be sent through theygibbsy\_\{\\text\{gibbs\}\}instead of the STE classification\. This can help the decoder adapt to the Gibbs\-sampled states it will actually receive, but it makes the STE training more difficult\. This also has the interesting effect of yielding Gibbs\-sampled accuracies which outstrip feed\-forward accuracies; this is not generally expected otherwise\.However, we have added two new loss terms, the Gibbs fixed\-point lossℒFP\\mathcal\{L\}\_\{\\text\{FP\}\}and the magnitude regularizerℒmag\\mathcal\{L\}\_\{\\text\{mag\}\}, which we will describe shortly\. Three curriculum schedules prevent the Gibbs stability objective from destabilizing training:

1. \(a\)Backward coupling ramp:We introduce a parameterqqwhich linearly increases from 0 to 1\.0 over the first 10 epochs\. This parameter modulates the backward coupling in the calculation of the fixed\-point loss, allowing the model to slowly adjust to the influence of backward couplings\.
2. \(b\)Gibbs sweep ramp:The number of Gibbs sweeps per block increases from 2 toGtrainG\_\{\\text\{train\}\}over 10 epochs\. Early epochs see only short chains, which differ less from the feed\-forward initialization\. Later epochs require stability under longer Gibbs chains\.
3. \(c\)FP loss weight ramp:λFP\\lambda\_\{\\text\{FP\}\}ramps from 0 to its target value over 3 epochs, preventing the fixed\-point loss from dominating before the model has begun adapting\.

We also use gradient clipping at5\.05\.0\. Data augmentation includes random crops, horizontal flips, AutoAugment, and Cutout\-16, though the exact augmentation protocol should be dataset\-specific\. The batch size is set to 64\. We have observed that larger batch sizes result in much worse training performance, suggesting that the natural stochasticity of small batch sizes may be important for proper adaptation\.

### 3\.1The Fixed\-Point Loss

The fixed\-point loss is the mechanism by which the Gibbs equilibrium informs the gradient update\. It consists of two terms per block, measuring whether the sign\-forward outputs are self\-consistent under backward coupling from the Gibbs consensus\. We begin by sweeping each block forGtrainG\_\{\\text\{train\}\}sweeps, then extracting the empirical averagea=⟨s2⟩Gtraina=\\langle s\_\{2\}\\rangle\_\{G\_\{\\text\{train\}\}\}and settingygibbs=sign⁡ay\_\{\\text\{gibbs\}\}=\\operatorname\{sign\}a\. To reduce noise, the initial state for Gibbs sampling is usually computed using the feed\-forward activations, as described in[Section3\.2](https://arxiv.org/html/2607.00170#S3.SS2)\. The valueygibbsy\_\{\\text\{gibbs\}\}does not have any gradients, but it can be used as a target value to penalize the network for drifting far from its feed\-forward state under Gibbs sampling\. It gives the model an empirical target for the block’s equilibrium behavior\.

The total fixed\-point loss for the block isℒFP=ℒFP\(1\)\+ℒFP\(2\)\\mathcal\{L\}\_\{\\text\{FP\}\}=\\mathcal\{L\}\_\{\\text\{FP\}\}^\{\(1\)\}\+\\mathcal\{L\}\_\{\\text\{FP\}\}^\{\(2\)\}\. The two terms measure Gibbs consistency on the hidden and output layers, respectively, and are defined as follows\.

The first part of the fixed\-point loss is, roughly speaking, the error between the feed\-forward hidden layer and a rough mean\-field estimate for the hidden layer factoring in the backward coupling from an output layer set equal toygibbsy\_\{\\text\{gibbs\}\}\. It checks whether the sign\-forward hidden state remains self\-consistent once backward coupling from the Gibbs consensus is introduced\. We define it as

z1\\displaystyle z\_\{1\}=b\+K1​Ain​x,\\displaystyle=b\+K\_\{1\}A\_\{\\text\{in\}\}x,\(12\)z^1\\displaystyle\\hat\{z\}\_\{1\}=q​δ⋅amid⊙K2⊤​ygibbs,\\displaystyle=q\\delta\\cdot a\_\{\\text\{mid\}\}\\odot K\_\{2\}^\{\\top\}y\_\{\\text\{gibbs\}\},\(13\)ℒFP\(1\)\\displaystyle\\mathcal\{L\}\_\{\\text\{FP\}\}^\{\(1\)\}=‖sign⁡z1−tanh⁡\(z1\+z^1\)‖2\.\\displaystyle=\\left\\\|\\operatorname\{sign\}z\_\{1\}\-\\tanh\(z\_\{1\}\+\\hat\{z\}\_\{1\}\)\\right\\\|^\{2\}\.\(14\)This formula may be somewhat puzzling to readers, who might expect‖tanh⁡z1−tanh⁡\(z1\+q​δ​amid​K2⊤​a\)‖2\\\|\\tanh z\_\{1\}\-\\tanh\(z\_\{1\}\+q\\delta a\_\{\\text\{mid\}\}K\_\{2\}^\{\\top\}a\)\\\|^\{2\}\. However, this more obvious form does not work as well in practice\. Usingsign⁡a\\operatorname\{sign\}ainstead ofaaandsign⁡z1\\operatorname\{sign\}z\_\{1\}instead oftanh⁡z1\\tanh z\_\{1\}produces much better results\. Empirically, the signed target produces a stronger stability signal than the corresponding tanh target\. We interpret this as an implicit margin\-enlarging objective: the loss penalizes not only disagreement with the Gibbs consensus, but also small output fields that would be vulnerable to finite\-sweep fluctuations\. Additionally, the stronger loss signal improves training speed\. We should also note that setting the target to be a sign rather than a tanh also implicitly helps the magnitude regularizer, as it tries to push the spin averages estimated by the feed\-forward pre\-activation closer to±1\\pm 1\.

The second term propagates that coupled hidden estimate forward and checks whether the resulting output remains aligned with the original feed\-forward output\. We assume that the mean\-field value of the hidden layer is in facttanh⁡\(z1\+z^1\)\\tanh\(z\_\{1\}\+\\hat\{z\}\_\{1\}\)and propagate it to the output layer, checking it against the inference\-mode STE forward output\. We set

z2\\displaystyle z\_\{2\}=K3​Ain​x\+K2​Amid​tanh⁡z1,\\displaystyle=K\_\{3\}A\_\{\\text\{in\}\}x\+K\_\{2\}A\_\{\\text\{mid\}\}\\tanh z\_\{1\},\(15\)z^2\\displaystyle\\hat\{z\}\_\{2\}=K3​Ain​x\+K2​Amid​tanh⁡\(z1\+z^1\),\\displaystyle=K\_\{3\}A\_\{\\text\{in\}\}x\+K\_\{2\}A\_\{\\text\{mid\}\}\\tanh\(z\_\{1\}\+\\hat\{z\}\_\{1\}\),\(16\)ℒFP\(2\)\\displaystyle\\mathcal\{L\}\_\{\\text\{FP\}\}^\{\(2\)\}=‖sign⁡z2−tanh⁡\(δ​z^2\)‖2\.\\displaystyle=\\left\\\|\\operatorname\{sign\}z\_\{2\}\-\\tanh\(\\delta\\hat\{z\}\_\{2\}\)\\right\\\|^\{2\}\.\(17\)Again, the reader may have expected‖tanh⁡\(δ​z2\)−tanh⁡\(δ​z^2\)‖2\\\|\\tanh\(\\delta z\_\{2\}\)\-\\tanh\(\\delta\\hat\{z\}\_\{2\}\)\\\|^\{2\}, and again we have found that this more obvious form does not perform well in practice\. Usingsign⁡z2\\operatorname\{sign\}z\_\{2\}increases the loss, exaggerating it in the correct direction, and encouraging spin\-averages further from zero\. This extra magnitude regularization ends up being helpful in practice\.

Note that all the gradients flow through the tanh, as thesign⁡z1\\operatorname\{sign\}z\_\{1\}andsign⁡z2\\operatorname\{sign\}z\_\{2\}terms have no gradients\. Note also thatygibbsy\_\{\\text\{gibbs\}\}has no gradient, so this design yields clean and low\-variance gradients with a short backpropagation tree: the stochastic Gibbs chain provides the*target*\(what the equilibrium looks like\), while only smooth, differentiable operations contribute to the gradient\.

### 3\.2Feed\-Forward Spin Initialization

In the preceding section, we did not describe how the Markov chain is initialized\. As mentioned in[Section2\.3](https://arxiv.org/html/2607.00170#S2.SS3), in the inference mode we use the hardware\-realistic initialization of i\.i\.d\. Rademacher\. However, for generating the fixed\-point loss, it is more computationally efficient to initialize as close as possible to the equilibrium state, since we can get reliable low\-noise estimates of the long\-term Gibbs consensus with far fewer sweeps\. Therefore, when using the Gibbs sampler in training, we initialize the spins using the feed\-forward activations of the block, setting

s1\(0\)\\displaystyle s\_\{1\}^\{\(0\)\}=sign⁡\(b\+K1​Ain​x\),\\displaystyle=\\operatorname\{sign\}\(b\+K\_\{1\}A\_\{\\text\{in\}\}x\),\(18\)s2\(0\)\\displaystyle s\_\{2\}^\{\(0\)\}=sign⁡\(K3​Ain​x\+K2​Amid​tanh⁡\(b\+K1​Ain​x\)\)\.\\displaystyle=\\operatorname\{sign\}\(K\_\{3\}A\_\{\\text\{in\}\}x\+K\_\{2\}A\_\{\\text\{mid\}\}\\tanh\(b\+K\_\{1\}A\_\{\\text\{in\}\}x\)\)\.\(19\)This gets the system close enough to equilibrium to provide a good training signal\. In practice, with this initialization evenGtrain=10G\_\{\\text\{train\}\}=10is sufficient to provide good stability information for training, with only diminishing returns up toGtrain=40G\_\{\\text\{train\}\}=40\.

As a variation, we sometimes do use random Gibbs initialization, as opposed to feed\-forward, and/or eliminate burn\-in by settingGb=0G\_\{b\}=0\. The former encourages convergence toward the desired state from a broader set of initial conditions, while the latter increases the sample size at the cost of some noise\. These modifications may help in some cases, but generally come at the cost of increasingGtrainG\_\{\\text\{train\}\}and thus slowing down training\.

### 3\.3Training\-Time Temperature Gradient

The fixed\-point loss described in[Section3\.1](https://arxiv.org/html/2607.00170#S3.SS1)takes the temperature gradient parameterδ\\deltaas a hyperparameter\. Naively, one might assume that we should train with the same value ofδ\\deltawith which we intend to perform inference\. However, this is not the case\. In general, we will perform very few sweeps during training \(say,Gtrain=10G\_\{\\text\{train\}\}=10\), and many more sweeps during inference\. We want to train with aδ\\deltathat is high enough that lowGGgives good information about the equilibrium state; lowGGand lowδ\\deltawill produce aygibbsy\_\{\\text\{gibbs\}\}which is essentially random noise\. On the other hand, we will perform inference with much lower values ofδ\\deltaand higher values ofGGto extract high classification performance\.

Additionally, we have observed that networks trained with a particular value ofδtrain\\delta\_\{\\text\{train\}\}perform well onδ<δtrain\\delta<\\delta\_\{\\text\{train\}\}, moderately well atδ≈δtrain\\delta\\approx\\delta\_\{\\text\{train\}\}, and very poorly forδ\>δtrain\\delta\>\\delta\_\{\\text\{train\}\}\. It follows that if we want a network capable of flexibly operating across different regimes, we need a higherδtrain\\delta\_\{\\text\{train\}\}\. However,δtrain=1\\delta\_\{\\text\{train\}\}=1is suboptimal for training as well:δtrain\\delta\_\{\\text\{train\}\}needs to be low enough that the Gibbs equilibrium is similar to the Gibbs equilibrium in the low\-δ\\deltainference regime\. Through experimentation, we have determined thatδtrain=0\.7\\delta\_\{\\text\{train\}\}=0\.7is usually a good compromise\.

### 3\.4Magnitude Regularizer

Pre\-activations which are too far from zero or too close to zero are both problematic\. Pre\-activations too close to zero require large numbers of sweeps for convergence \(see[Section5](https://arxiv.org/html/2607.00170#S5)\)\. When small pre\-activations from the input layer occur in the hidden layer, the spin becomes more susceptible to being changed by backward coupling from the output layer, which is also undesirable\. Pre\-activations which are too far from zero result in longer autocorrelation times, and thus also in more sweeps required for convergence\. Additionally, large pre\-activations saturate the gradient of tanh, resulting in slower learning\. A bilateral penalty on pre\-activation magnitudes prevents pathological field strengths\.ℒmag\\mathcal\{L\}\_\{\\text\{mag\}\}is calculated as the average over all spins of

ReLU​\(m2−\|z\|\)2\+0\.1⋅ReLU​\(\|z\|−3​m\)2,\\text\{ReLU\}\\big\(\\tfrac\{m\}\{2\}\-\|z\|\\big\)^\{2\}\+0\.1\\cdot\\text\{ReLU\}\\big\(\|z\|\-3m\\big\)^\{2\},\(20\)wherezzis the feed\-forward pre\-activation of the spin andmmis the target magnitude, a training hyperparameter which we usually set to32\\frac\{3\}\{2\}\.

## 4Results

In[Table1](https://arxiv.org/html/2607.00170#S4.T1)we present selected results for the basic classification datasets MNIST, FashionMNIST, CIFAR\-10, and CIFAR\-100\. The table shows classification accuracy under thermodynamic inference for each dataset\. We also show several variations to demonstrate how cost and accuracy can be traded off\.

Table 1:Results on major image classification datasets\. ‘Accuracy’ refers to the model’s classification performance under thermodynamic operation, i\.e\. with each thermodynamic block run according to the Gibbs sampling algorithm, while ‘FF Accuracy’ refers to the model’s classification accuracy in feed\-forward mode\. ‘Parameters’ denotes the number of trainable parameters in the feed\-forward STE model, while ‘% Classical’ denotes the percentage of FLOPs in the feed\-forward model not offloaded to Ising systems\.Two things may jump out to the reader about this table: first of all, our better results often come simply from better feed\-forward accuracy\. This is a result of weakeningλFP\\lambda\_\{\\text\{FP\}\}\. This change increases feed\-forward accuracy, but requires more sweeps and lowerδ\\deltato get close to it\. The other thing that might surprise the reader is the fact that FMNIST accuracy exceeds the feed\-forward accuracy\. This is due to the use of a slightly different version of the training algorithm using algorithm tricks described as variations in[Section3](https://arxiv.org/html/2607.00170#S3)\.

We include a measurement of “% Classical” for each model\. This refers to the number of classical FLOPs required for thermodynamic inference as a percentage of the FLOPs required to run the model in classical feed\-forward mode\. Even in the smallest model, over 99% of the FLOPs have been off\-loaded to the thermodynamic hardware, whereas for the larger models essentially all of the computation is thermodynamic\. For our CIFAR\-10 model, featured in most of our analytical experiments in[Section6](https://arxiv.org/html/2607.00170#S6)and[Section5](https://arxiv.org/html/2607.00170#S5), 99\.99% of the computation is offloaded to the thermodynamic units\. This shows that if thermodynamic hardware achieves significant power\-consumption reductions, our model could take advantage of them\.

### 4\.1Parameter Counting

The parameter numbers we listed in[Table1](https://arxiv.org/html/2607.00170#S4.T1)count the number of classical parameters in the feed\-forward version of the model\. This refers to the number of parameters that need to be trained\. It is not the same as the number of spins required to implement the model using Ising chips, nor the number of couplings present on\-chip\. Implementing the thermodynamic block on an Ising machineunrollsthe convolutions into sparse linear connections and uses a spin for each channel of each pixel\. The block’s affine transforms are folded into the linear couplings and biases\. The number of biases for each block is equal to the number of spins\. Performing the accounting for our 31M\-parameter CIFAR\-10 model, we obtain the total spin and coupling counts shown in[Table2](https://arxiv.org/html/2607.00170#S4.T2)\. Note that these counts refer to the model before any quantization or sparsification is applied\. Sparsification techniques can significantly compress the number of spin\-spin couplings required\. We discuss this in more detail in[Section7\.2](https://arxiv.org/html/2607.00170#S7.SS2)\.

Table 2:Total number of spins and spin\-spin couplings for each thermodynamic block in the 6\-block CIFAR\-10 classifier model, before any post\-training sparsification techniques\.

## 5Theoretical Model

Previous theory\[moore2025correspondenceisingmachinesneural\]proves that low enoughδ\\deltaand high enoughGGwill give good performance\. In[Section3](https://arxiv.org/html/2607.00170#S3), we presented algorithms designed to improve performance at moderateδ\\deltaandGG\. However, in order to practically guide experiments and inference, we need to know how performance actually scales with these two parameters: for a certain performance target, what is the scaling behavior ofGGasδ→0\\delta\\rightarrow 0? And how does the error rate behave as sweeps progress? This is a critical practical question, as a thermodynamic computing system which takes forever to compute is not much use to anyone\. Broadly speaking, it turns out that the error rate is a simple power law in the sweep number, that the critical prefactor is the density of pre\-activations near zero, and that sweep cost is roughly proportional to the inverse square ofδ\\delta\.

In the following sections, we will answer these questions rigorously and precisify these heuristic statements\. We will precisely pin down the sources of error, provide approximate formulae for calculating how error relates toδ\\deltaandGG, and demonstrate that the mixing times of our thermodynamic blocks are controllable\. We will build an analytical model that allows us to predict and control error rates, and exhibit experiments demonstrating its power\. The key observation in the following analysis is that our systems are quite large and weakly coupled—therefore, simply ignoring spin\-spin correlations and assuming spatial independence allows us to predict system\-scale statistical quantities with a high degree of accuracy\. The deviations from this assumption can be cleanly identified by comparing theory and practice; as it turns out, our theory is accurate precisely where it needs to be to find the optimal inference hyperparameters for high classification performance\.

In[Section5\.1](https://arxiv.org/html/2607.00170#S5.SS1), we investigate and define the error sources in our model\. In[Section5\.2](https://arxiv.org/html/2607.00170#S5.SS2), we provide an effective and simple model using Gaussian mixtures for the distribution of activations, which is the critical quantity in statistical error estimation\. In[Section5\.3](https://arxiv.org/html/2607.00170#S5.SS3), we use this model to estimate the effect of choosingGGusing numerical integration, then provide a simplified closed\-form mathematical solution\. Finally, in[AppendixC](https://arxiv.org/html/2607.00170#A3), we show that we can keep the mixing times under control, tightening the validity of the rest of the theory\. Experimental investigations applying this theory to the analysis and tuning of our models follow in[Section6](https://arxiv.org/html/2607.00170#S6)\.

### 5\.1Sources of Error

Before we can analyze the thermodynamic performance parameters, it is necessary to discuss conceptually the taxonomy of error sources and their precise definition\. Thermodynamic operation introduces errors relative to feed\-forward operation at each block\. Therefore, understanding the etiology of these errors and keeping them under control is vital to obtaining good inference performance by intelligent selection ofGGandδ\\delta\.

As we have discussed, a thermodynamic block operating on an Ising machineapproximatesits feed\-forward STE output\. In general, when our models are run in thermodynamic mode by running each block with time\-averaged Gibbs sampling, the classification accuracy of the network approximates its feed\-forward classification performance\. The lower the value ofδ\\deltaand the higher the value ofGG, the closer the thermodynamically sampled accuracy comes to the feed\-forward accuracy\. In general, however, the feed\-forward accuracy is an upper bound to the thermodynamic performance, seeing as errors introduced by Gibbs sampling are essentially random noise and therefore extremely unlikely to help classification accuracy rather than simply introduce entropy\.333It should be noted that algorithmic variations such as introducingygibbsy\_\{\\text\{gibbs\}\}path into the cross\-entropy signal can violate this assumption\. However, we will focus on the simpler case where the STE training only uses the feed\-forward signal\.Examples of this phenomenon can be see in our main benchmark results,[Table1](https://arxiv.org/html/2607.00170#S4.T1)\. We observe, indeed, that the thermodynamic accuracy is slightly lower than the feed\-forward accuracy\.

There are three parts to the error analysis\. At the level of the binary tensors passed between blocks, we need to understand the errors introduced by the choice ofGG, and the errors introduced by the choice ofδ\\delta\. At the level of the model as a whole, we need to understand the statistical relationship between these error sources and the rate of misclassification\. The first two factors, which comprise the underlying sources of error, can be analyzed mathematically\. The third factor, which tackles how the per\-spin error rates translate to decreases classification accuracy, is fundamentally empirical, and will be measured and modeled in[Section6\.2](https://arxiv.org/html/2607.00170#S6.SS2),[Section6\.3](https://arxiv.org/html/2607.00170#S6.SS3), and[Section6\.4](https://arxiv.org/html/2607.00170#S6.SS4)\. For now, we turn to defining the underlying error sources precisely\.

#### 5\.1\.1Stochastic Error

Stochastic error is the error introduced by the finiteness ofGG\. For fixedδ\\delta, asG→∞G\\rightarrow\\infty, we know that⟨s⟩G,δ→μδ\\langle s\\rangle\_\{G,\\delta\}\\rightarrow\\mu\_\{\\delta\}, where⟨⋅⟩G,δ\\langle\\cdot\\rangle\_\{G,\\delta\}indicates empirical average afterGGsweeps andμδ\\mu\_\{\\delta\}indicates the true Boltzmann average ofss\. This is a well\-known consequence of the theory of Markov chains\. For finiteGG, there will generally be a discrepancy:⟨s⟩G,δ≠μδ\\langle s\\rangle\_\{G,\\delta\}\\neq\\mu\_\{\\delta\}\. In the output layer, if the sign is not changed, the error is tolerated\. If the sign is changed, an error is introduced into the block output\. This source of error is measured by thestochastic error rate:

R\(G,δ\)=1\|s2\|∑i=1\|s2\|𝟏\{sign⟨s2,i⟩G,δ≠signμi,δ\}\.\\displaystyle R\(G,\\delta\)=\\frac\{1\}\{\|s\_\{2\}\|\}\\sum\_\{i=1\}^\{\|s\_\{2\}\|\}\\bm\{1\}\\\{\\operatorname\{sign\}\\langle s\_\{2,i\}\\rangle\_\{G,\\delta\}\\neq\\operatorname\{sign\}\\mu\_\{i,\\delta\}\\\}\.\(21\)In other words, the stochastic error rateRRof a particular block is the fraction of spins in the output layer of that block whose finite time averages differ in sign from their true long\-term averages\. We will see in[Section5\.3\.1](https://arxiv.org/html/2607.00170#S5.SS3.SSS1)that

R​\(G,δ\)∼Cδ​G,\\displaystyle R\(G,\\delta\)\\sim\\frac\{C\}\{\\delta\\sqrt\{G\}\},\(22\)whereCCis a constant dependent only on the weights of the block and the input\. Practically speaking, stochastic error is combated by increasing the sweep countGG\.

#### 5\.1\.2Thermal Error

Thermal error is the error introduced by the choice ofδ\\delta\. By the correspondence theorem\[moore2025correspondenceisingmachinesneural\], we know thatμδ/δ→z2\\mu\_\{\\delta\}/\\delta\\rightarrow z\_\{2\}asδ→0\\delta\\rightarrow 0, and thereforeygibbs=sign⁡\(μδ\)→sign⁡\(z2\)=ySTEy\_\{\\text\{gibbs\}\}=\\operatorname\{sign\}\(\\mu\_\{\\delta\}\)\\rightarrow\\operatorname\{sign\}\(z\_\{2\}\)=y\_\{\\text\{STE\}\}whenever the corresponding component ofz2z\_\{2\}is nonzero, whereySTEy\_\{\\text\{STE\}\}is the output of the block in feed\-forward mode\. Due to the effect of backward couplings, the Boltzmann spin average of the output spinss2s\_\{2\}, which we refer to asμδ\\mu\_\{\\delta\}, will differ fromtanh⁡\(δ​z2\)\\tanh\(\\delta z\_\{2\}\)\. If it does not differ in sign, the discrepancy is tolerated\. If it does differ in sign, an error is introduced into the block output\. Error introduced by the discrepancy at finiteδ\\deltais measured by thethermal error rate:

E\(δ\)=1\|s2\|∑i=1\|s2\|𝟏\{signμi,δ≠sign\(z2\)i\},\\displaystyle E\(\\delta\)=\\frac\{1\}\{\|s\_\{2\}\|\}\\sum\_\{i=1\}^\{\|s\_\{2\}\|\}\\bm\{1\}\\\{\\operatorname\{sign\}\\mu\_\{i,\\delta\}\\neq\\operatorname\{sign\}\(z\_\{2\}\)\_\{i\}\\\},\(23\)whereμi,δ\\mu\_\{i,\\delta\}is the Boltzmann average of spinii\. This measures the fraction of spins whose sign\-outputs will not converge to their feed\-forward values no matter how long the system is sampled\. Thermal error is combated by lowering the temperature gradientδ\\delta\. Larger values ofλFP\\lambda\_\{\\text\{FP\}\}result in lower thermal error at any given value ofδ\\delta\. Generally, forδ\>δtrain\\delta\>\\delta\_\{\\text\{train\}\}, thermal error is high enough to cause accuracy collapse\.

#### 5\.1\.3Overall Error

Thermodynamic sampling introduces random error on a per\-block basis in two ways, stochastic and thermal error\. At the scale of the whole model, errors in the internal activations translate statistically into errors in the output\. Because errors are introduced at each block, the total error rate between the thermodynamically operated model and the feed\-forward operated model increases with depth as errors accumulate\. Due to this accumulation, deeper models require a combination of increased per\-block error tolerance, lowerδ\\delta, and higher sweep count: increasing tolerance lowers the accuracy penalty of errors, decreasingδ\\deltalowers thermal error, and increasingGGdecreases stochastic error\.

### 5\.2Modeling the Output Spin Distribution

![Refer to caption](https://arxiv.org/html/2607.00170v1/figures/gaussian_mixture.png)Figure 3:Empirical\|z\|\|z\|distribution and 3\-component Gaussian mixture fit for the first block of our CIFAR\-10 model\.In an Ising block, we have a very large number of output spins\. We can calculate per\-spin statistical quantities\. However, the behavior of our systems is governed by a large\-scale aggregation over these individual quantities\. The aggregate behavior can differ qualitatively from the behavior of any individual component\. For instance, we observe that each spin converges exponentially, but the overall error rate converges according to a power law\. Given the very large number of spins, and the fact that the statistical behavior of each spin depends only on its long\-term average activation, it is both justified and mathematically convenient to assume that a block contains an infinite number of spins governed by an absolutely continuous distribution\.

We begin by analyzing the distributionρ​\(\|z\|\)\\rho\(\|z\|\)of absolute pre\-activations of output spins in the feed\-forward model\. Empirically, this can be measured by using a histogram\. In[Figure3](https://arxiv.org/html/2607.00170#S5.F3)we show a real example of a\|z\|\|z\|\-histogram\. This shape is very typical\. The peak is caused by the magnitude regularizer; almost all blocks trained with our method look like this\. The next step is to choose an analytical approximationρ^\\hat\{\\rho\}of the distributionρ\\rhothat fits well and is convenient to work with\. A Gaussian mixture is effective, natural, and renders the analysis quite easy \(empirically,N=3N=3mixture components is enough for a good fit\)\. It is most important thatρ^\\hat\{\\rho\}fitsρ\\rhowell nearz=0z=0, as this region turns out to govern the asymptotics of the stochastic error rate, but qualitative global fit is still desirable\. To preserve behavior nearz=0z=0, when fitting the Gaussian mixture we constrain one of the spikes to be centered onz=0z=0and require thatρ^​\(0\)=ρ​\(0\)\\hat\{\\rho\}\(0\)=\\rho\(0\)\. To maintain a good global fit in theL1L^\{1\}sense, we perform the mixture fit in linear space, rather than log\-probability space: asymptotic errors in the tail distribution are simply not important\.

We will introduce some notation to refer to our Gaussian mixtures\. LetNNbe the number of mixture components\. Each componentkk, for1≤k≤N1\\leq k\\leq N, has a weightwkw\_\{k\}, a standard deviationσk\\sigma\_\{k\}, and a meanμk\\mu\_\{k\}\. Our approximation is therefore

ρ^​\(x\)=∑k=1Nwk​ρ^\(k\)​\(x\)=∑k=1NwkZk​exp⁡\(−\(x−μk\)22​σk2\)\\displaystyle\\hat\{\\rho\}\(x\)=\\sum\_\{k=1\}^\{N\}w\_\{k\}\\hat\{\\rho\}^\{\(k\)\}\(x\)=\\sum\_\{k=1\}^\{N\}\\frac\{w\_\{k\}\}\{Z\_\{k\}\}\\exp\\left\(\-\\frac\{\(x\-\\mu\_\{k\}\)^\{2\}\}\{2\\sigma\_\{k\}^\{2\}\}\\right\)\(24\)where we use the truncated Gaussian normalizer

Zk=∫0∞exp⁡\(−\(z−μk\)22​σk2\)​𝑑z=π2​σk​erfc⁡\(−μkσk​2\)\\displaystyle Z\_\{k\}=\\int\_\{0\}^\{\\infty\}\\exp\\left\(\-\\frac\{\(z\-\\mu\_\{k\}\)^\{2\}\}\{2\\sigma\_\{k\}^\{2\}\}\\right\)\\,dz=\\sqrt\{\\frac\{\\pi\}\{2\}\}\\sigma\_\{k\}\\operatorname\{erfc\}\\left\(\-\\frac\{\\mu\_\{k\}\}\{\\sigma\_\{k\}\\sqrt\{2\}\}\\right\)\(25\)
#### 5\.2\.1Factoring in Thermal Error

As we will see in[Section6\.2](https://arxiv.org/html/2607.00170#S6.SS2), thermal error can be effectively modeled as an independent zero\-mean Gaussian perturbation of the output\-spin pre\-activations\. Empirical modeling shows that the variance of the perturbationsϵδ2\\epsilon^\{2\}\_\{\\delta\}is governed by a power law,CT​δαTC\_\{T\}\\delta^\{\\alpha\_\{T\}\}\. The Gaussian mixture approximation for the feed\-forward pre\-activations makes it very easy to estimate the distribution of spin averages under thermal error perturbation: it is simplytanh∗⁡ρ^δ\\tanh\_\{\*\}\\hat\{\\rho\}\_\{\\delta\}, whereρ^δ\(k\)∼𝒩​\(μk,σk2\+CT​δαT\)\\hat\{\\rho\}^\{\(k\)\}\_\{\\delta\}\\sim\\mathcal\{N\}\(\\mu\_\{k\},\\sigma^\{2\}\_\{k\}\+C\_\{T\}\\delta^\{\\alpha\_\{T\}\}\)\. For computing error rates at finiteδ\\delta, we will use the blurred ‘pre\-activation’ distributionρ^δ\\hat\{\\rho\}\_\{\\delta\}, which is also a Gaussian mixture\. We will also use the notationσk,δ2=σk2\+CT​δαT\\sigma\_\{k,\\delta\}^\{2\}=\\sigma\_\{k\}^\{2\}\+C\_\{T\}\\delta^\{\\alpha\_\{T\}\}for brevity\.

### 5\.3Modeling the Stochastic Error Rate

Let us first consider the probability of error on a single output spinσi\\sigma\_\{i\}\. After discarding burn\-in, writeGGfor the number of retained sweeps and let

⟨σi⟩G:=1G​∑t=1Gσi\(t\),\\displaystyle\\langle\\sigma\_\{i\}\\rangle\_\{G\}:=\\frac\{1\}\{G\}\\sum\_\{t=1\}^\{G\}\\sigma\_\{i\}^\{\(t\)\},\(26\)whereσi\(t\)\\sigma\_\{i\}^\{\(t\)\}denotes the value of the spiniiat sweeptt, and letμi\\mu\_\{i\}denote the equilibrium mean ofσi\\sigma\_\{i\}under the Boltzmann distribution, i\.e\.⟨σi⟩∞\\langle\\sigma\_\{i\}\\rangle\_\{\\infty\}\. An error occurs whensign⟨σi⟩G≠sign\(μi\)\\operatorname\{sign\}\\langle\\sigma\_\{i\}\\rangle\_\{G\}\\neq\\operatorname\{sign\}\(\\mu\_\{i\}\)\. Sinceσi∈\{−1,\+1\}\\sigma\_\{i\}\\in\\\{\-1,\+1\\\}, its equilibrium variance is1−μi21\-\\mu\_\{i\}^\{2\}\. We use the MCMC convention for the integrated autocorrelation time,

τi,int:=1\+2​∑ℓ=1∞Cov⁡\(σi\(0\),σi\(ℓ\)\)Var⁡\(σi\)\.\\displaystyle\\tau\_\{i,\\mathrm\{int\}\}:=1\+2\\sum\_\{\\ell=1\}^\{\\infty\}\\frac\{\\operatorname\{Cov\}\(\\sigma\_\{i\}^\{\(0\)\},\\sigma\_\{i\}^\{\(\\ell\)\}\)\}\{\\operatorname\{Var\}\(\\sigma\_\{i\}\)\}\.\(27\)With this convention, the Markov chain central limit theorem gives

G\(⟨σi⟩G−μi\)\\ext@arrow0099\\arrowfill@\-\-⟶𝒟𝒩\(0,\(1−μi2\)τi,int\),\\displaystyle\\sqrt\{G\}\\left\(\\langle\\sigma\_\{i\}\\rangle\_\{G\}\-\\mu\_\{i\}\\right\)\\ext@arrow 0099\\arrowfill@\\relbar\\relbar\\longrightarrow\{\}\{\\mathcal\{D\}\}\\mathcal\{N\}\\\!\\left\(0,\\,\(1\-\\mu\_\{i\}^\{2\}\)\\tau\_\{i,\\mathrm\{int\}\}\\right\),\(28\)or equivalently, for largeGG,

⟨σi⟩G≈𝒩​\(μi,1G​\(1−μi2\)​τi,int\)\.\\displaystyle\\langle\\sigma\_\{i\}\\rangle\_\{G\}\\approx\\mathcal\{N\}\\\!\\left\(\\mu\_\{i\},\\,\\frac\{1\}\{G\}\(1\-\\mu\_\{i\}^\{2\}\)\\tau\_\{i,\\mathrm\{int\}\}\\right\)\.\(29\)Therefore, by the Markov chain central limit theorem,

ℙ\[sign⟨σi⟩G≠sign\(μi\)\]≈12erfc\(\|μi\|G2​\(1−μi2\)​τi,int\)\.\\displaystyle\\mathbb\{P\}\\\!\\left\[\\operatorname\{sign\}\\langle\\sigma\_\{i\}\\rangle\_\{G\}\\neq\\operatorname\{sign\}\(\\mu\_\{i\}\)\\right\]\\approx\\frac\{1\}\{2\}\\operatorname\{erfc\}\\\!\\left\(\|\\mu\_\{i\}\|\\sqrt\{\\frac\{G\}\{2\(1\-\\mu\_\{i\}^\{2\}\)\\tau\_\{i,\\mathrm\{int\}\}\}\}\\right\)\.\(30\)
#### 5\.3\.1Asymptotic Total Error Rate

![Refer to caption](https://arxiv.org/html/2607.00170v1/figures/R_prediction.png)Figure 4:Measured stochastic error rate compared against numerically evaluated integral for the first block of our CIFAR\-10 model\. The fit is very close, but the integral does not capture the period\-two oscillation caused by the chromatic Gibbs sampling\.We have calculated the stochastic error rate for a single output spin as a function of its true mean and integrated autocorrelation time\. Now, we extend our calculation of the stochastic error rate to the entire system, calculatingRRby integrating our per\-spin stochastic error rate over the distributionρ\\rho\. We will assume for simplicity that all spins have the same integrated autocorrelation timeτ\\tau\. Then, we will start by writing the integral approximation,

𝔼​\[R​\(G,δ\)\]≈12​∫0∞erfc​\(G2​τ​tanh⁡\(δ​z\)1−tanh\(δz\)2\)​ρ^δ​\(z\)​𝑑z=12​∫0∞erfc​\(θ​sinh⁡\(δ​z\)\)​ρ^δ​\(z\)​𝑑z,\\displaystyle\\mathbb\{E\}\[R\(G,\\delta\)\]\\approx\\frac\{1\}\{2\}\\int\_\{0\}^\{\\infty\}\\text\{erfc\}\\left\(\\sqrt\{\\frac\{G\}\{2\\tau\}\}\\frac\{\\tanh\(\\delta z\)\}\{\\sqrt\{1\-\\tanh\(\\delta z\)^\{2\}\}\}\\right\)\\hat\{\\rho\}\_\{\\delta\}\(z\)\\,dz=\\frac\{1\}\{2\}\\int\_\{0\}^\{\\infty\}\\text\{erfc\}\(\\theta\\sinh\(\\delta z\)\)\\hat\{\\rho\}\_\{\\delta\}\(z\)\\,dz,\(31\)where we make the substitution

θ:=G2​τ\.\\displaystyle\\theta:=\\sqrt\{\\frac\{G\}\{2\\tau\}\}\.\(32\)An experimental validation, plotted in[Figure4](https://arxiv.org/html/2607.00170#S5.F4), shows that this integral is an excellent approximation of the real stochastic error rate—though it cannot capture the period\-two oscillation in error rates, which is an artifact of the two\-color Gibbs sampling, a feature which is not detected by the CLT\-style error estimate\. At this point, we can easily find the asymptotic approximation of the tail for anyδ\\delta: By substitution, we get

12​δ​θ​∫0∞erfc​\(s\)​ρ^δ​\(1δ​arsinh​\(sθ\)\)​d​s1\+s2/θ2=1π​∫0∞e−s2​P^δ​\(1δ​arsinh​\(sθ\)\)​𝑑s\\displaystyle\\frac\{1\}\{2\\delta\\theta\}\\int\_\{0\}^\{\\infty\}\\text\{erfc\}\(s\)\\hat\{\\rho\}\_\{\\delta\}\\left\(\\frac\{1\}\{\\delta\}\\text\{arsinh\}\\left\(\\frac\{s\}\{\\theta\}\\right\)\\right\)\\frac\{ds\}\{\\sqrt\{1\+s^\{2\}/\\theta^\{2\}\}\}=\\frac\{1\}\{\\sqrt\{\\pi\}\}\\int\_\{0\}^\{\\infty\}e^\{\-s^\{2\}\}\\hat\{P\}\_\{\\delta\}\\left\(\\frac\{1\}\{\\delta\}\\text\{arsinh\}\\left\(\\frac\{s\}\{\\theta\}\\right\)\\right\)\\,ds\(33\)whereP^δ\\hat\{P\}\_\{\\delta\}is the CDF ofρ^δ\\hat\{\\rho\}\_\{\\delta\}\. Lettingu=δ​θu=\\delta\\theta, Taylor expansion then gives us the leading\-order asymptotic formula

𝔼​R=ρ^δ​\(0\)2​u​π\+𝒪​\(u−2\)\.\\displaystyle\\mathbb\{E\}R=\\frac\{\\hat\{\\rho\}\_\{\\delta\}\(0\)\}\{2u\\sqrt\{\\pi\}\}\+\\mathcal\{O\}\(u^\{\-2\}\)\.\(34\)Ifρ^δ\\hat\{\\rho\}\_\{\\delta\}has an even differentiable extension through the origin, thenρ^δ′​\(0\)=0\\hat\{\\rho\}\_\{\\delta\}^\{\\prime\}\(0\)=0and the remainder improves to𝒪​\(u−3\)\\mathcal\{O\}\(u^\{\-3\}\)\. The primary asymptotic quantity isu=δ​θu=\\delta\\theta, which is proportional toδ​G\\delta\\sqrt\{G\}, obtaining the first\-order power law approximation we promised in[Section5\.1\.1](https://arxiv.org/html/2607.00170#S5.SS1.SSS1)\. Note that whileρ^δ​\(0\)\\hat\{\\rho\}\_\{\\delta\}\(0\)does depend onδ\\delta, it is only weakly dependent at small values ofδ\\delta, and this effect can be ignored for the purposes of understandingδ→0\\delta\\rightarrow 0behavior\.

#### 5\.3\.2Low\-δ\\deltaClosed Form

We would like to understand how fast the rateRRconverges to its leading power law\. Introducing the magic approximation\[LETHER1993573\]

erfc⁡\(x\)≈e−q​\(x\),\\displaystyle\\operatorname\{erfc\}\(x\)\\approx e^\{\-q\(x\)\},q​\(x\)=1623​x2\+2π​x,\\displaystyle q\(x\)=\\frac\{16\}\{23\}x^\{2\}\+\\frac\{2\}\{\\sqrt\{\\pi\}\}x,\(35\)and using the fact that forδ<14\\delta<\\frac\{1\}\{4\}we can essentially ignore thesinh\\sinhand obtain a nearly tight upper bound, we can approximate the integral from[Equation31](https://arxiv.org/html/2607.00170#S5.E31)as

∫0∞erfc​\(θ​sinh⁡\(δ​z\)\)​ρ^δ​\(z\)​𝑑z≤∫0∞erfc​\(u​z\)​ρ^δ​\(z\)​𝑑z\\displaystyle\\int\_\{0\}^\{\\infty\}\\text\{erfc\}\(\\theta\\sinh\(\\delta z\)\)\\hat\{\\rho\}\_\{\\delta\}\(z\)\\,dz\\leq\\int\_\{0\}^\{\\infty\}\\text\{erfc\}\\left\(uz\\right\)\\hat\{\\rho\}\_\{\\delta\}\(z\)\\,dz\(36\)≈∑k=1NwkZk,δ​∫0∞exp⁡\(−1623​u2​z2−2π​u​z−\(z−μk\)22​σk,δ2\)​𝑑z\.\\displaystyle\\approx\\sum\_\{k=1\}^\{N\}\\frac\{w\_\{k\}\}\{Z\_\{k,\\delta\}\}\\int\_\{0\}^\{\\infty\}\\exp\\left\(\-\\frac\{16\}\{23\}u^\{2\}z^\{2\}\-\\frac\{2\}\{\\sqrt\{\\pi\}\}uz\-\\frac\{\(z\-\\mu\_\{k\}\)^\{2\}\}\{2\\sigma\_\{k,\\delta\}^\{2\}\}\\right\)dz\.\(37\)This is now a straightforward, if messy, application of the error function\. By introducing the following shorthands,

ak​\(u\):=6423​u2\+2σk,δ2,\\displaystyle a\_\{k\}\(u\):=\\frac\{64\}\{23\}u^\{2\}\+\\frac\{2\}\{\\sigma\_\{k,\\delta\}^\{2\}\},bk​\(u\):=2π​u−μkσk,δ2,\\displaystyle b\_\{k\}\(u\):=\\frac\{2\}\{\\sqrt\{\\pi\}\}u\-\\frac\{\\mu\_\{k\}\}\{\\sigma\_\{k,\\delta\}^\{2\}\},\(38\)we can evaluate the integral for the Gaussian mixture componentkkas

πak​exp⁡\(−μk22​σk,δ2\)​erfcx⁡\(bkak\)\.\\sqrt\{\\frac\{\\pi\}\{a\_\{k\}\}\}\\exp\\left\(\-\\frac\{\\mu\_\{k\}^\{2\}\}\{2\\sigma\_\{k,\\delta\}^\{2\}\}\\right\)\\operatorname\{erfcx\}\\left\(\\frac\{b\_\{k\}\}\{\\sqrt\{a\_\{k\}\}\}\\right\)\.\(39\)We then obtain a closed\-form approximation for𝔼​R\\mathbb\{E\}R,

𝔼​R≈∑k=1Nwk​Rk,\\displaystyle\\mathbb\{E\}R\\approx\\sum\_\{k=1\}^\{N\}w\_\{k\}R\_\{k\},Rk=ρ^δ\(k\)​\(0\)​π16​ak​erfcx⁡\(bkak\)\.\\displaystyle R\_\{k\}=\\hat\{\\rho\}\_\{\\delta\}^\{\(k\)\}\(0\)\\sqrt\{\\frac\{\\pi\}\{16a\_\{k\}\}\}\\operatorname\{erfcx\}\\left\(\\frac\{b\_\{k\}\}\{\\sqrt\{a\_\{k\}\}\}\\right\)\.\(40\)The large\-uubehavior ofRkR\_\{k\}is governed byak−1/2a\_\{k\}^\{\-1/2\}, which is asymptotically proportional tou−1u^\{\-1\}\. The erfcx provides a small\-uucorrection factor, and becomes akk\-independent constant at largeuu\. It should be noted that the normalization ensures thatRk​\(0\)=12R\_\{k\}\(0\)=\\frac\{1\}\{2\}\.

#### 5\.3\.3Convergence Rate

![Refer to caption](https://arxiv.org/html/2607.00170v1/figures/spike.png)Figure 5:Theoretical convergence rate for a single Gaussian spike away fromz=0z=0, with time of transition to power\-law convergence visualized\.The closed\-form solution to the integral allows us to inspect analytically how each spike in the distributionρ\\rhocontributes to the convergence rate\. The error rate contribution for a single Gaussian spike, as a function ofuu, follows a predictable structure: For very smalluu, the error rate decays from1/21/2to1/51/5nearly linearly\. Then, asuugrows, the error rate decays roughly exponentially, as the bulk of the erfc kernel passes through the spike, and we see very fast convergence\. At some point a crossover is reached where the erfc kernel is concentrated around the intersection of the spike with theyyaxis, and past this point, the error rate stabilizes and decays like a1/u1/upower law\. On a log\-log plot, this results in a sharp fall\-off in error which suddenly stabilizes into a straight\-line descent\. This behavior, along with the example spike in question, is depicted in the graph in[Figure5](https://arxiv.org/html/2607.00170#S5.F5)\. In the case that the spike is centered on zero, there is no exponential\-decay phase, and the error rate becomes a power law immediately\.

We would like to know how largeuumust be before theRRcurve converges to its long\-term power law behavior\. Since we can approximately decomposeρ\\rhointo a sum of spikes, we can answer this question independently for each Gaussian spike\. The key to understanding the cliff in the log\-log plot is the contribution of the erfcx function: for largeuuit converges to a constant, but for smalluuit can contribute a massive pre\-factor\. This is because fors\>0s\>0,erfcx⁡\(−s\)≈2​es2\\operatorname\{erfcx\}\(\-s\)\\approx 2e^\{s^\{2\}\}: thus, whenb​\(u\)<0b\(u\)<0, the erfcx term is decaying nearly exponentially\. It follows, then, that this decay stops whenb=0b=0, i\.e\.

u=μ​π2​σ2\.\\displaystyle u=\\frac\{\\mu\\sqrt\{\\pi\}\}\{2\\sigma^\{2\}\}\.\(41\)This implies that the contribution to the error rate of a Gaussian spike begins to closely approximate the power law when

G≥τ​μ2​π2​σ4​δ2\.\\displaystyle G\\geq\\frac\{\\tau\\mu^\{2\}\\pi\}\{2\\sigma^\{4\}\\delta^\{2\}\}\.\(42\)Beyond this point, the error rate contribution of the spike converges to the power law more slowly, with relative error orderu−1u^\{\-1\}\. Before this point, the error drops off much more quickly with increasingGG\. This crossover point can be quite large for a concentrated spike located away from the origin\. From this, we can derive the overall qualitative behavior ofRRin terms of the spike components ofρ^\\hat\{\\rho\}: for very smalluu, each spike contributes nearly equally\. Asuugrows, however, concentrated spikes far from the origin decay very quickly, while wide spikes near the origin decay much more slowly\. In other words, the concentrated spike far from the origin quickly ceases to matter relative to the wide bulk near the origin\. Eventually, all spikes contribute a power law, but the power\-law contribution is proportional to the spike’s mass atz=0z=0; it follows that concentrated spikes far from the origin contribute almost nothing to the long\-term power law rate\. This means that the bulk of the error rate is contributed by the spikes with large intersection with theyy\-axis: the bulk spike centered atz=0z=0is predominant, as it has no exponential\-decay phase and is a power law from the start, and any spike withσ\>μ\\sigma\>\\mucontributes significantly\.

## 6Experiments

Now that we have developed theoretical tools for analyzing our thermodynamic blocks, we will apply them to analytical experiments demonstrating how to characterize and tune the trained blocks\. We begin by showing the results of naive hyperparameter sweeping, exhibiting a heatmap of classification accuracy as a function of uniformly chosenGGandδ\\delta, in[Section6\.1](https://arxiv.org/html/2607.00170#S6.SS1)\. This should give the reader a good qualitative picture of the capabilities of our thermodynamic models\. Then, we move on to finer analysis\. In[Section6\.2](https://arxiv.org/html/2607.00170#S6.SS2), we fit the Gaussian power law model described in[Section5\.2\.1](https://arxiv.org/html/2607.00170#S5.SS2.SSS1)to experimentally measured thermal error\. Then, in[Section6\.3](https://arxiv.org/html/2607.00170#S6.SS3), we use our model of stochastic error to experimentally measure the power\-law relationship between block\-level stochastic error and overall classification accuracy\. Then, we combine our analytical tools and modeled quantities in[Section6\.4](https://arxiv.org/html/2607.00170#S6.SS4), demonstrating how we can predict model performance as a function ofδ\\deltaandRR, and how a desired error rate can be converted to an implementable sweep\-count schedule\. Finally, we investigate the importance of the Gibbs regularization algorithm in shaping the performance heatmaps by performing an ablation experiment in[Section6\.5](https://arxiv.org/html/2607.00170#S6.SS5)\.

### 6\.1Uniform\-Sweep Inference Performance

![Refer to caption](https://arxiv.org/html/2607.00170v1/figures/exp3.png)Figure 6:Gibbs accuracy on a\(δ,G\)\(\\delta,G\)grid with uniformGGper block, random initialization\. The color map is centered on the STE baseline, so cells rendered in dark blue match STE accuracy and cells in dark red collapse to random\. The dashed ridge marks the accuracy\-optimalδ\\deltaat eachGG; it drifts from largeδ\\deltaat smallGG\(short\-chain regime\) toward smallδ\\deltaat largeGG\(equilibrium regime\)\.Before we engage in more sophisticated methods for selectingGGgivenδ\\delta, it is worth seeing how the two hyperparameters determine classification accuracy directly\. The simplest possible inference schedule is to setGGuniformly for all blocks\. To visualize the qualitative effects of different choices ofGGandδ\\delta, we plot accuracy against these two variables in[Figure6](https://arxiv.org/html/2607.00170#S6.F6), for the case of our 6\-block CIFAR\-10 model\. This plot visually summarizes several of the key features that we have been discussing: increasing cost\-to\-converge asδ→0\\delta\\rightarrow 0, convergence to correctness forδ<δtrain\\delta<\\delta\_\{\\text\{train\}\}, increasing performance ceiling asδ→0\\delta\\rightarrow 0, and poor performance forδ\>δtrain\\delta\>\\delta\_\{\\text\{train\}\}orGGnot large enough\. We would like to observe that mid\-rangeδ\\deltacan give surprisingly good performance at very low sweep count\. An accuracy of89%89\\%can be achieved withδ=0\.45\\delta=0\.45using only 102 total Gibbs sweeps\. This shows that if inference\-cost minimization is the highest priority, we can still get pretty good results at very low sweep cost; the corresponding power consumption will depend on the hardware implementation\.

#### 6\.1\.1Explaining the High\-δ\\deltaOverfit

One feature of[Figure6](https://arxiv.org/html/2607.00170#S6.F6)which immediately draws attention and hasn’t been addressed by the theory discussed so far is the right\-hand side of the diagram: forδ\>δtrain\\delta\>\\delta\_\{\\text\{train\}\}, as the Markov chain progresses we see first an increase in accuracy, then a decrease back to random\. It is probably not a coincidence that the peak accuracy in this regime occurs roughly atG=GtrainG=G\_\{\\text\{train\}\}\. In other words, the Gibbs chains here have ‘overfit’, i\.e\. they have learned only at theGGthey were trained for, rather than generalizing to the equilibrium state\. This is a satisfactory moral explanation, but it still does not mechanically explain this phenomenon, which is rather odd: how can a system reliably go from random, to close\-to\-correct, and then move to its real equilibrium, which is totally wrong? Somehow, both the random initial state and the final equilibrium are bad, but the system reliably passes through a ‘good’ region on the way from one to the other\.

Atδ=0\.75\\delta=0\.75, when we sweep long enough \(G\>1000G\>1000, say\) that the system is essentially converged, the accuracy is terrible\. It follows that this poor accuracy is caused by excessive thermal error\. Since the system produces better accuracy at 10 sweeps per block, it follows that the better accuracy iscausedby stochastic error\. In other words, at highδ\\delta, stochastic error can partially cancel out the effects of thermal error\. However, as we will see in[Section6\.4](https://arxiv.org/html/2607.00170#S6.SS4), this cannot be due to broad statistical effects, as injecting i\.i\.d\. thermal and stochastic error cannot replicate the effect\. In other words, this anomaly may be caused by spin\-spin correlation patterns, something which our theory cannot capture\. The fact that the effect occurs at highδ\\delta, and therefore higher coupling strengths, strongly supports this inference\. However, it should also be noted that this anomaly occurs away from the real region of interest: to get good performance, we want lowδ\\delta, and that is precisely where our independent\-spin theory functions as expected\.

### 6\.2Modeling Thermal Error

![Refer to caption](https://arxiv.org/html/2607.00170v1/figures/thermal_power_law.png)Figure 7:Experimentally measured variance of per\-spin thermal perturbation, as a function ofδ\\delta; a power law fits the data cleanly\. The dashed line indicates the estimated stochastic noise floor\.It is difficult to isolate thermal error because we cannot access the true spin averagesμi\\mu\_\{i\}\. However, we can just estimateμ\\muby using a very large number of Gibbs sweeps, so that the stochastic error is very small, and assess the differences betweenμi\\mu\_\{i\}andtanh⁡\(δ​z2\)\\tanh\(\\delta z\_\{2\}\), with the understanding that some noise has entered into the measurement\. While this relationship is quite complex and difficult to capture theoretically, it turns out that we have a pretty good approximation using independent Gaussian noise,

arctanh⁡\(μi,δ\)−δ​zi≈𝒩​\(0,ϵδ2\)\.\\displaystyle\\operatorname\{arctanh\}\(\\mu\_\{i,\\delta\}\)\-\\delta z\_\{i\}\\approx\\mathcal\{N\}\(0,\\epsilon\_\{\\delta\}^\{2\}\)\.\(43\)In other words, we can essentially model thermal error as a Gaussian perturbation of pre\-activations in the feed\-forward model\. This approximation is not perfect\. In practiceϵ\\epsilonisn’t actually uniform inzz, and the distribution isn’t actually Gaussian \(it has some non\-negligible excess kurtosis\), but this is a good enough approximate model\. We mostly care about thermal error when it affects spins whose pre\-activation is close to zero: these spins are the most likely to introduce sign errors due to thermal perturbation, while spins with strong pre\-activation are not likely to change sign under thermal perturbation\.

In order to use this model, we need to actually measureϵδ2\\epsilon^\{2\}\_\{\\delta\}\. We experimentally measure the pre\-activation perturbation defined in[Equation43](https://arxiv.org/html/2607.00170#S6.E43)for various values ofδ\\delta, specifically at spins where\|δ​zi\|<14\|\\delta z\_\{i\}\|<\\frac\{1\}\{4\}\. We then obtain the variance of the perturbation, and fit a power law\. The result, shown in[Figure7](https://arxiv.org/html/2607.00170#S6.F7), shows that a power\-law fit is very close over the measured range\. Using this power law, we can now inject simulated thermal error at anyδ\\deltaby perturbing pre\-activations according to random Gaussian noise\. This will be our model forsynthetic thermal errorand is used in the following experiments\.

### 6\.3Global Effects of Stochastic Error

![Refer to caption](https://arxiv.org/html/2607.00170v1/figures/exp1.png)Figure 8:Accuracy versus synthetic stochastic error at zero thermal error fits well to a power law\.It will be useful to see what the direct relationship between stochastic error and overall accuracy is, assuming thermal error is zero and ignoring the specifics of the Gibbs dynamics\. We will synthetically inject stochastic errorrrat each block and run inference in feed\-forward mode, then measure accuracy\. Naively, we could just flip spins i\.i\.d\. uniformly, but our analysis shows that stochastic error is not evenly distributed: it affects small\-mean spins more frequently, and the weighting depends onδ\\delta\. Suppose we have selectedδ\\deltaand the target error raterr\. Then, using binary search, we findγ\\gammasuch that

r=1\|s2\|​∑i=1\|s2\|12​erfc​\(γ​sinh⁡\(δ​\|zi\|\)\),\\displaystyle r=\\frac\{1\}\{\|s\_\{2\}\|\}\\sum\_\{i=1\}^\{\|s\_\{2\}\|\}\\frac\{1\}\{2\}\\text\{erfc\}\\left\(\\gamma\\sinh\(\\delta\|z\_\{i\}\|\)\\right\),\(44\)whereziz\_\{i\}is the feed\-forward pre\-activation of output spinii\. In the case thatδ=0\\delta=0, we use theδ→0\\delta\\rightarrow 0asymptotic limit by instead performing binary search on

r=1\|s2\|​∑i=1\|s2\|12​erfc​\(γ​\|zi\|\)\.\\displaystyle r=\\frac\{1\}\{\|s\_\{2\}\|\}\\sum\_\{i=1\}^\{\|s\_\{2\}\|\}\\frac\{1\}\{2\}\\text\{erfc\}\\left\(\\gamma\|z\_\{i\}\|\\right\)\.\(45\)In either case, once the rightγ​\(r\)\\gamma\(r\)is determined, the terms inside the sum are the probabilities that each spin should be flipped\. We can easily sample these flips independently\.

We perform an experiment on our CIFAR\-10 model by injecting stochastic error at rateRRinto the feed\-forward model at zero thermal error\. The results are depicted in[Figure8](https://arxiv.org/html/2607.00170#S6.F8)\. Uniform flipping is actually very optimistic: in fact, correctly weighted injected stochastic error degrades performance much more quickly than uniform error\. This leads us to the understanding that low\-\|μ\|\|\\mu\|output spins are disproportionately important for classification performance\. This explains the results in other experiments which show that accuracy is essentially worthless asGGincreases, up until a sharp phase transition to near\-perfect accuracy, which we infer must occur when the low\-\|μ\|\|\\mu\|spins begin to converge\. Interestingly, we see clear power\-law behavior relating the stochastic error rate to the accuracy gap at low values ofRR\. A linear regression shows that the accuracy loss is roughly proportional toR2\.55R^\{2\.55\}\.

### 6\.4Optimal Selection of Sweep Count

![Refer to caption](https://arxiv.org/html/2607.00170v1/figures/synthetic_heatmap.png)Figure 9:Theoretical performance of CIFAR\-10 classifier at varyingRRandδ\\delta, based on synthetic stochastic and thermal noise\.How do we select the number of sweeps to run our models with? Naively, we might use the same number of sweeps for each block, and simply increase the sweep count until we get good performance\. This is the approach depicted in[Figure6](https://arxiv.org/html/2607.00170#S6.F6)\. However, this is not efficient, nor is it sensitive to the problem: some blocks are more expensive to run than others, and ideally we would allocate more sweeps where they are needed, thus minimizing cost while maximizing performance\. Instead of specifyingGGupfront and seeing what performance comes out, we would like to specify performance up front and see whatGGcomes out\.

We can do exactly that using our models of stochastic and thermal error\. To begin with, we can use our idealized synthetic noise models to theoretically predict network performance at varying levels ofRRandδ\\delta\. To do this, we run the model in feed\-forward mode, injecting noise while we do so\. First, we inject synthetic thermal noise based onδ\\delta, perturbing the pre\-activations of the output layers based on[Equation43](https://arxiv.org/html/2607.00170#S6.E43)\. Then, we flip spins at a rateRRusing the relationship[Equation44](https://arxiv.org/html/2607.00170#S6.E44)\. We plot the resulting classification accuracy from this experiment in[Figure9](https://arxiv.org/html/2607.00170#S6.F9)\. This prediction captures the overall performance characteristics of the network, though it is not sensitive to more complex and subtle behavior like the large\-δ\\deltaoverfit: note that at largeδ\\delta, it fails to predict the bad\-good\-bad transition with decreasingRR\. This is likely due to the fact that the overfitting behavior is caused by underlying spin\-spin correlations and delicate structure of the energy function which cannot be captured by the independent\-spin synthetic noise model\.

Now, we will test these predicted accuracies with an experiment\. In order to do this, we need to convert a target value ofRRto a prescription for sweep count on each block\. First, we fit a Gaussian mixture model to the empirical histogram of pre\-activations, as described in[Section5\.2](https://arxiv.org/html/2607.00170#S5.SS2)\. Then, we computeG​\(R\)G\(R\)by performing binary search on the numerical integral[Equation31](https://arxiv.org/html/2607.00170#S5.E31), with the additional correction for the average integrated autocorrelation timeτ=1\+0\.71​δ2\\tau=1\+0\.71\\delta^\{2\}, based on our eigendecomposition theory \(see[AppendixC](https://arxiv.org/html/2607.00170#A3)\)\. This converts a pair\(δ,R\)\(\\delta,R\)to a sweep\-count schedule for each block\. The results of performing the experiment on our CIFAR\-10 model are plotted in[Figure10](https://arxiv.org/html/2607.00170#S6.F10)\. We show both the accuracy heatmap and a total sweep cost heatmap\. As expected, we see that lowerδ\\deltaand lowerRRresult in better performance and higher cost\.

![Refer to caption](https://arxiv.org/html/2607.00170v1/figures/offline.png)Figure 10:Experimental performance for varyingRRandδ\\delta, with sweep counts selected by our theory\. The left plot shows classification accuracy; the right plot shows total sweep cost\.While the experimental results in[Figure10](https://arxiv.org/html/2607.00170#S6.F10)overall confirm the theoretical prediction depicted in[Figure9](https://arxiv.org/html/2607.00170#S6.F9), there are some discrepancies to take note of\. First of all, we note that the experimentally measured good\-bad boundary forδ<0\.7\\delta<0\.7is rougher and less uniform than the theoretical prediction\. For instance, there is an unexpected bump in accuracy aroundδ=0\.4,R=0\.07\\delta=0\.4,R=0\.07\. We note, however, that the total sweep cost for this measurement is only 24, meaning an average of 4 sweeps per block\. By comparison with[Figure6](https://arxiv.org/html/2607.00170#S6.F6), we see that for small sweep counts, one more or less sweep can have a very significant effect on accuracy, especially at mid\-rangeδ\\delta\. Therefore, the relationship of estimated error to realized error becomes quantized at low costs, resulting in a significantly more jagged pattern: we may desire a certain stochastic error level, but we are forced to only do an integer number of sweeps\.

The second discrepancy between the two plots is more fundamental: the shape of the blue region beyondδ=δtrain=0\.7\\delta=\\delta\_\{\\text\{train\}\}=0\.7does not curve down smoothly as in the theoretical plot, but instead bumps outwards before curving back\. This is precisely the high\-δ\\delta‘overfitting’ effect that we discussed in[Section6\.1\.1](https://arxiv.org/html/2607.00170#S6.SS1.SSS1)\. The fact that it appears in the experiment, but not in the theory, is very telling regarding the etiology of the overfitting effect: the theoretical model does not consider spin\-spin correlation patterns, instead treating the statistical aggregate as composed of independent errors\. This assumption causes the overfitting effect to disappear\. We take this as evidence consistent with the overfitting effect being caused by spin\-spin correlation patterns and the non\-random interaction of thermal and stochastic error effects\. This makes sense as a hypothesis from a theoretical point of view, since the high\-δ\\deltaregime is precisely that mode of operation in which the spin\-spin couplings are strongest, and the independent\-spin mean field assumption that we used to construct the theory is least accurate\. However, on a practical note, we note that the best accuracy is achieved at lowδ\\delta, which is precisely where our theory is most effective\. In other words, our theory works where we need it to work\.

### 6\.5Regularizer Strength:λFP\\lambda\_\{\\mathrm\{FP\}\}Tradeoff

![Refer to caption](https://arxiv.org/html/2607.00170v1/figures/exp5_lambdafp_tradeoff.png)Figure 11:Log\-scaled plot showing error rates againstδ\\deltafor a network trained with varying values ofλFP\\lambda\_\{\\text\{FP\}\}\. As expected, lower values ofλFP\\lambda\_\{\\text\{FP\}\}lead to better performance at smallδ\\delta, but faster accuracy collapse asδ\\deltaincreases\. ForλFP\\lambda\_\{\\text\{FP\}\}too high, training collapses completely\.The parameterλFP\\lambda\_\{\\text\{FP\}\}controls how strong the Gibbs regularization is\. Since this regularization term is the most critical part of the training algorithm, the selection ofλFP\\lambda\_\{\\text\{FP\}\}is vital to achieving the desired results\. The key tradeoff is between STE accuracy and Gibbs robustness: smallλFP\\lambda\_\{\\text\{FP\}\}yields a higher accuracy ceiling, but demands higherGGand lowerδ\\deltato approximate it\. LargeλFP\\lambda\_\{\\text\{FP\}\}, on the other hand, produces a robust network which can operate at higherδ\\deltaand fewer sweeps, but has a lower maximum performance cap\.

We can view this tradeoff quite cleanly by doing an ablation experiment: we will train our 6\-block CIFAR\-10 model for several different values ofλFP\\lambda\_\{\\text\{FP\}\}, and plot its performance and cost when inference is performed with adaptive sweeping at a fixed value ofRR\. The results are plotted in[Figure11](https://arxiv.org/html/2607.00170#S6.F11)\. The experiment usedR=0\.01R=0\.01andδtrain=0\.7\\delta\_\{\\text\{train\}\}=0\.7\. We observe the expected behavior: lowerλFP\\lambda\_\{\\text\{FP\}\}leads to better maximum performance at smallδ\\delta, but also a quicker falloff in accuracy asδ\\deltaincreases\. LargerλFP\\lambda\_\{\\text\{FP\}\}leads to flatter performance, i\.e\. better resilience to largerδ\\delta, but also lower overall performance\. Accuracy always falls off beyondδ=δtrain\\delta=\\delta\_\{\\text\{train\}\}in any case, as expected\. The extreme values ofλFP\\lambda\_\{\\text\{FP\}\}are also worth noting:λFP=0\\lambda\_\{\\text\{FP\}\}=0\(regularizer disabled\) does produce great accuracy at very smallδ\\delta, leaning entirely on the magnitude regularizer, but this performance evaporates entirely with even slightly largerδ\\delta\. Extremely largeλFP\\lambda\_\{\\text\{FP\}\}, on the other hand, predictably causes total training collapse\.

## 7How Close Are We to Hardware?

This work is a hardware\-aware software paper\. Its purpose is to develop algorithms for an emerging class of probabilistic and thermodynamic devices, not to deploy them on currently available hardware\. Our experiments address a key gap in algorithmic research: can we scale the training of models capable of running on Ising hardware to practical sizes? Our results show that pure backpropagation training is capable of producing large\-scale deep models which function under Gibbs sampling and binary inter\-block communication, fulfilling key requirements for the development of software for Ising\-type systems\. If Ising hardware systems continue to become bigger and more powerful, our results indicate that the model\-training algorithms can scale alongside that hardware advancement\.

However, stability under Gibbs sampling and binary I/O are not the only requirements common to all existing Ising hardware\. In practice, almost every Ising hardware system also imposes some restrictions on weight quantization and sparsification, often requiring specific graph topologies\. Currently, a wide variety of experimental hardware options exist\. Some technologies promise dense graph topology, but extreme weight quantization, while some promise high weight precision at the cost of very restrictive graph topologies\. The approaches differ greatly in technology readiness level \(TRL\), roadmap, and forecasted advantages\. In short, a single “Ising hardware” standard simply does not exist: what does exist is a diverse family of prototypes, each with its unique caveats and tradeoffs\. We will first briefly review the landscape of currently available options, then discuss broadly how quantization and sparsification can be added to our models\.

### 7\.1Review of Current Hardware

At present, there is no single standard Ising\-machine architecture\. The phrase “Ising machine” covers several families of systems whose capabilities are not directly comparable\. Some are physical annealers, some are optical or spintronic stochastic systems, and others are quantum\-inspired digital or FPGA/GPU solvers\. Some operate on stochastic silicon, while some require more exotic materials or environments\. This diversity is encouraging, but it also means that “running on an Ising machine” is not a well\-defined hardware target\. The details of many machines are trade secrets, and the landscape is not always clear\. However, we will review some of the existing publicly\-known systems to understand the general sense of scale and the available tradeoffs\.

#### 7\.1\.1Superconducting Quantum Annealers

D\-Wave systems are the most mature commercial examples of physical Ising/QUBO\-oriented hardware\. The current Advantage2 system has thousands of qubits, tens of thousands of couplings, and a sparse Zephyr topology with bounded local connectivity\[dwave\_advantage2\_docs,dwave\_advantage2\_release\]\. These machines are highly relevant to the broader Ising\-computing ecosystem, but they are not the most relevant target for the current paper\. Quantum annealers are expensive and energy intensive, leaving little economic justification for running AI inference on them instead of a classical GPU\. Besides that constraint, they are also designed, as the name implies, forannealingto low energy states—a mode of operation incompatible with the high\-temperature sampling discussed in this paper\. On the other hand, a direct connection between quantum computing and AI models may indeed yet be useful, and the sampling\-temperature issue may not be unresolvable\.

#### 7\.1\.2Digital and CMOS Ising Systems

Digital annealing systems, including Fujitsu’s Digital Annealer and Hitachi’s CMOS/Momentum Annealing line, provide a different tradeoff\. They use conventional semiconductor hardware or hybrid software\-hardware systems to solve large QUBO or Ising optimization problems with high parallelism\[fujitsu\_da\_benchmark\_2025,hitachi\_relaxed\_ma\_2024\]\. Digital platforms can handle much larger logical problem sizes than sparse physical quantum annealers, and some support dense or effectively dense QUBO formulations\. This is attractive from the perspective of variable count and connectivity\. These platforms, as the name suggests, typically target annealing problems for combinatorial optimization, a typical traditional application of the Ising model for computation\. As digital systems, the underlying technology is more malleable, and likely could be modified to support our models’ sampling requirements\. However, the exact technical details of such systems are not fully known: these systems provide a useful gauge of scale for existing hardware, but specific discussions of compatibility are speculative\.

#### 7\.1\.3Simulated Bifurcation Machines

Toshiba’s Simulated Bifurcation Machine and SQBM\+ platform represent another important quantum\-inspired direction\. These systems implement Ising\-like optimization through highly parallelizable simulated bifurcation algorithms on GPUs, FPGAs, and cloud infrastructure, and current public materials describe support for 10 million spins, and up to 2,000 spins with fully connected interactions\[toshiba\_sqmplus\_2026,toshiba\_sbm\_edgechaos\_2026\]\. This makes the platform significant for large\-scale Ising optimization\. From the perspective of thermodynamic neural inference, the same caveat applies: simulated bifurcation is a solver for optimization dynamics, not a drop\-in replacement for the fixed\-temperature Gibbs time averages used in our blocks\. Nonetheless, these machines are large\-scale implementations of the Ising system, and are implemented on malleable substrates like FPGAs\. In other words, this technology likely can support our models, but needs modifications to its harness to do so\.

#### 7\.1\.4Coherent Ising Machines

Coherent Ising machines use optical degrees of freedom, such as degenerate optical parametric oscillator pulses, to represent spins\. Large measurement feedback systems have demonstrated very large spin counts and dense effective couplings, including100,512100\{,\}512spins and more than1010B pairwise interactions in an NTT coherent Ising machine\[ntt\_cim100k\]\. This scale is impressive and relevant to the long\-term hardware story; in fact, dense connectivity at these scales is far more than we need\. However, current research focuses on sampling a single amplified spin configuration—in other words, Boltzmann time averaging would require a significant modification\. Additionally, weight quantization requirements are extreme; coupling strengths on published machines are required to be in\{−1,0,1\}\\\{\-1,0,1\\\}\. Engineering a model to fit into this ternary quantization is not impossible, but would require tradeoffs and aggressive quantization\-aware training algorithms\.

#### 7\.1\.5Probabilistic\-Bit and Spintronic Ising Machines

The closest conceptual match to our model is the probabilistic\-bit, or p\-bit, family of machines\. These systems implement tunable stochastic binary units and update them using coupling logic, making them naturally aligned with Gibbs\-like sampling and Boltzmann\-machine\-style computation\. Recent magnetic tunnel junction demonstrations show device\-level progress, but current experimental systems remain much smaller than the networks studied here\[yang2026mtj\_pim\]\. For this reason, p\-bit and related spintronic probabilistic architectures may be the most natural long\-term substrate for thermodynamic neural inference, even though they are not yet at the required scale\.

#### 7\.1\.6FPGA and GPU Gibbs Samplers

Finally, classical FPGA and GPU implementations of Gibbs sampling are important baselines\. They are not thermodynamic physical devices in the same sense as spintronic or optical systems, but they can implement the relevant Markov chain directly and can therefore directly test the algorithmic interface assumed in this paper\. Prior Gibbs\-sampled Ising classifiers such as\[Niazi\_2024\]are especially useful in this respect\. They show that sweep\-based Ising inference can be accelerated substantially by dedicated hardware, bringing a large number of sweeps within a reasonable cost range both in terms of energy expenditure and inference time\. We discuss this comparison in more detail in[Section7\.3](https://arxiv.org/html/2607.00170#S7.SS3)\.

### 7\.2Sparsification and Quantization

Practical Ising hardware and finite\-precision Ising emulators typically require quantized couplings\. Many also require strong sparsity or structured topology\. As we have discussed, there is no standard for exactly how much quantization or sparsification is required\. Therefore, we have no exact target to shoot for\. However, since both ingredients are a key part of making thermodynamic AI systems practical on Ising machines, we will discuss in some generality how sparsification and quantization can factor into our algorithms\.

Although convolutional layers are sparse compared with fully connected layers, large channel counts still produce high local degree after unrolling\. In the last block of our CIFAR\-10 model, for example, the hidden layer has784784channels\. A hidden spin connected by3×33\\times 3filters to both the input and output layers has local degree

2⋅32⋅784=14,112,2\\cdot 3^\{2\}\\cdot 784=14\{,\}112,which is the largest local degree in the model\. Reducing this to maximum degree141141would require roughly99%99\\%sparsity, while reducing it to maximum degree1414would require roughly99\.9%99\.9\\%sparsity\.

These are aggressive targets, beyond what should be expected from a purely post\-training pruning pass\. Such low spin degrees would entail integrating training techniques such as quantization\-aware training, sparsity\-aware training, structured channel sparsity, grouped convolutions, low\-rank factorizations, and lottery\-ticket\-style dynamic sparsification\[evci2021rigginglotterymakingtickets\]\. These options are all natural extensions of the present pipeline\. In the current paper, however, we intentionally keep the main algorithm simple and use a post\-training compression experiment to demonstrate that our trained thermodynamic models are, in fact, amenable to quantization and sparsification, and not excessively brittle\.

We apply Wanda\[sun2024simpleeffectivepruningapproach\]to the moderate\-λFP\\lambda\_\{\\mathrm\{FP\}\}six\-block CIFAR\-10 model\. Before compression, this model has feed\-forward STE accuracy92\.75%92\.75\\%and Gibbs inference accuracy92\.13%92\.13\\%atδ=0\.1\\delta=0\.1\. With77\-bit weight quantization and80%80\\%sparsity, the feed\-forward STE accuracy falls to88\.59%88\.59\\%, while Gibbs inference accuracy atδ=0\.1\\delta=0\.1falls to86\.7%86\.7\\%\. The maximum per\-spin degree is reduced to approximately2\.82\.8k\. This suggests that while simple post\-training approaches have their limitations, our models have some resilience to compression\. Our Gibbs\-stability tuning methods appear reasonably robust and do not decay much faster than overall feed\-forward accuracy does under compressive pressure\. While this result is not yet aggressive enough to fit our CIFAR\-10 model onto current hardware, it is a good sign that more advanced methods are worth implementing, and that there is a clear roadmap for the training of compressed hardware\-compatible models using our methods\.

### 7\.3Sweep Count

Inference cost needs to be conceptually separated from hardware capacity\. Capacity asks whether a device can represent the required spins, couplings, precision, and I/O pattern\. Inference cost, however, refers to how much must be expended to obtain a result from running that device\. In our model, inference cost is measured in terms of Gibbs sweeps, though hardware implementing other dynamics would have some other analogous measure of system simulation cost\. Our convergence theory was motivated by the need to control inference cost\. A model which works in theory but requires impractically high sweep count is not very useful\. Therefore, we have quantified our sweep cost and shown theoretically that it does not diverge unexpectedly\. However, the reader may still wonder how realistic the sweep counts listed in[Table1](https://arxiv.org/html/2607.00170#S4.T1)actually are in the current hardware context\.

For comparison,\[Niazi\_2024\]uses10510^\{5\}sweeps for MNIST and FashionMNIST\-scale Gibbs\-sampled Ising classification on an FPGA device, claiming that this number of sweeps can be executed within a few milliseconds and at low power consumption\. Our smallest models use three orders of magnitude fewer sweeps, while our larger models still use substantially fewer sweeps than the previous result\. This does not mean that any current device can run our large models within a few milliseconds, but it does show that sweep counts of order10510^\{5\}can be practical on at least one specialized implementation\. Whether that timing transfers to our much larger models remains a hardware question, and power expenditure should be expected to depend on model size and architecture\. In short, this comparison shows that sweep count does not by itself rule out practical implementation, and may be a point of improvement over existing results\. However, without a fair and direct comparison, we cannot make any concrete and precise conclusions about the time or energy cost of running our models in real hardware\.

### 7\.4Implications for Hardware Design

The goal of thermodynamic AI research is to create useful models that can run on emerging probabilistic hardware\. We have developed algorithms that can scale alongside Ising hardware as its capabilities expand and provided evidence that training algorithms need not be the immediate limiting factor for implementing AI models on thermodynamic systems\. However, our algorithmic work is also important because it provides a target for hardware to shoot for\. Thermodynamic AI will not become a reality at scale unless hardware research considers it a possibility and works towards moving Ising hardware in a direction that can support it\. A common theme in our review of existing hardware is that most approaches have targeted the traditional application of Ising systems, annealing combinatorial optimization problems\. This can lead to the development of more powerful Ising systems, but thermodynamic AI will require different adaptations of the underlying system to become a reality\. To this end, we will elaborate the implications of our algorithmic research for the development of Ising hardware\.

The first implication is that future thermodynamic inference hardware should be developed to support convolutional graph topology\. Convolutional inference, being naturally sparse, is a natural direction for Ising hardware\. It also presents many promising applications for low\-power edge inference, from implantable medical sensors to miniature surveillance systems, due to its natural applications in image, sensor, and audio processing\. However, most existing hardware was not designed with this workload in mind\. Hardware approaches which require specific graph topologies could be designed using an unrolled\-convolutional network topology\. On the other hand, digital systems could exploit weight sharing and translational symmetry to compress couplings in convolutional Ising layers in the same way that GPUs do\. Since convolutional networks are a natural application of the technology, and training algorithms for programming convolutional Ising hardware now exist, it makes sense for convolutional weights to become a hardware\-level priority\.

The second implication is that average\-and\-sign readout should be implemented as a hardware\-native operation\. Much Ising\-machine work focuses on either ground\-state readout or instantaneous Boltzmann samples\. Our computation uses a different observable: the sign of the time average of an output spin\. This can be computed in host software by collecting repeated Boltzmann samples, but doing so wastes bandwidth and moves a simple reduction away from the device\. A thermodynamic neural chip should attach lightweight accumulation circuits to designated output spins, track their running averages over the sweep window, and export only the final signs or low\-precision average estimates\.

The third implication is that input spin fixing should be a first\-class hardware feature\. Throughout the paper, each thermodynamic block receives a binary input tensor and treats the corresponding input spins as fixed while the hidden and output spins evolve\. Existing platforms can sometimes approximate this by applying large local fields, but this is not the same as a clean hardware\-level clamping mechanism\. Our neural inference naturally uses Boltzmann distributions conditional on the input spins, so a device intended for this workload should support input pinning and conditional Ising simulation\. On many platforms, especially digital implementations, such a change could be straightforward, but it has to be kept in mind during the design process\.

## 8Conclusion

In this paper, we have exhibited a backpropagation\-based training algorithm for creating large convolutional networks which perform under block\-wise Ising\-system Gibbs sampling\. The ability to train such a network at scale is a necessary building block for implementing AI models on Ising\-type thermodynamic hardware\. We have furthermore shown experiments and calculations demonstrating the relatively fast convergence and high performance of these models, as well as their ability to handle sizable benchmarks like CIFAR\-100 with millions of spins\. While not yet at the level of being uploaded directly to hardware, we have de\-risked the production of Ising\-compatible models and shifted the natural next focus towards enforcing device\-specific hardware compatibility constraints\.

In other words, deep thermodynamic neural inference is now algorithmically viable enough to become a hardware\-software codesign target\. With the right training procedure and hardware implementation, Ising\-compatible stochastic dynamics can support modern convolutional inference, and the cost of that inference can be analyzed, predicted, and designed rather than guessed\. This paper takes the step from a theoretical correspondence between neural inference and stochastic systems toward a practical framework for scaling thermodynamic AI\. That said, in this paper we have focused on arguing the core question of viability more than elaborating on every possible direction or pushing for detailed hardware implementation\. Therefore, there are many natural directions open for future research that aims to flesh out, exploit, and realize the above results\.

Implementing a thermodynamic computer in real hardware requires two things\. The software must be designed to fit the hardware specifications, and the hardware must support the demands of the software\. To see the models we have discussed running on real Ising hardware, both sides must be addressed\. First, a specific hardware platform must be picked as a codesign target\. Then, the training algorithms must be refined to target that platform and its particular specifications\. In reality, this would likely look like a modification to the existing procedure which integrates the right flavor of quantization, sparsification, and graph topology restriction\. In short, specialization and algorithmic refinement are needed for implementing our models on any particular hardware platform\.

Not all thermodynamic computing hardware is based on the Ising model\. Some systems operate according to completely different principles, but another large class consists of hardware devices that implement a physical model closely related to but not identical to the Ising model\. A common example is the Kuramoto model for coupled oscillators, implemented for example by coupled LC oscillator arrays with sub\-harmonic injection locking\[11186164,shil\]\. In the low\-temperature \(i\.e\. strong locking\) limit, the Kuramoto Hamiltonian has the same ground states as the analogous Ising Hamiltonian\. However, more weakly coupled Kuramoto systems, possibly with injected stochastic perturbations or possibly with modified oscillators\[khan2025analyzingparametricoscillatorising\], approximate the Ising model inexactly\. Nonetheless, the existence of prototype Kuramoto hardware indicates that it may be worth extending our mean\-field theoretical correspondence\. If we can replicate our methods for training thermodynamic hardware with cousins of the Ising model, the number of possible hardware application targets will increase\.

In the paper, we have focused only on using our convolutional networks for image classification\. This was mostly for simplicity and uniformity\. In fact, our methods should extend to a broad class of convolutional networks\. One\-dimensional convolutional networks could do classification problems for audio or biomedical sensor data\. On the other hand, mostly\-convolutional models like object detection and signal detection are a natural fit for extending the sophistication of our models without greatly changing the training algorithm\. Really, many models that are mostly convolutional networks could be adapted to a thermodynamic system\. Image segmentation and diffusion\-based image generation are also possible directions\. However, more sophisticated application of thermodynamic models may require more sophisticated host data processing\. The fact that only binary data is passed between our blocks seriously interferes with the performance of a diffusion model, for instance\. On the other hand, adding lightweight host processing or residual data streams between the primary thermodynamic blocks could open up a whole family of powerful thermodynamic models trained with simple variations or extensions of the framework presented here\.

## 9Acknowledgments

A\. G\. M\. would like to thank Professor Thomas Chen of the Department of Mathematics at the University of Texas at Austin for his supervision and guidance\. He also thanks Professor Adam Klivans and the National Science Foundation through the Institute for Foundations of Machine Learning \(IFML\) for their support\.

## References

## Appendix APretraining

The purpose of this sections is to provide supplemental algorithmic detail explaining how models are prepared for Gibbs regularization \([Section3](https://arxiv.org/html/2607.00170#S3)\)\. In short, we train a standard model with bounded activations \(Phase 1\), then convert it to the right shape \(Phase 2\), then train it as an STE across block boundaries \(Phase 3\)\. This gets the system close enough to an Ising\-compatible solution that our main training algorithm becomes effective\. We refer to this process as ‘pretraining’ because it produces a trained model as input for our main training algorithm\.

We begin with Phase 1, training a conventional residual convolutional model, the WideResNet, whose block structure is conveniently thermodynamically friendly\. Every block contains an internal activation after the first convolution and an output activation after the residual sum\. Both are implemented with the penalized tanh activation\[xu2016revisesaturatedactivationfunctions\]

ptanh​\(x\)=\{tanh⁡\(x\),x≥0,0\.25​tanh⁡\(x\),x<0\.\\mathrm\{ptanh\}\(x\)=\\begin\{cases\}\\tanh\(x\),&x\\geq 0,\\\\ 0\.25\\,\\tanh\(x\),&x<0\.\\end\{cases\}The sharp corner in the activation function provides a ReLU\-like decision boundary which helps the network learn, while keeping overall activations bounded in a roughly sigmoid shape\.

The role of this phase is simply to learn a strong bounded\-activation teacher model before any binarization or Gibbs\-stability constraints are imposed\. Because every inter\-block activation is already bounded, the model starts in a regime that is structurally close to the binary communication pattern required by thermodynamic hardware, while still being easy to optimize with ordinary deep\-learning tools\.

In the implementation we train this model with standard large\-scale image\-classification practice: SGD with Nesterov momentum, cosine decay with warmup, label smoothing, dropout in the hidden path, and dataset\-appropriate augmentation including random crop, horizontal flips where appropriate, AutoAugment, and Cutout\.

### A\.1Adaptation & Folding

Phase 2 of pretraining converts the bounded continuous teacher into an exact tanh\-plus\-affine parameterization that is compatible with later binarization and Gibbs interpretation\. We first replace every penalized tanh by a fixed affine\-transformed tanh approximation

f​\(x\)=a​tanh⁡\(b​x\+c\)\+df\(x\)=a\\tanh\(bx\+c\)\+dwith

a=0\.625,b=1\.239,c=−0\.523,d=0\.375\.a=0\.625,\\qquad b=1\.239,\\qquad c=\-0\.523,\\qquad d=0\.375\.These constants were chosen to minimizeL2L^\{2\}loss against ptanh\. We continue the training process until we recover the accuracy lost by the activation function replacement\.

Finally, we perform an exact algebraic fold of the inner affine parameters\(b,c\)\(b,c\)into the preceding BatchNorm layers\. After this fold, each activation function has the form

u↦a⊙tanh⁡\(u\)\+d\.u\\mapsto a\\odot\\tanh\(u\)\+d\.This fold is exact; it is not an approximation\. We then move theaaandddparameters of the output\-layer activations to the input layer of the next block\. This results in a network whose inter\-block data is in the range\[−1,1\]\[\-1,1\]and is of the formtanh⁡\(⋯\)\\tanh\(\\cdots\)\. This network is now ready for binarization: by replacing the final tanh with sign, we produce a network based exclusively on affine transforms and tanh which passes only binary tensors with values in\{±1\}\\\{\\pm 1\\\}between blocks\. This is exactly the sort of raw material that we need for our Ising hardware implementation\.

### A\.2Straight\-Through Estimator

Once the continuous pretraining is complete, we switch the architecture to a straight\-through estimator \(STE\) and continue to train\. The goal is to adapt the folded tanh network to binary inter\-block communication\. The binarization only happens between blocks: within a block, we still use tanh, because the block is conceptualized as a high\-temperature mean\-field Ising system\. The block structure produced by phase 2 \(after batch\-norm folding\) is

y​\(x\)=tanh⁡\(K3​Ain​x\+K2​Amid​tanh⁡\(b\+K1​Ain​x\)\)\\displaystyle y\(x\)=\\tanh\(K\_\{3\}A\_\{\\text\{in\}\}x\+K\_\{2\}A\_\{\\text\{mid\}\}\\tanh\(b\+K\_\{1\}A\_\{\\text\{in\}\}x\)\)\(46\)We now replace the final activation function withζ\\zeta, the surrogate sign function:

y​\(x\)=ζ​\(K3​Ain​x\+K2​Amid​tanh⁡\(b\+K1​Ain​x\)\)\\displaystyle y\(x\)=\\zeta\(K\_\{3\}A\_\{\\text\{in\}\}x\+K\_\{2\}A\_\{\\text\{mid\}\}\\tanh\(b\+K\_\{1\}A\_\{\\text\{in\}\}x\)\)\(47\)During training, the boundary activation functionζ\\zetais a straight\-through Bernoulli estimator ofsign\\operatorname\{sign\}\. The forward pass uses Bernoulli spins drawn from the corresponding local fields, while the backward pass uses the derivative oftanh\\tanhas a smooth surrogate\. We use88stochastic Bernoulli forward samples per step, adding noise tolerance and reducing gradient variance\. During inference, on the other hand,ζ\\zetais set to the hardsign\\operatorname\{sign\}function\. In other words,

Inference:​ζ​\(z\)=sign⁡\(z\),\\displaystyle\\text\{Inference: \}\\zeta\(z\)=\\operatorname\{sign\}\(z\),Training:​ζ​\(z\)=ℬ​\(σ​\(2​z\)\),\\displaystyle\\text\{Training: \}\\zeta\(z\)=\\mathcal\{B\}\(\\sigma\(2z\)\),d​ζd​z=1−tanh2⁡\(z\)\\displaystyle\\frac\{d\\zeta\}\{dz\}=1\-\\tanh^\{2\}\(z\)\(48\)As usual for straight\-through estimators, the forward and backward definitions ofζ\\zetaare intentionally different: the forward pass uses a discrete stochastic estimator, while the backward pass uses the derivative of a smooth surrogate\. Note that all inter\-block data is now in\{±1\}\\\{\\pm 1\\\}, a critical part of our final model criteria\. The internal activations are still continuous, and based on tanh\. This is not an issue because the intra\-block activations now represent mean\-field averages of the fully binary Gibbs process\.

Rather than switching abruptly from continuous to binary behavior, the script uses a short temperature\-anneal period in which the boundary activation interpolates between tanh and Bernoulli; after this warmup, training proceeds fully in Bernoulli/sign mode\. To reduce variance, the training step averages multiple stochastic forward samples of the same minibatch\.

The STE training also uses two teachers\. The first is the original ptanh network, which provides a standard logit\-level distillation target\. The second is a frozen copy of the tanh\-adapted network, which provides feature\-level distillation on the outputs of the three residual stages\. In addition, we regularize the per\-channel offsets toward zero in order to encourage centered binary activations, and we penalize excessively large BatchNorm scales so that the subsequent folding step remains well behaved\. We train this STE as far as it will go before we start adding Gibbs regularization\.

Finally, we eliminate elements of the network which are not affine or tanh by a folding process\. The parametersbbandccare algebraically absorbed into the preceding batch normalization layer by the transformationBN′=b⋅BN\+c\\text\{BN\}^\{\\prime\}=b\\cdot\\text\{BN\}\+c, and the batch normalization is then folded into the convolution weights and biases\. After folding, each activation site has the forma⊙tanh⁡\(z\)\+da\\odot\\tanh\(z\)\+d, wherezzis the raw convolution output\.

## Appendix BAlgorithm Variations for Gibbs Regularization

### B\.1Block\-Weighted Magnitude Regularization

It can sometimes help to adjust the target magnitudes of the magnitude regularizer according to block position\. At blockbbofnbn\_\{b\}total blocks, we can scale the target magnitude by

1\+15​\(1−bnb\)\.\\displaystyle 1\+\\frac\{1\}\{5\}\\left\(1\-\\frac\{b\}\{n\_\{b\}\}\\right\)\.\(49\)This essentially tries to make earlier blocks more confident, which prevents the pile\-up of stochastic error across blocks that can cause problems in deeper models\. It should be noted that because the magnitude regularizer roughly determines the position of the main spike in theρ\\rhodistribution, applying a block\-weighted magnitude regularizer in this way will extend the period of fast convergence for earlier blocks\.

### B\.2Flip\-Rate\-Weighted Fixed\-Point Loss

In multi\-block models, the per\-block fixed\-point lossesℒFP,i\\mathcal\{L\}\_\{\\text\{FP\},i\}of[Section3\.1](https://arxiv.org/html/2607.00170#S3.SS1)are aggregated to form a single training signal\. The default recipe used in our experiments simply sums them,

ℒFPtotal=∑i=1nbℒFP,i,\\mathcal\{L\}\_\{\\text\{FP\}\}^\{\\text\{total\}\}\\;=\\;\\sum\_\{i=1\}^\{n\_\{b\}\}\\mathcal\{L\}\_\{\\text\{FP\},i\},\(50\)treating every block as equally important at every training step\. A natural alternative is to allocate gradient mass to each block in proportion to how unstable that block currently is, so that whichever block currently lags behind in Gibbs stability receives the largest share of each update\. We describe this variant here as a drop\-in replacement for \([50](https://arxiv.org/html/2607.00170#A2.E50)\); on shallow models the two aggregations behave essentially identically in practice, but on deeper or wider networks one or two blocks can dominate the residual instability and benefit from preferential updates\. We define the per\-block flip ratefi∈\[0,1\]f\_\{i\}\\in\[0,1\]as the fraction of spin sites at blockiiat which the Gibbs consensusygibbs\(i\)y\_\{\\text\{gibbs\}\}^\{\(i\)\}disagrees with the feed\-forward STE sign of the block\-output fieldz2\(i\)z\_\{2\}^\{\(i\)\}:

fi=⟨𝟏​\[sign⁡z2\(i\)≠ygibbs\(i\)\]⟩,f\_\{i\}\\;=\\;\\left\\langle\\bm\{1\}\\\!\\left\[\\,\\operatorname\{sign\}z\_\{2\}^\{\(i\)\}\\,\\neq\\,y\_\{\\text\{gibbs\}\}^\{\(i\)\}\\,\\right\]\\right\\rangle,\(51\)where⟨⋅⟩\\langle\\cdot\\rangledenotes the empirical mean over the batch and over the spatial–channel positions of blockii\. The flip rate is read off the Gibbs chain that already runs to produce the fixed\-point target of[Section3\.1](https://arxiv.org/html/2607.00170#S3.SS1), and is treated as a scalar with no gradient\. It is large when the block’s feed\-forward sign and its Gibbs equilibrium disagree, and shrinks toward zero as the fixed\-point loss drives them into agreement\. Given\{fi\}i=1nb\\\{f\_\{i\}\\\}\_\{i=1\}^\{n\_\{b\}\}, form unnormalised weights with a fixed baselineα\>0\\alpha\>0,

w~i=α\+fi,\\tilde\{w\}\_\{i\}\\;=\\;\\alpha\+f\_\{i\},\(52\)and renormalise so the weights sum tonbn\_\{b\}, preserving the overall scale of the unweighted aggregate \([50](https://arxiv.org/html/2607.00170#A2.E50)\):

wi=nb​w~i∑j=1nbw~j=nb​\(α\+fi\)nb​α\+∑j=1nbfj\.w\_\{i\}\\;=\\;\\frac\{n\_\{b\}\\,\\tilde\{w\}\_\{i\}\}\{\\sum\_\{j=1\}^\{n\_\{b\}\}\\tilde\{w\}\_\{j\}\}\\;=\\;\\frac\{n\_\{b\}\\,\(\\alpha\+f\_\{i\}\)\}\{n\_\{b\}\\,\\alpha\+\\sum\_\{j=1\}^\{n\_\{b\}\}f\_\{j\}\}\.\(53\)The flip\-rate\-weighted aggregate that replaces \([50](https://arxiv.org/html/2607.00170#A2.E50)\) is then

ℒFPtotal=∑i=1nbwi​ℒFP,i\.\\mathcal\{L\}\_\{\\text\{FP\}\}^\{\\text\{total\}\}\\;=\\;\\sum\_\{i=1\}^\{n\_\{b\}\}w\_\{i\}\\,\\mathcal\{L\}\_\{\\text\{FP\},i\}\.\(54\)We useα=0\.2\\alpha=0\.2in all of our experiments with this variant\. The construction has three convenient properties\. First, when all blocks have the same flip rate,fi≡cf\_\{i\}\\equiv c, the weights collapse towi=1w\_\{i\}=1and \([54](https://arxiv.org/html/2607.00170#A2.E54)\) reduces exactly to the unweighted aggregate \([50](https://arxiv.org/html/2607.00170#A2.E50)\); only*relative*differences in flip rate redistribute mass between blocks\. Second, the baselineα\\alphalower\-bounds each block’s share: asα→0\\alpha\\to 0, a fully\-stable block \(fi=0f\_\{i\}=0\) would receivewi=0w\_\{i\}=0and stop receiving fixed\-point gradient at all, while asα→∞\\alpha\\to\\inftythe weights flatten back to uniform and we recover \([50](https://arxiv.org/html/2607.00170#A2.E50)\)\. Third, the maximum amplification any single block can receive is bounded:

maxi⁡wi≤nb​\(1\+α\)1\+α\+\(nb−1\)​α,\\max\_\{i\}w\_\{i\}\\;\\leq\\;\\frac\{n\_\{b\}\\,\(1\+\\alpha\)\}\{1\+\\alpha\+\(n\_\{b\}\-1\)\\,\\alpha\},\(55\)attained when the worst block hasfk=1f\_\{k\}=1and the others havef=0f=0\. Fornb=3n\_\{b\}=3andα=0\.2\\alpha=0\.2this caps amplification atwk≤2\.25w\_\{k\}\\leq 2\.25, with the corresponding stable blocks receivingw≈0\.375w\\approx 0\.375\. Mechanistically,wiw\_\{i\}acts as a per\-block, per\-step learning\-rate multiplier that concentrates gradient mass on whichever block currently sits at the bottleneck of Gibbs stability\. Once the fixed\-point loss drives that block’s flip rate down, the weighting redistributes to the next bottleneck\. The unweighted aggregate \([50](https://arxiv.org/html/2607.00170#A2.E50)\) is the limitα→∞\\alpha\\to\\inftyof the same scheme, and the two aggregations agree exactly whenever the per\-block flip rates are equal, which is empirically the regime our shallowest models operate in for most of Gibbs regularization training\.

## Appendix CIntegrated Autocorrelation Time

The missing piece in our convergence analysis is an estimate of the integrated autocorrelation time for each output spin\. Numerical experiments show that in our convolutional Ising blocks, the measured per\-spin integrated autocorrelation time is typically within\[1,1\.5\]\[1,1\.5\]—essentially, our chains are nearly i\.i\.d\. We would like to justify, understand, or predict this theoretically\.

### C\.1Dobrushin Approximation

To demonstrate that we can at least control the autocorrelation time, it is sufficient to find an upper bound on the maximum integrated autocorrelation timeτint∗\\tau\_\{\\text\{int\}\}^\{\*\}\. We will start with the easy case: Glauber \(not red\-black Gibbs\) dynamics andδ→0\\delta\\rightarrow 0\. Consider the natural approximation of the Dobrushin constant

α~​\(δ\):=maxi​∑j≠itanh⁡\(δ​\|Ji​j\|\)\.\\displaystyle\\tilde\{\\alpha\}\(\\delta\):=\\max\_\{i\}\\sum\_\{j\\neq i\}\\tanh\(\\delta\|J\_\{ij\}\|\)\.\(56\)Usingtanh⁡x≤x\\tanh x\\leq xforx≥0x\\geq 0, define

κ:=maxi​∑j≠i\|Ji​j\|\.\\displaystyle\\kappa:=\\max\_\{i\}\\sum\_\{j\\neq i\}\|J\_\{ij\}\|\.\(57\)Thenα~​\(δ\)≤δ​κ\\tilde\{\\alpha\}\(\\delta\)\\leq\\delta\\kappa, andα~​\(δ\)=δ​κ\+𝒪​\(δ3\)\\tilde\{\\alpha\}\(\\delta\)=\\delta\\kappa\+\\mathcal\{O\}\(\\delta^\{3\}\)asδ→0\\delta\\rightarrow 0when the maximizing row is stable\. A Dobrushin\-type argument gives, forδ​κ<1\\delta\\kappa<1, a bound of the form

τint∗≤C01−α~​\(δ\)≤C01−δ​κ,\\displaystyle\\tau\_\{\\text\{int\}\}^\{\*\}\\leq\\frac\{C\_\{0\}\}\{1\-\\tilde\{\\alpha\}\(\\delta\)\}\\leq\\frac\{C\_\{0\}\}\{1\-\\delta\\kappa\},\(58\)whereC0C\_\{0\}is independent ofδ\\deltain this regime\. Combining this estimate with the leading stochastic\-error scaling gives

R​\(G\)≲G−1/2δ​1−δ​κ\.\\displaystyle R\(G\)\\lesssim\\frac\{G^\{\-1/2\}\}\{\\delta\\sqrt\{1\-\\delta\\kappa\}\}\.\(59\)Within this loose bound, we can optimize the choice ofδ\\delta:

−1δ\+κ2​\(1−δ​κ\)=0⟹δ=23​κ\.\\displaystyle\-\\frac\{1\}\{\\delta\}\+\\frac\{\\kappa\}\{2\(1\-\\delta\\kappa\)\}=0\\quad\\Longrightarrow\\quad\\delta=\\frac\{2\}\{3\\kappa\}\.\(60\)
This analysis has a major and unavoidable limitation:δ\\deltacontrols not only the convergence speed, but also the correctness of the final distributionμ\\mu\. For finitely many nonzero outputs, we know that there exists someΔ\\Deltasuch that0<δ<Δ0<\\delta<\\Deltaimplies that the signs ofμ\\muagree with the feed\-forward signs\. However, in practice, increasingδ\\deltacomes with an accuracy tradeoff because the stable distribution drifts away from the low\-δ\\delta‘correct’ distribution\. In other words, the optimalδ\\deltafor convergence rate to the stable distribution may not be the optimalδ\\deltafor achieving accurate inference\.

However, in the context of the above estimate, this caveat is not actually very relevant: the fact is that the Dobrushin bound is simply far too loose to be useful, and theδ\\deltait predicts is so small that accuracy distortions aren’t relevant in that regime\. In the case of our CIFAR\-10 network, the measuredα~​\(1\)\\tilde\{\\alpha\}\(1\)values of the blocks are in the region of 50 to 250, and we observe that optimal convergence speeds are occurring forδ\\deltalarger than1/101/10\.

### C\.2Linearized Red\-Black Mode Theory

The Dobrushin analysis shows that autocorrelation does not explode asδ→0\\delta\\rightarrow 0\. However, it is not a particularly useful bound, and only applies toδ\\deltasmall enough that it’s no longer practically relevant\. Now, we develop a more refined theory which is capable of giving us much tighter estimates, at the cost of only being an approximation\. We will also present numerical experiments showing that, in practice, it is accurate and slightly conservative\.

The approach is to linearize the red\-black Gibbs sampling, giving us a linearized update rule for each layer\. We expand around the Jacobian of naive mean\-field theory\. This gives us a symmetric positive semidefinite operator\. We then use eigendecomposition to determine the dominant modes, i\.e\. those with the largest eigenvalues\. The resulting eigenvalues directly control the autocorrelation times of those modes\. Our goal is to use this mode decomposition to obtain a practical upper bound on the per\-spin integrated autocorrelation time\. Since each coordinate observable is a convex combination of the modal observables, the top eigenvalue yields an immediate upper bound on all per\-spin integrated autocorrelation times\.

Consider a single thermodynamic block with two spin layers, denoteds1s\_\{1\}ands2s\_\{2\}, updated by alternating red\-black Gibbs sweeps\. LetKKbe the convolution coupling between the two layers, and

A:=diag⁡\(amid\)\.A:=\\operatorname\{diag\}\(a\_\{\\mathrm\{mid\}\}\)\.The conditional mean maps have the form

𝔼​\[s1\+∣s2\]\\displaystyle\\mathbb\{E\}\[s\_\{1\}^\{\+\}\\mid s\_\{2\}\]=tanh⁡\(b1\+δ​A​K⊤​s2\),\\displaystyle=\\tanh\\\!\\big\(b\_\{1\}\+\\,\\delta\\,AK^\{\\top\}s\_\{2\}\\big\),\(61\)𝔼​\[s2\+∣s1\]\\displaystyle\\mathbb\{E\}\[s\_\{2\}^\{\+\}\\mid s\_\{1\}\]=tanh⁡\(δ​\(K​\(A​s1\+dmid\)\+b2\)\),\\displaystyle=\\tanh\\\!\\big\(\\delta\\,\(K\(As\_\{1\}\+d\_\{\\mathrm\{mid\}\}\)\+b\_\{2\}\)\\big\),\(62\)whereb1b\_\{1\}andb2b\_\{2\}collect the image\-dependent biases given by the coupling with the input layer and the fixed bias, andδ\\deltais the post\-training temperature scale\.

Letm1,m2m\_\{1\},m\_\{2\}denote the stationary means of the two layers, and define the diagonal susceptibility matrices

D1:=diag⁡\(1−m12\),D2:=diag⁡\(1−m22\)\.\\displaystyle D\_\{1\}:=\\operatorname\{diag\}\(1\-m\_\{1\}^\{2\}\),\\qquad D\_\{2\}:=\\operatorname\{diag\}\(1\-m\_\{2\}^\{2\}\)\.\(63\)These are the coordinatewise derivatives of the conditional mean maps at stationarity\. We then introduce normalized centered fluctuations

z1:=D1−1/2​\(s1−m1\),z2:=D2−1/2​\(s2−m2\)\.\\displaystyle z\_\{1\}:=D\_\{1\}^\{\-1/2\}\(s\_\{1\}\-m\_\{1\}\),\\qquad z\_\{2\}:=D\_\{2\}^\{\-1/2\}\(s\_\{2\}\-m\_\{2\}\)\.\(64\)Symmetrically linearizing the two half\-steps around stationarity gives

z1\+\\displaystyle z\_\{1\}^\{\+\}≈H1←2​z2,\\displaystyle\\approx H\_\{1\\leftarrow 2\}z\_\{2\},H1←2\\displaystyle H\_\{1\\leftarrow 2\}=δ​D11/2​A​K⊤​D21/2,\\displaystyle=\\delta\\,D\_\{1\}^\{1/2\}AK^\{\\top\}D\_\{2\}^\{1/2\},\(65\)z2\+\\displaystyle z\_\{2\}^\{\+\}≈H2←1​z1,\\displaystyle\\approx H\_\{2\\leftarrow 1\}z\_\{1\},H2←1\\displaystyle H\_\{2\\leftarrow 1\}=δ​D21/2​K​A​D11/2\.\\displaystyle=\\delta\\,D\_\{2\}^\{1/2\}KAD\_\{1\}^\{1/2\}\.\(66\)Composing the two half\-steps yields the same\-layer full\-sweep linearizations

z1\(t\+1\)\\displaystyle z\_\{1\}^\{\(t\+1\)\}≈S1​z1\(t\),\\displaystyle\\approx S\_\{1\}z\_\{1\}^\{\(t\)\},S1\\displaystyle S\_\{1\}:=H1←2​H2←1,\\displaystyle:=H\_\{1\\leftarrow 2\}H\_\{2\\leftarrow 1\},\(67\)z2\(t\+1\)\\displaystyle z\_\{2\}^\{\(t\+1\)\}≈S2​z2\(t\),\\displaystyle\\approx S\_\{2\}z\_\{2\}^\{\(t\)\},S2\\displaystyle S\_\{2\}:=H2←1​H1←2\.\\displaystyle:=H\_\{2\\leftarrow 1\}H\_\{1\\leftarrow 2\}\.\(68\)Explicitly,

S1\\displaystyle S\_\{1\}=δ2​D11/2​A​K⊤​D2​K​A​D11/2,\\displaystyle=\\,\\delta^\{2\}\\,D\_\{1\}^\{1/2\}AK^\{\\top\}D\_\{2\}KAD\_\{1\}^\{1/2\},\(69\)S2\\displaystyle S\_\{2\}=δ2​D21/2​K​A​D1​A​K⊤​D21/2\.\\displaystyle=\\,\\delta^\{2\}\\,D\_\{2\}^\{1/2\}KAD\_\{1\}AK^\{\\top\}D\_\{2\}^\{1/2\}\.\(70\)These operators are symmetric positive semidefinite\. Indeed, if we define

L:=δ​D21/2​K​A​D11/2,\\displaystyle L:=\\;\\delta\\;D\_\{2\}^\{1/2\}KAD\_\{1\}^\{1/2\},\(71\)then

S1=L⊤​L,S2=L​L⊤\.\\displaystyle S\_\{1\}=L^\{\\top\}L,\\qquad S\_\{2\}=LL^\{\\top\}\.\(72\)Now, let\(λr,ur\)\(\\lambda\_\{r\},u\_\{r\}\)denote an orthonormal eigensystem ofS2S\_\{2\}:

S2​ur=λr​ur\\displaystyle S\_\{2\}u\_\{r\}=\\lambda\_\{r\}u\_\{r\}\(73\)Supposeδ\\deltais small enough that0≤λr<10\\leq\\lambda\_\{r\}<1\. The nonzero eigenvalues ofS1S\_\{1\}andS2S\_\{2\}agree\. For the modal observable

ar​\(t\):=⟨ur,z2\(t\)⟩,\\displaystyle a\_\{r\}\(t\):=\\langle u\_\{r\},z\_\{2\}^\{\(t\)\}\\rangle,\(74\)we have the approximate relation

𝔼​\[ar​\(t\+1\)∣ar​\(t\)\]≈λr​ar​\(t\),\\displaystyle\\mathbb\{E\}\[a\_\{r\}\(t\+1\)\\mid a\_\{r\}\(t\)\]\\approx\\lambda\_\{r\}a\_\{r\}\(t\),\(75\)so the modal lag\-ggautocorrelation is roughly

ρr​\(g\)≈λrg\.\\displaystyle\\rho\_\{r\}\(g\)\\approx\\lambda\_\{r\}^\{g\}\.\(76\)Therefore the infinite integrated autocorrelation time of moderris

τr=1\+2​∑g=1∞λrg=1\+λr1−λr\.\\displaystyle\\tau\_\{r\}=1\+2\\sum\_\{g=1\}^\{\\infty\}\\lambda\_\{r\}^\{g\}=\\frac\{1\+\\lambda\_\{r\}\}\{1\-\\lambda\_\{r\}\}\.\(77\)Thus the eigenvalues have a direct dynamical meaning: smallλr\\lambda\_\{r\}correspond to nearly i\.i\.d\. modes, while eigenvalues close to11correspond to persistent slow modes\.

A useful metric of the spinwise concentration of a mode is the participation ratio

PR​\(ur\)=1∑iur​\(i\)4,\\displaystyle\\mathrm\{PR\}\(u\_\{r\}\)=\\frac\{1\}\{\\sum\_\{i\}u\_\{r\}\(i\)^\{4\}\},\(78\)which distinguishes localized modes from broad collective modes\. Broad slow modes are precisely the dangerous case for correlated many\-spin errors, while localized slow modes affect only a small subset of coordinates\.

#### C\.2\.1Dependence onδ\\delta

The operator depends onδ\\deltaboth explicitly and through the stationary susceptibilitiesD1​\(δ\),D2​\(δ\)D\_\{1\}\(\\delta\),D\_\{2\}\(\\delta\)\. Nearδ=0\\delta=0, however, the structure is especially simple\. Since the second bank is driven by a field scaled byδ\\delta, we havem2​\(δ\)=O​\(δ\)m\_\{2\}\(\\delta\)=O\(\\delta\), henceD2​\(δ\)=I\+O​\(δ2\)D\_\{2\}\(\\delta\)=I\+O\(\\delta^\{2\}\)\. The first bank has a nonzero forward field already atδ=0\\delta=0, som1​\(δ\)=m1​\(0\)\+O​\(δ2\)m\_\{1\}\(\\delta\)=m\_\{1\}\(0\)\+O\(\\delta^\{2\}\), henceD1​\(δ\)=D1​\(0\)\+O​\(δ2\)D\_\{1\}\(\\delta\)=D\_\{1\}\(0\)\+O\(\\delta^\{2\}\)\. It follows that

S1​\(δ\)\\displaystyle S\_\{1\}\(\\delta\)=δ2​S¯1\+O​\(δ4\),\\displaystyle=\\delta^\{2\}\\,\\bar\{S\}\_\{1\}\+O\(\\delta^\{4\}\),\(79\)S2​\(δ\)\\displaystyle S\_\{2\}\(\\delta\)=δ2​S¯2\+O​\(δ4\),\\displaystyle=\\delta^\{2\}\\,\\bar\{S\}\_\{2\}\+O\(\\delta^\{4\}\),\(80\)for appropriateδ\\delta\-independent limiting operatorsS¯1,S¯2\\bar\{S\}\_\{1\},\\bar\{S\}\_\{2\}\. Consequently each eigenvalue satisfies

λr​\(δ\)=δ2​λ¯r\+O​\(δ4\),\\displaystyle\\lambda\_\{r\}\(\\delta\)=\\delta^\{2\}\\bar\{\\lambda\}\_\{r\}\+O\(\\delta^\{4\}\),\(81\)and hence

τr​\(δ\)=1\+λr​\(δ\)1−λr​\(δ\)=1\+2​λ¯r​δ2\+O​\(δ4\)\.\\displaystyle\\tau\_\{r\}\(\\delta\)=\\frac\{1\+\\lambda\_\{r\}\(\\delta\)\}\{1\-\\lambda\_\{r\}\(\\delta\)\}=1\+2\\bar\{\\lambda\}\_\{r\}\\delta^\{2\}\+O\(\\delta^\{4\}\)\.\(82\)Thus the modal integrated autocorrelation times remain close to11asδ→0\\delta\\to 0\. In particular, there is no blow\-up of autocorrelation time in the low\-δ\\deltaregime\.

![Refer to caption](https://arxiv.org/html/2607.00170v1/figures/s2_mode_scatter_many.png)Figure 12:Predicted vs\. measured lag\-1 autocorrelation function for dominant eigenmodes of theL⊤​LL^\{\\top\}LandL​L⊤LL^\{\\top\}operators\.
#### C\.2\.2Experimental Validation on Modal Observables

The theory above was tested on its natural observables, namely the eigenmodes themselves\. For each retained modeuru\_\{r\}, the empirical red\-black Gibbs trajectory was projected ontouru\_\{r\}, and the resulting scalar mode process was used to estimate the lag\-1 autocorrelation\. When compared against the theoretical predictionλr\\lambda\_\{r\}, the agreement was very strong: across the top9696modes, the Pearson correlation between predicted and measured lag\-1 modal autocorrelation was approximately0\.9800\.980for banks1s\_\{1\}and approximately0\.9750\.975for banks2s\_\{2\}\. Weighted correlations were slightly higher still\. This is strong evidence that theL⊤​LL^\{\\top\}Ltheory captures the persistent same\-bank modal dynamics of full red\-black sweeps with high accuracy\.

It is important to emphasize that modal validation is the correct comparison for this theory\. Coordinatewise truncated integrated autocorrelation can be affected by short\-lag oscillatory effects and by the truncation of the empirical sum, whereas the present PSD theory is intentionally a theory of the persistent same\-bank envelope\. For the purpose of obtaining upper bounds on infinite integrated autocorrelation times, however, this is the correct level of description\.

#### C\.2\.3Per\-Spin Integrated Autocorrelation from Modes

Within the same linearized PSD model, the integrated autocorrelation time of a coordinate is a convex combination of the modal values\. Ifur​\(i\)u\_\{r\}\(i\)denotes theiith coordinate of the normalized eigenvectoruru\_\{r\}, then

τi=∑rur​\(i\)2​1\+λr1−λr\.\\displaystyle\\tau\_\{i\}=\\sum\_\{r\}u\_\{r\}\(i\)^\{2\}\\frac\{1\+\\lambda\_\{r\}\}\{1\-\\lambda\_\{r\}\}\.\(83\)This immediately implies the global bound

τi≤1\+λ11−λ1,\\displaystyle\\tau\_\{i\}\\leq\\frac\{1\+\\lambda\_\{1\}\}\{1\-\\lambda\_\{1\}\},\(84\)whereλ1\\lambda\_\{1\}is the top eigenvalue\. Hence the maximum per\-spin integrated autocorrelation time satisfies

τint∗≤1\+λ11−λ1\.\\displaystyle\\tau\_\{\\mathrm\{int\}\}^\{\*\}\\leq\\frac\{1\+\\lambda\_\{1\}\}\{1\-\\lambda\_\{1\}\}\.\(85\)Sinceλ1​\(δ\)=O​\(δ2\)\\lambda\_\{1\}\(\\delta\)=O\(\\delta^\{2\}\), we conclude that

τint∗​\(δ\)=1\+O​\(δ2\)\.\\displaystyle\\tau\_\{\\mathrm\{int\}\}^\{\*\}\(\\delta\)=1\+O\(\\delta^\{2\}\)\.\(86\)This is the central theoretical conclusion needed for the convergence analysis: in the low\-δ\\deltaregime, the difficulty of sign recovery is controlled primarily by the smallness of the meansμi\\mu\_\{i\}, not by any growth in the integrated autocorrelation time\.

A sharper coordinatewise upper bound is available from the topkkmodes\. Define

ϕ​\(λ\):=1\+λ1−λ\.\\displaystyle\\phi\(\\lambda\):=\\frac\{1\+\\lambda\}\{1\-\\lambda\}\.\(87\)Then monotonicity ofϕ\\phigives

τi≤∑r=1kur​\(i\)2​ϕ​\(λr\)\+\(1−∑r=1kur​\(i\)2\)​ϕ​\(λk\)\.\\displaystyle\\tau\_\{i\}\\leq\\sum\_\{r=1\}^\{k\}u\_\{r\}\(i\)^\{2\}\\phi\(\\lambda\_\{r\}\)\+\\Big\(1\-\\sum\_\{r=1\}^\{k\}u\_\{r\}\(i\)^\{2\}\\Big\)\\phi\(\\lambda\_\{k\}\)\.\(88\)Thus the top modes do not merely explain the large\-scale dynamics qualitatively; they provide explicit rigorous upper bounds on every per\-spin integrated autocorrelation time within the linearized theory\.

#### C\.2\.4Consequences for Convergence of Time Averages

The purpose of introducingτi,int\\tau\_\{i,\\mathrm\{int\}\}is to control the effective sample size of the time average of spinii\. If

mi\(n\):=1n​∑t=1nσi\(t\),\\displaystyle m\_\{i\}^\{\(n\)\}:=\\frac\{1\}\{n\}\\sum\_\{t=1\}^\{n\}\\sigma\_\{i\}^\{\(t\)\},\(89\)then under the convention

τi,int=1\+2​∑g=1∞ρi​\(g\),\\displaystyle\\tau\_\{i,\\mathrm\{int\}\}=1\+2\\sum\_\{g=1\}^\{\\infty\}\\rho\_\{i\}\(g\),\(90\)the asymptotic variance of the time average is

Var⁡\(mi\(n\)\)≈\(1−μi2\)​τi,intn\.\\displaystyle\\operatorname\{Var\}\(m\_\{i\}^\{\(n\)\}\)\\approx\\frac\{\(1\-\\mu\_\{i\}^\{2\}\)\\tau\_\{i,\\mathrm\{int\}\}\}\{n\}\.\(91\)Therefore the number of sweeps required to determine the correct sign ofmi\(n\)m\_\{i\}^\{\(n\)\}with high confidence is controlled by

n≳1−μi2μi2​\(confidence factor\)​τi,int\.\\displaystyle n\\gtrsim\\frac\{1\-\\mu\_\{i\}^\{2\}\}\{\\mu\_\{i\}^\{2\}\}\\,\(\\text\{confidence factor\}\)\\,\\tau\_\{i,\\mathrm\{int\}\}\.\(92\)Using the modal upper bound yields the sufficient condition

n≳1−μi2μi2​\(confidence factor\)​1\+λ11−λ1\.\\displaystyle n\\gtrsim\\frac\{1\-\\mu\_\{i\}^\{2\}\}\{\\mu\_\{i\}^\{2\}\}\\,\(\\text\{confidence factor\}\)\\,\\frac\{1\+\\lambda\_\{1\}\}\{1\-\\lambda\_\{1\}\}\.\(93\)In particular, sinceλ1​\(δ\)=O​\(δ2\)\\lambda\_\{1\}\(\\delta\)=O\(\\delta^\{2\}\), the autocorrelation correction stays bounded and close to11asδ→0\\delta\\to 0\. The dominant low\-δ\\deltascaling therefore comes from the factorμi−2\\mu\_\{i\}^\{\-2\}, not from any growth of autocorrelation\.

#### C\.2\.5Summary

The red\-black full\-sweep dynamics admits a natural linearization in normalized coordinates whose same\-bank persistent dynamics is governed by the symmetric PSD operatorsS1=L⊤​LS\_\{1\}=L^\{\\top\}LandS2=L​L⊤S\_\{2\}=LL^\{\\top\}\. Their eigenvalues determine modal autocorrelation times through

τr=1\+λr1−λr,\\displaystyle\\tau\_\{r\}=\\frac\{1\+\\lambda\_\{r\}\}\{1\-\\lambda\_\{r\}\},\(94\)and each per\-spin integrated autocorrelation time is a convex combination of these modal values\. Hence

τi≤1\+λ11−λ1,\\displaystyle\\tau\_\{i\}\\leq\\frac\{1\+\\lambda\_\{1\}\}\{1\-\\lambda\_\{1\}\},\(95\)with the sharper top\-kkrefinement above\. Empirically, as seen in[Figure12](https://arxiv.org/html/2607.00170#A3.F12), the modal lag\-1 autocorrelations are predicted with Pearson correlation approximately0\.980\.98, validating the theory at the level of its natural observables\. Sinceλ1​\(δ\)=O​\(δ2\)\\lambda\_\{1\}\(\\delta\)=O\(\\delta^\{2\}\), the resulting autocorrelation bound remains close to11in the low\-δ\\deltaregime\.