1 Introduction↩︎

Thermodynamic computing has recently re-emerged as a promising physical route for information processing beyond conventional deterministic digital architectures [1]–[12]. Thermodynamic computers do not regard thermal fluctuations as sources of error to be suppressed, but rather seek to transform the random motion itself into a computational resource. Along this direction, Whitelam and Casert introduced a nonlinear thermodynamic computing framework: within this framework, fluctuating degrees of freedom constrained by a quartic potential and coupled to a thermal bath act as thermodynamic neurons capable of performing nonlinear computations [13].

In the Whitelam–Casert paradigm, the elementary computing unit is a continuously fluctuating physical coordinate whose activity is shaped by a nonlinear energy landscape. A driven quartic potential converts an external input into an input-dependent stochastic response, so that the observed state of the system at a finite time \(t_f\) implements an effective nonlinear activation. This idea provides a physically grounded alternative to abstract artificial neurons: the activation function is not prescribed analytically, but emerges from the nonequilibrium evolution of a thermally driven dynamical system. The resulting architecture establishes a direct connection between stochastic thermodynamics, nonlinear dynamical systems, and neural computation.

However, the representational capacity of such a Langevin computer depends critically on how the fluctuating state is read out. If the final-time distribution of a thermodynamic neuron is reduced only to its mean displacement, the readout corresponds to a first-order coarse graining of the conditional probability density \(P(x,t_f|I)\). This is sufficient to define a deterministic activation-like response, but it discards much of the physical information carried by the nonequilibrium ensemble. In a nonlinear potential, the input does not merely shift the center of the distribution; it can also reshape its width, tail weight, asymmetry, and non-Gaussian structure. These changes are especially important for a quartic landscape, where the fourth-order confinement term directly controls the high-amplitude excursions of the stochastic coordinate.

This observation motivates the central question of the present work: can a nonequilibrium thermodynamic computer be made more expressive by reading out not only the mean response, but an ordered set of distributional response statistics? We answer this question by introducing moment-resolved readout, a moment-based extension of Langevin computing. For each driven thermodynamic reservoir, instead of using only the first moment \(\mathbb{E}[\boldsymbol{x}]\), we construct an energy-aligned polynomial response vector \[\phi_{\boldsymbol{K}}(\boldsymbol{I}) = \left[ \mathbb{E}[\boldsymbol{x}], \mathbb{E}[\boldsymbol{x}^{\odot 2}], \mathbb{E}[\boldsymbol{x}^{\odot 4}] \right],\] where \(\boldsymbol{x}^{\odot p}\) denotes the elementwise \(p\)-th power. The first moment captures the mean displacement response. The second and fourth raw moments provide progressively higher-order polynomial summaries of the finite-time distribution and are naturally aligned with the quadratic and quartic terms of the local confinement. Because they are raw rather than central moments, however, they mix changes in the distribution center with changes in its width, asymmetry, and higher-order shape. We therefore interpret them as energy-aligned polynomial response observables, rather than as pure measures of variance or non-Gaussianity.

Recent progress has also shown that Langevin thermodynamic computers can be used for generative modeling, where structured data are synthesized from noise by the natural time evolution of a trained thermodynamic system [14]. In that setting, the central quantity is the probability of generating the reverse of a noising trajectory, and the learning rule admits a physical interpretation in terms of heat emission and entropy production. This generative perspective demonstrates that nonequilibrium Langevin dynamics can encode data-producing transformations in an energy landscape. The present work addresses a complementary discriminative problem. Rather than asking how a thermodynamic computer can generate structured samples from noise, we ask how a finite-time nonequilibrium distribution should be read when the task is pattern recognition.

The second limitation addressed in this paper concerns the use of a single physical reservoir. A single Langevin network can generate rich nonlinear responses, but the geometry of its responses is constrained by fixed input projections, coupling topology, training pathways, and noise scale [15]. Simply increasing the size of a homogeneous reservoir does not guarantee a corresponding increase in useful computational diversity, as additional coordinates may explore regions of the response space that are highly correlated [16], [17]. For physical computing systems, this can lead to redundancy: two reservoirs may perform excellently as independent classifiers, but if their error patterns are too similar, combining them will not improve generalization ability.

To overcome this limitation, we introduce a heterogeneous multi-reservoir architecture [3], [18]–[21]. Instead of using identical replicas of one trained thermodynamic computer, we construct three reservoirs with different dynamical origins: a reservoir \(K_1\) trained directly using the polynomial-response surrogate, a second reservoir \(K_2\) initialized from a mean-response-trained configuration and subsequently refined using the polynomial surrogate, and a third reservoir \(K_3\) initialized independently with a stronger internal coupling scale [22], [23]. These reservoirs are designed to generate partially distinct response bases under the same input stimulus. The final representation is the feature-level concatenation \[\Phi(\boldsymbol{I}) = \left[ \phi_{\boldsymbol{K}_1}(\boldsymbol{I}), \phi_{\boldsymbol{K}_2}(\boldsymbol{I}), \phi_{\boldsymbol{K}_3}(\boldsymbol{I}) \right],\] followed by a single standardized linear readout. The purpose of this construction is to examine whether reservoirs with distinct dynamical and training histories can provide complementary response coordinates and classification errors. Because the multi-reservoir model also increases feature dimension and computational resources, the present comparison should not by itself be interpreted as a complete separation of heterogeneity effects from capacity or sampling-budget effects. This work therefore extends nonlinear thermodynamic computing out of equilibrium in two complementary directions. First, it replaces mean-only readout with moment-resolved readout. Second, it replaces a single-reservoir representation with a heterogeneous multi-reservoir response basis, allowing weaker reservoirs with partially complementary error patterns to contribute to the joint representation [24]–[26].

2 Methodology↩︎

2.1 Stochastic Thermodynamics of Quartic-Confined Langevin Neurons↩︎

We model the thermodynamic computer as a collection of \(N\) coupled physical nodes driven by nonequilibrium stochastic dynamics. Given an external pattern vector \(\boldsymbol{I}\in\mathbb{R}^{D}\), the microstate of the \(r\)-th reservoir, \(r\in\{1,\dots,R\}\), is described by the time-dependent vector \(\boldsymbol{x}(t)=[x_1(t),x_2(t),\dots,x_N(t)]^\top\). Each nonlinear degree of freedom evolves according to the overdamped Langevin equation [27], [28] \[\label{eq:langevin} dx_i=\mu F_i(\boldsymbol{x},\boldsymbol{I})\,dt+\sqrt{2\mu k_B T}\,dB_i(t),\tag{1}\] where \(\mu\) is the mobility parameter setting the intrinsic relaxation time scale, \(k_B T\) denotes the thermal noise intensity of the background heat bath, and \(B_i(t)\) are mutually independent standard Brownian motions. The deterministic force acting on the \(i\)-th node combines the local quartic confinement, internal coupling, static bias, and external input projection: \[\label{eq:force} F_i(\boldsymbol{x},\boldsymbol{I}) = -2J_2x_i-4J_4x_i^3 +b_i+\sum_{j=1}^{N}W_{ij}x_j +\sum_{k=1}^{D}K_{ik}I_k .\tag{2}\] Here \(J_2\) and \(J_4\) specify the quadratic and stabilizing quartic components of the intrinsic potential, respectively. The matrices \(\boldsymbol{W}\) and \(\boldsymbol{K}\) define the internal recurrent interaction topology and the external input projection, while \(b_i\) is a local static bias. For an isolated node, or when discussing only the onsite contribution to the coupled dynamics, the input and bias act through the effective scalar drive \[h_i(\boldsymbol{I})=b_i+\sum_{k=1}^{D}K_{ik}I_k,\] leading to the local quartic potential \[\label{eq:local95potential} U_i(x_i;h_i)=J_2x_i^2+J_4x_i^4-h_i x_i .\tag{3}\]

In the numerical implementation, the recurrent matrix is directed and is not constrained to be symmetric. At initialization, its off-diagonal elements are independently sampled according to \[W_{ij}\sim \mathcal{N}\!\left(0,\frac{s_W^2}{N}\right), \qquad i\neq j, \qquad W_{ii}=0 . \label{eq:w95initialization}\tag{4}\] The input and bias parameters are initialized as \[K_{ik}\sim \mathcal{N}\!\left(0,\frac{s_K^2}{D}\right), \qquad b_i\sim\mathcal{N}(0,s_b^2), \label{eq:kb95initialization}\tag{5}\] with the values of \(s_W\), \(s_K\), and \(s_b\) specified in the numerical protocol below.

A global scalar potential for the coupled system would require the integrability condition \[\frac{\partial F_i}{\partial x_j} = \frac{\partial F_j}{\partial x_i}, \qquad\text{equivalently}\qquad W_{ij}=W_{ji}. \label{eq:integrability95condition}\tag{6}\] If \(\boldsymbol{W}=\boldsymbol{W}^{\mathsf T}\), the deterministic force can be derived from \[U(\boldsymbol{x};\boldsymbol{I}) = \sum_{i=1}^{N} \left( J_2x_i^2+J_4x_i^4-h_i(\boldsymbol{I})x_i \right) - \frac{1}{2}\boldsymbol{x}^{\mathsf T}\boldsymbol{W}\boldsymbol{x}. \label{eq:global95symmetric95potential}\tag{7}\] The matrices used here are generally asymmetric, and therefore the full coupled dynamics are not assumed to derive from a global scalar energy landscape. In this paper, the term “energy-sensitive” refers specifically to local polynomial observables aligned with the onsite quadratic and quartic confinement, while the recurrent interaction acts as a generally nonconservative nonequilibrium drive.

At the ensemble level, the finite-time probability density \(P(\boldsymbol{x},t|\boldsymbol{I})\) obeys the associated Fokker–Planck equation [29] \[\label{eq:fokker95planck} \frac{\partial P}{\partial t} = -\mu\sum_{i=1}^{N} \frac{\partial}{\partial x_i} \left[ F_i(\boldsymbol{x},\boldsymbol{I})P \right] + \mu k_B T \sum_{i=1}^{N} \frac{\partial^2 P}{\partial x_i^2}.\tag{8}\] Because the physical network is observed at a finite readout time \(t_f\), we do not assume that \(P(\boldsymbol{x},t_f|\boldsymbol{I})\) has relaxed to an equilibrium Boltzmann distribution.

Rather than representing this finite-time nonequilibrium distribution only by the first moment \(\mathbb{E}[\boldsymbol{x}]\), we construct a finite-dimensional moment-resolved response representation from raw moment observables. For the coupled nonlinear Langevin system considered here, an exact analytical expression for the transient density \(P(\boldsymbol{x},t_f|\boldsymbol{I})\) is generally unavailable. We therefore approximate its response statistics by Monte Carlo sampling of the stochastic dynamics. The theoretical \(p\)-th raw moment of the \(i\)-th coordinate is \[\label{eq:theoretical95moment} \mathbb{E}\!\left[x_i^p(t_f)\mid \boldsymbol{I}\right] = \int x_i^p P(\boldsymbol{x},t_f|\boldsymbol{I})\,d\boldsymbol{x}.\tag{9}\] Using \(M\) independently reset trajectories initialized from \(\boldsymbol{x}(0)=\boldsymbol{0}\), this expectation is estimated by the empirical ensemble average \[\label{eq:empirical95moment} m_{p,i}(\boldsymbol{I}) = \frac{1}{M}\sum_{m=1}^{M} \left[ x_i^{(m)}(t_f;\boldsymbol{I}) \right]^p, \qquad p\in\{1,2,4\}.\tag{10}\] This estimator is the standard Monte Carlo approximation to the corresponding finite-time moment of the Langevin ensemble; its sampling uncertainty decreases with the number of independent replicas.

The choice \(p\in\{1,2,4\}\) provides a low-order polynomial basis aligned with the linear drive and the quadratic and quartic terms of the onsite confinement. To state the physical meaning precisely, let \[\begin{align} \bar{x}_i &= \mathbb{E}[x_i], & \sigma_i^2 &= \mathbb{E}\!\left[ (x_i-\bar{x}_i)^2 \right], \\ \mu^{\mathrm c}_{q,i} &= \mathbb{E}\!\left[ (x_i-\bar{x}_i)^q \right]. \end{align} \label{eq:central95moment95definitions}\tag{11}\] The second raw moment satisfies \[m_{2,i} = \mathbb{E}[x_i^2] = \sigma_i^2+\bar{x}_i^2, \label{eq:raw95second95decomposition}\tag{12}\] and the fourth raw moment can be decomposed as \[m_{4,i} = \mathbb{E}[x_i^4] = \mu^{\mathrm{c}}_{4,i} + 4\bar{x}_i\mu^{\mathrm{c}}_{3,i} + 6\bar{x}_i^2\sigma_i^2 + \bar{x}_i^4. \label{eq:raw95fourth95decomposition}\tag{13}\] All moments in Eqs. 11 – 13 are evaluated at \(t_f\) and conditioned on the input \(\boldsymbol{I}\); these arguments are suppressed for compactness.

Consequently, \(m_{2,i}\) is not a pure variance observable, and \(m_{4,i}\) is not a pure measure of kurtosis or non-Gaussian tail weight. Instead, they are raw polynomial observables that combine displacement and central-shape information. Their relevance here follows from their direct alignment with the quadratic and quartic powers appearing in the local dynamics. Higher raw moments could in principle be included, but their Monte Carlo estimators become increasingly sensitive to rare samples and generally exhibit larger finite-\(M\) variance.

For a reservoir with input projection \(\boldsymbol{K}\), the local moment-resolved response representation is then defined as \[\label{eq:moment95representation} \phi_{\boldsymbol{K}}(\boldsymbol{I}) = \left[ \boldsymbol{m}_1(\boldsymbol{I}), \boldsymbol{m}_2(\boldsymbol{I}), \boldsymbol{m}_4(\boldsymbol{I}) \right] \in\mathbb{R}^{3N}.\tag{14}\] Equivalently, \(\boldsymbol{m}_p(\boldsymbol{I})\) denotes the elementwise \(p\)-th raw moment vector of the finite-time response distribution. This construction maps the transient input-conditioned density \(P(\boldsymbol{x},t_f|\boldsymbol{I})\) to a finite-dimensional moment-resolved response representation for subsequent classification. The mapping is intentionally incomplete: it retains selected nodewise polynomial moments but does not reconstruct the full probability density or the inter-node correlations \(\mathbb{E}[x_i x_j]\).

2.2 Nonequilibrium Trajectory Dynamics and Finite-Time Moment Readout↩︎

As shown in Fig. 1, the quartic-confined Langevin neuron converts an input drive into a tilted nonlinear energy landscape and finite-time stochastic trajectories. In the quartic-confined single-well response regime, the deterministic component of the force tends to relax the coordinate toward the global minimum, whereas the thermal bath continuously perturbs the trajectory away from purely deterministic descent. Consequently, different reset-sampled stochastic replicas undergo finite-time thermal excursions that reshape the nonequilibrium distribution within the observation time \(t_f\), as illustrated schematically in Fig. 1 (b).

In this work, “nonequilibrium” is used in the operational finite-time sense: the system is reset to \(\boldsymbol{x}(0)=\boldsymbol{0}\), driven by an input-dependent force, and read at a prescribed time \(t_f\) without assuming relaxation to a stationary distribution. The present classification experiment does not by itself quantify the distance between \(P(\boldsymbol{x},t_f|\boldsymbol{I})\) and the corresponding long-time stationary distribution. Accordingly, the term finite-time nonequilibrium response refers to the reset-and-read protocol rather than to a separately established steady-state entropy-production regime. Let \(\mathcal{L}_{\boldsymbol{I}}^{\dagger}\) denote the Fokker–Planck operator associated with Eq. 8. For the reset-sampling protocol used in this work, the input-conditioned probability density can be written formally as \[\label{eq:propagator} P(\boldsymbol{x},t_f|\boldsymbol{I}) = \exp\!\left(t_f\mathcal{L}_{\boldsymbol{I}}^{\dagger}\right) \delta(\boldsymbol{x}),\tag{15}\] where the delta distribution represents the reset initial condition \(\boldsymbol{x}(0)=\boldsymbol{0}\). This expression emphasizes that the readout probes the finite-time propagator generated by the driven Langevin dynamics. It does not require the system to have relaxed to a stationary Boltzmann distribution at the observation time.

Numerically, each stochastic replica is advanced using an Euler–Maruyama discretization of Eq. 1  [30]: \[\label{eq:euler95maruyama} x_{i,n+1}^{(m)} = x_{i,n}^{(m)} + \mu F_i(\boldsymbol{x}_{n}^{(m)},\boldsymbol{I})\Delta t + \sqrt{2\mu k_B T\Delta t}\, \xi_{i,n}^{(m)},\tag{16}\] where \(\xi_{i,n}^{(m)}\sim\mathcal{N}(0,1)\) are independent Gaussian random variables across node index \(i\), time step \(n\), and replica index \(m\). The final-time samples \(\{\boldsymbol{x}^{(m)}(t_f;\boldsymbol{I})\}_{m=1}^{M}\) provide a Monte Carlo representation of \(P(\boldsymbol{x},t_f|\boldsymbol{I})\).

For numerical safety, the updated state is projected after each Euler–Maruyama step onto the interval \[x_{i,n+1}^{(m)} \leftarrow \Pi_{[-x_{\mathrm c},x_{\mathrm c}]} \!\left( x_{i,n+1}^{(m)} \right), \label{eq:hard95boundary}\tag{17}\] where \[\Pi_{[-x_{\mathrm c},x_{\mathrm c}]}(y) = \min \left\{ x_{\mathrm c}, \max(-x_{\mathrm c},y) \right\}. \label{eq:hard95boundary95operator}\tag{18}\] where \(x_{\mathrm c}=4\) in the reported calculations. The implementation records the number of attempted updates that exceed this interval, so that the influence of the numerical boundary can be assessed independently from the classification result.

The statistical readout of these finite-time trajectories is summarized separately in Fig. 2. Rather than compressing the distribution into the first-order response alone, the proposed moment-resolved thermodynamic response representation uses raw moment channels to sample distinct physical aspects of the nonequilibrium ensemble.

The first moment \(\mathbb{E}[x_i]\) measures the mean displacement response. The second raw moment \(\mathbb{E}[x_i^2]\) combines the variance with the squared mean displacement, while the fourth raw moment \(\mathbb{E}[x_i^4]\) combines fourth-order central-shape information with lower-order location and fluctuation contributions, as shown in Eqs. 12 and 13 . These observables are therefore interpreted as progressively higher-order polynomial summaries of the finite-time response, rather than as pure variance and kurtosis measurements. The response vector in Eq. 14 retains information beyond the mean while remaining directly aligned with the polynomial structure of the onsite force [31].

3 Results and Discussion↩︎

3.1 Experimental Protocol↩︎

We evaluate the proposed models on the standard MNIST handwritten-digit classification benchmark. The official training set contains \(60000\) samples and the official test set contains \(10000\) samples. The results reported here follow the fixed \(60000/10000\) protocol used during model development: reservoir construction and final readout fitting use the complete training set, and the reported accuracies are evaluated on the official test set.

During an earlier exploratory stage, a \(2000\)-sample subset of the official test set was inspected for preliminary comparisons among feature combinations and fusion settings. Consequently, the reported \(10000\)-sample results should be interpreted as fixed-protocol reproduction results rather than as an untouched confirmatory estimate. For the final public reproduction, all numerical parameters, random seeds, reservoir definitions, feature channels, and readout hyperparameters will be fixed before regenerating the complete \(60000/10000\) results. No subsequent parameter adjustment will be made from the reproduced test predictions.

Each \(28\times28\) grayscale image is flattened into a vector \(\boldsymbol{I}\in\mathbb{R}^{784}\), with pixel intensities scaled according to \[I_k=\frac{u_k}{255}, \qquad u_k\in\{0,1,\ldots,255\}. \label{eq:input95scaling}\tag{19}\] No pixelwise test-set statistics are used. Standardization is applied only to the extracted response features, using statistics computed from the training feature matrix.

The finite-time Langevin simulations use \[\begin{align} N&=256, & M&=128, & J_2&=1, & J_4&=1, \\ \mu&=1, & k_BT&=0.1. \end{align} \label{eq:primary95dynamics95parameters}\tag{20}\] together with \[\begin{align} \Delta t&=10^{-3}, & t_f&=1, & N_t&=\frac{t_f}{\Delta t}=1000, \\ \boldsymbol{x}(0)&=\boldsymbol{0}. \end{align} \label{eq:integration95parameters}\tag{21}\] A hard numerical boundary \(x_{\mathrm c}=4\) is applied according to Eq. 18 . The primary stochastic trajectory seed is \(20260609\). Noise streams are generated independently across input samples, reservoir index, replica index, node index, and integration step. The same extracted feature realization is reused across nested moment-order ablations, so that the comparisons between \(\boldsymbol{m}_1\), \((\boldsymbol{m}_1,\boldsymbol{m}_2)\), and \((\boldsymbol{m}_1,\boldsymbol{m}_2,\boldsymbol{m}_4)\) are paired at the feature-extraction level.

All reservoirs use \(N=256\) thermodynamic nodes. For each input image, \(M=128\) independent Langevin trajectories are sampled to estimate the raw response moments \(m_1\), \(m_2\), and \(m_4\). Thus, each single-reservoir moment-resolved response representation has dimension \(3N=768\). The full heterogeneous multi-reservoir response representation concatenates three reservoir blocks and has dimension \(9N=2304\).

The three reservoirs are constructed using fixed but distinct initialization and surrogate-training histories. The recurrent matrices remain fixed during input-projection training. Reservoirs \(K_1\) and \(K_2\) use \(s_W=0.1\), whereas \(K_3\) uses \(s_W=0.2\). All input projections are initialized with \(s_K=1.2\), and the local biases are initialized with \(s_b=0.02\).

The input projections are trained using deterministic local-response surrogates rather than by differentiating through the full stochastic Langevin simulation. For an input \(\boldsymbol{I}\), the local surrogate response \(x_i^\star\) is defined as the real stationary solution of \[\begin{align} 2J_2x_i^\star + 4J_4(x_i^\star)^3 &= h_i(\boldsymbol{I}), \\ h_i(\boldsymbol{I}) &= b_i+\sum_{k=1}^{D}K_{ik}I_k. \end{align} \label{eq:local95surrogate95fixed95point}\tag{22}\] The mean-response surrogate uses \[\boldsymbol{\psi}_{\mathrm{mean}}(\boldsymbol{I}) = \boldsymbol{x}^\star(\boldsymbol{I}), \label{eq:mean95surrogate}\tag{23}\] whereas the polynomial surrogate uses \[\boldsymbol{\psi}_{\mathrm{poly}}(\boldsymbol{I}) = \left[ \boldsymbol{x}^\star(\boldsymbol{I}), (\boldsymbol{x}^\star(\boldsymbol{I}))^{\odot2}, (\boldsymbol{x}^\star(\boldsymbol{I}))^{\odot4} \right]. \label{eq:polynomial95surrogate}\tag{24}\] These surrogates do not include thermal noise, finite-time integration, or recurrent coupling and are used only to initialize and optimize the input projections. All classification features reported in Tables ¿tbl:tab:ablation?–¿tbl:tab:moment95order95ablation? are subsequently regenerated using the full finite-time coupled Langevin dynamics of Eq. 1 .

Reservoir \(K_1\) is trained directly with the polynomial surrogate for \(60\) epochs using learning rate \(10^{-3}\) and seed \(20260611\). Reservoir \(K_2\) is first trained with the mean-response surrogate for \(100\) epochs using learning rate \(3\times10^{-3}\), and is then refined for \(60\) epochs with the polynomial surrogate using learning rate \(10^{-3}\) and seed \(20260612\). Reservoir \(K_3\) is independently initialized and trained with the polynomial surrogate for \(80\) epochs using learning rate \(10^{-3}\), \(s_W=0.2\), and seed \(20260613\). The surrogate-training batch size is \(128\), and the regularization coefficient on the trained input projection is \(10^{-5}\).

The main comparison includes single-reservoir representations, pairwise concatenated representations, and the proposed three-reservoir feature-level model. Equal-weight decision-level logit averaging is included as a fixed diagnostic baseline. This comparison evaluates one simple form of late fusion and is not intended as an exhaustive comparison with calibrated, validation-weighted, or learned stacking methods. We denote a single-reservoir moment-resolved representation by SR-MRR, a pairwise concatenated representation by PR-MRR, the heterogeneous multi-reservoir feature-level model by HMR-FL, and the equal-weight decision-level logit-fusion baseline by HMR-LF.

The full HMR-FL model reduces the number of full-test errors from 349, 318, and 397 for \(K_1\), \(K_2\), and \(K_3\), respectively, to 305.

3.2 Error-Set Complementarity in Heterogeneous Reservoirs↩︎

The individual and pairwise results are summarized in Table ¿tbl:tab:ablation?. Among the single-reservoir models, \(K_2\) gives the highest observed test accuracy, \(96.82\%\). Concatenating response features does not produce a monotonic improvement: the \(K_1+K_2\) model reaches \(96.58\%\), below the standalone \(K_2\) result. This observation shows that increasing feature dimension alone is not sufficient to guarantee improved test performance under the fixed readout protocol.

The wrong-set Jaccard coefficient provides a descriptive measure of overlap between the samples misclassified by two models [23]–[26]. The \(K_1\)–\(K_2\) pair has the largest observed wrong-set overlap, \(J=0.7015\), whereas the \(K_2\)–\(K_3\) pair has the lowest, \(J=0.5376\). These values indicate differences in classification-error overlap, but they do not directly measure feature-space correlation, canonical angles, or linear collinearity. The reduced \(K_1+K_2\) accuracy is therefore consistent with partially redundant classification behavior, but the Jaccard analysis alone does not establish redundancy as its cause. Feature dimensionality, optimization, and regularization may also contribute.

Relative to \(K_2\), the full \(K_1+K_2+K_3\) feature model corrects \(68\) samples that \(K_2\) misclassifies and introduces \(55\) errors on samples that \(K_2\) classifies correctly. The exact two-sided McNemar test on these discordant pairs gives \(p\simeq0.279\). The \(0.13\)-percentage-point difference is therefore reported as the best observed accuracy among the compared configurations, not as a statistically established improvement over \(K_2\).

Although \(K_3\) is the weakest standalone model, its error set overlaps less with that of \(K_2\) than does the error set of \(K_1\). The full feature-level model reaches \(96.95\%\), suggesting that a weaker reservoir may still provide useful coordinates when its classification errors differ from those of stronger reservoirs. This relationship is descriptive rather than causal: wrong-set overlap does not demonstrate that the corresponding feature blocks are statistically decorrelated.

3.3 Moment-Order Ablation of Thermodynamic Response Representations↩︎

To isolate the physical contribution of the different response channels, we perform a formal moment-order ablation, as shown in Table ¿tbl:tab:moment95order95ablation?, using the same computational readout pipeline, optimization schedule, and \(L_2\) penalty as the main HMR-FL result.

The mean-only representation reaches \(93.77\%\). Adding the second raw moment increases the observed accuracy to \(96.82\%\), while adding the fourth raw moment gives \(96.95\%\). Because \(m_2=\sigma^2+\bar{x}^2\), the improvement obtained by adding \(\boldsymbol{m}_2\) cannot be attributed exclusively to fluctuation variance: the channel also contains nonlinear information about the squared mean displacement. Likewise, \(\boldsymbol{m}_4\) combines fourth-order central-shape information with lower-order location and fluctuation contributions.

The difference between \((\boldsymbol{m}_1,\boldsymbol{m}_2)\) and \((\boldsymbol{m}_1,\boldsymbol{m}_2,\boldsymbol{m}_4)\) is \(13\) correctly classified samples, or \(0.13\) percentage points, in the reported run. This is a small additional observed improvement and is not described here as a statistically established or systematic effect. The ablation supports the usefulness of progressively higher-order raw polynomial observables under the fixed protocol, but further trajectory and readout repetitions are needed to distinguish a robust fourth-order contribution from finite-sampling variability.

3.4 Information Bottleneck in Decision-Level Fusion: Feature-Level Fusion vs. Logits Ensembling↩︎

To further assess whether cross-reservoir information is better exploited before or after class-level compression, we compare the proposed unified feature-level readout with an equal-weight decision-level logits fusion baseline. In the decision-level setting, each reservoir response block is first mapped by an independent local readout to a \(10\)-dimensional class-logit vector, and the final logits are obtained as \[\label{eq:late95fusion} \boldsymbol{z}_{\rm LF} = \frac{\boldsymbol{z}_1+\boldsymbol{z}_2+\boldsymbol{z}_3}{3}.\tag{25}\] This late-fusion baseline achieves a full-test accuracy of \(96.84\%\), slightly above the strongest single-reservoir model but below the \(96.95\%\) accuracy obtained by the proposed feature-level fusion model.

The observed difference is consistent with an early compression constraint, although the comparison does not establish that compression is the sole cause. In decision-level fusion, each \(768\)-dimensional reservoir block is first mapped by an independently trained readout to a \(10\)-dimensional logit vector. The final averaging operation can use only these compressed class-level representations. In feature-level fusion, a single classifier is instead optimized over all \(2304\) input features simultaneously.

Importantly, the feature-level classifier remains linear and does not explicitly form cross-reservoir products or covariance features. Its advantage is therefore not that it directly evaluates quantities such as \(\phi_{K_r,i}\phi_{K_s,j}\), but that the class-feature weights of all reservoir blocks are optimized jointly under a common loss before reservoir-wise compression.

To formalize the feature-level construction, the response blocks extracted from the three biased dynamical configurations are concatenated into \[\label{eq:joint95representation} \begin{align} \boldsymbol{\Phi}(\boldsymbol{I}) &= \left[ \phi_{\boldsymbol{K}_1}(\boldsymbol{I}), \phi_{\boldsymbol{K}_2}(\boldsymbol{I}), \phi_{\boldsymbol{K}_3}(\boldsymbol{I}) \right] \in \mathbb{R}^{9N},\\ N&=256,\qquad 9N=2304 . \end{align}\tag{26}\] Before linear classification, the joint representation is standardized using training-set statistics: \[\begin{align} \mu_j &= \frac{1}{S} \sum_{s=1}^{S} \Phi_j(\boldsymbol{I}_s), \\ \widehat{\sigma}_j &= \left[ \frac{1}{S-1} \sum_{s=1}^{S} \left( \Phi_j(\boldsymbol{I}_s)-\mu_j \right)^2 \right]^{1/2} \end{align} \label{eq:training95feature95statistics}\tag{27}\] and \[s_j = \begin{cases} \widehat{\sigma}_j, & \widehat{\sigma}_j\geq10^{-12}, \\[3pt] 1, & \widehat{\sigma}_j<10^{-12}, \end{cases} \label{eq:feature95scale}\tag{28}\] and \[\widetilde{\Phi}_j(\boldsymbol{I}) = \frac{\Phi_j(\boldsymbol{I})-\mu_j}{s_j}. \label{eq:standardization}\tag{29}\] All feature means and standard deviations are computed from the training feature matrix only and are subsequently applied unchanged to the test features. The fallback \(s_j=1\) prevents division by a numerically constant feature channel.

The classification logits are computed by a single linear readout, \[\boldsymbol{z}(\boldsymbol{I}) = \boldsymbol{W}_{\mathrm{out}} \widetilde{\boldsymbol{\Phi}}(\boldsymbol{I}) + \boldsymbol{g}_{\mathrm{out}}, \label{eq:linear95readout}\tag{30}\] where \(\boldsymbol{W}_{\mathrm{out}}\in\mathbb{R}^{C\times2304}\) is the readout matrix, \(\boldsymbol{g}_{\mathrm{out}}\in\mathbb{R}^{C}\) is the bias vector, and \(C=10\) is the number of classes. The predicted label is \[\widehat y(\boldsymbol{I}) = \operatorname*{arg\,max}_{c\in\{1,\ldots,C\}} z_c(\boldsymbol{I}). \label{eq:prediction}\tag{31}\]

Partitioning the readout matrix according to the three reservoir blocks, \[\boldsymbol{W}_{\mathrm{out}} = \left[ \boldsymbol{W}_1,\boldsymbol{W}_2,\boldsymbol{W}_3 \right], \label{eq:partitioned95readout95matrix}\tag{32}\] gives \[\boldsymbol{z}(\boldsymbol{I}) = \sum_{r=1}^{3} \boldsymbol{W}_r \widetilde{\phi}_{\boldsymbol{K}_r}(\boldsymbol{I}) + \boldsymbol{g}_{\mathrm{out}}. \label{eq:additive95feature95fusion}\tag{33}\] where \(\boldsymbol{W}_r\in\mathbb{R}^{C\times 3N}\) for \(r\in\{1,2,3\}\). Equation 33 makes clear that the classifier combines reservoir blocks additively. It jointly optimizes their class-feature coefficients, but does not explicitly introduce bilinear cross-reservoir interactions.

The readout parameters are optimized by minimizing the empirical cross-entropy with \(L_2\) regularization: \[\mathcal{L} = - \frac{1}{S} \sum_{s=1}^{S} \ln \left[ \frac{ \exp\!\left(z_{y_s}(\boldsymbol{I}_s)\right) }{ \sum_{c=1}^{C} \exp\!\left(z_c(\boldsymbol{I}_s)\right) } \right] + \frac{\lambda}{2} \left\| \boldsymbol{W}_{\mathrm{out}} \right\|_F^2. \label{eq:loss95function}\tag{34}\] The bias vector \(\boldsymbol{g}_{\mathrm{out}}\) is not included in the \(L_2\) penalty. For all results in the main comparison, \(S=60000\), \(\lambda=10^{-3}\), the learning rate is \(0.05\), the batch size is \(128\), the number of epochs is \(150\), and the readout seed is \(20260621\).

Because the readout is optimized on the uncompressed direct sum of the reservoir features, all class-feature coefficients are trained jointly under the same global loss. This construction relaxes the reservoir-wise compression constraint of independent local readouts, while remaining an additive linear classifier.

3.5 Relation to Generative Thermodynamic Computing↩︎

The present discriminative moment-resolved discriminative framework is complementary to recent generative thermodynamic computing. In generative thermodynamic computing, the trained Langevin system is required to transform noise into structured data, and the learning objective is naturally formulated in terms of reverse-trajectory likelihood and thermodynamic irreversibility. By contrast, the present work uses the Langevin system as a finite-time physical feature extractor. The input is already structured, and the computational question is how much discriminative information can be extracted from the transient nonequilibrium response distribution.

This distinction leads to different computational objectives. Generative thermodynamic computing optimizes stochastic dynamics to transform noise into structured samples, whereas the present work uses finite-time trajectories to construct discriminative polynomial-moment features. The multi-reservoir results further indicate that reservoirs with distinct initialization and training histories can exhibit different classification-error overlaps. Thus, trajectory-level generative likelihood and moment-resolved finite-time readout represent complementary uses of Langevin dynamics for physical machine learning.

4 Conclusion↩︎

In summary, this work extends finite-time Langevin computing from mean-only readout to a moment-resolved representation based on the raw polynomial observables \(\mathbb{E}[\boldsymbol{x}]\), \(\mathbb{E}[\boldsymbol{x}^{\odot2}]\), and \(\mathbb{E}[\boldsymbol{x}^{\odot4}]\). These raw moments should not be interpreted as pure measurements of variance and kurtosis: they combine displacement and central-shape contributions while remaining naturally aligned with the polynomial structure of the onsite dynamics.

Under the fixed MNIST reproduction protocol, the complete three-reservoir feature representation gives the best observed accuracy of \(96.95\%\). The difference from the strongest single-reservoir model is \(0.13\) percentage points and is not statistically significant under the reported exact McNemar test. The wrong-set Jaccard analysis shows that the reservoirs exhibit different degrees of classification-error overlap, but it does not by itself establish feature-space decorrelation or collinearity. The results therefore provide suggestive, rather than definitive, evidence that heterogeneous response bases can contribute complementary information.

The current comparison also increases feature dimension, sampling resources, and readout capacity when additional reservoirs are introduced. Future work should include equal-resource controls, repeated trajectory and readout seeds, raw-versus-central moment comparisons, and observation-time or temperature scans. Such analyses will be needed to determine how much of the observed behavior originates specifically from dynamical heterogeneity, higher-order response statistics, and finite-time nonequilibrium evolution.

5 DATA AND CODE AVAILABILITY↩︎

The source code, fixed configuration files, random seeds, numerical protocols, and scripts used for the MNIST \(60000/10000\) reproduction are openly available at https://github.com/djz0924/Nonlinear_Thermodynamic_Computer. The accompanying release archive provides the \(K_1\), \(K_2\), and \(K_3\) reservoir weights, SHA256 checksums, final readout parameters, prediction files, software-environment metadata, and machine-readable result summaries. Large precomputed feature files are provided as optional release artifacts. The MNIST dataset is publicly available through its standard distribution channels and is not redistributed in the source repository.

References↩︎