Digitally Optimized Initializations for Fast Thermodynamic Computing

Abstract

Thermodynamic computing harnesses the relaxation dynamics of physical systems to perform matrix operations. A key limitation of such approaches is the often long thermalization time required for the system to approach equilibrium with sufficient accuracy. Here, we introduce a hybrid digital-thermodynamic algorithm that substantially accelerates relaxation through optimized initializations inspired by the Mpemba effect. In the proposed scheme, a classical digital processor efficiently computes an initialization that suppresses slow relaxation modes, after which the physical system performs the remaining computation through its intrinsic relaxation dynamics. We focus on overdamped Langevin dynamics for quadratic energy landscapes, analyzing the spectral structure of the associated Fokker-Planck operator and identifying the corresponding optimal initial covariances. This yields a predictable reduction in thermalization time, determined by the spectrum of the encoded matrix. We derive analytic expressions for the resulting speedups and numerically analyze thermodynamic implementations of matrix inversion and determinant computation as concrete examples. Our results show that optimized initialization protocols provide a simple and broadly applicable route to accelerating thermodynamic computations.

Figures (3)

• Figure 1: Hybrid digital-thermodynamic computing protocol with Mpemba-inspired initialization. (a) A classical processor computes an optimized initialization for the thermodynamic computer. Given an input matrix $\textbf{A}$, a Lanczos routine efficiently extracts the $K$ smallest eigenpairs, which are used to prepare an initial condition $\mathbf{x}_0^{\mathrm{opt}}$. This initialization is encoded in a thermodynamic device consisting of coupled harmonic oscillators (implemented as LC electrical circuits) that relax toward equilibrium, from which various matrix functions $f(\mathbf{A})$ can be estimated. (b) Comparison between standard and optimized initialization. Starting from the trivial initial state $\mathbf{x}_0^{\mathrm{std}}=\mathbf{0}$ generally excites all relaxation modes, whereas the optimized initialization prethermalizes the $K$ slowest modes, leading to a faster thermalization following a Mpemba-type speedup.
• Figure 2: Accelerating thermodynamic matrix inversion. (a) Ensemble with linearly spaced eigenvalues and Haar-random eigenvectors \ref{['eq:spec:linear:spacing', 'eq:haar:random:eigenv']}. (b) Positive Wishart random matrices \ref{['eq:wishart:def']}. Both panels show the time-dependent Frobenius distance between the evolving covariance matrix and the exact inverse \ref{['eq:cov:error']} for $d=500$. The parameter $K$ denotes the number of slow modes that are prethermalized in the optimized initialization \ref{['eq:optimized:covariance:matrix:initialization:k:modes']}. Insets show the resulting speedups for $K=10$ as a function of the matrix dimension $d$, evaluated at three different thermalization error thresholds $\epsilon_t$. Black lines indicate the corresponding asymptotic speedup \ref{['eq:mpemba:speedup']}. We set $\delta=0.5$, $m=1.1d$, $\mu=k_BT=1$, and average over $50$ random matrices for each panel. The shaded regions indicate the full min-max spread over the sampled ensemble
• Figure 3: Accelerating thermodynamic matrix-determinant computation. (a) Matrices with a linearly spaced spectrum and Haar-random eigenvectors, \ref{['eq:spec:linear:spacing', 'eq:haar:random:eigenv']}. (b) Random positive Wishart matrices, \ref{['eq:wishart:def']}. The main panels show the relative error with respect to the exact determinant \ref{['eq:rel:log:error']} as a function of the total compute time $T$\ref{['eq:total:compute:time']}, comparing standard initializations ($K=0$) with Mpemba-optimized initializations ($K=10$) for different numbers of trajectories. Insets show the resulting speedup factor $S_\mathrm{det}$ as a function of the matrix dimension $d$, where the total time required to reach the target error $\epsilon$ is obtained by log-log interpolation of the median error curve between the two bracketing $N_{\mathrm{traj}}$ values. For the thermalization time $t_0$\ref{['eq:thermalization:time:det']}, we set the threshold $\epsilon_\mathrm{t} = 10^{-4}$. We use $\tau =2$, $\delta = 0.5$, $m = 1.5d$, $\mu = k_BT = 1$, and $dt = 10^{-4}$. In the main panel, $d = 100$. Results are averaged over $50$ independent random matrix realizations per panel; error bars indicate the spread across realizations and markers denote the median.