Thermodynamic Linear Algebra

Maxwell Aifer, Kaelan Donatella, Max Hunter Gordon, Samuel Duffield, Thomas Ahle, Daniel Simpson, Gavin Crooks, Patrick J. Coles

Normal Computing Corporation, New York, New York, USA

Linear algebra is central to many algorithms in engineering, science, and machine learning; hence, accelerating it would have tremendous economic impact. Quantum computing has been proposed for this purpose, although the resource requirements are far beyond current technological capabilities. We consider an alternative physics-based computing paradigm based on classical thermodynamics, to provide a near-term approach to accelerating linear algebra. At first sight, thermodynamics and linear algebra seem to be unrelated fields. Here, we connect solving linear algebra problems to sampling from the thermodynamic equilibrium distribution of a system of coupled harmonic oscillators. We present simple thermodynamic algorithms for solving linear systems of equations, computing matrix inverses, and computing matrix determinants. Under reasonable assumptions, we rigorously establish asymptotic speedups for our algorithms, relative to digital methods, that scale linearly in matrix dimension. Our algorithms exploit thermodynamic principles like ergodicity, entropy, and equilibration, highlighting the deep connection between these two seemingly distinct fields, and opening up algebraic applications for thermodynamic computers.

§ INTRODUCTION

Basic linear algebra primitives like solving linear systems and inverting matrices are present in many modern algorithms. Such primitives are relevant to a multitude of applications, for example optimal control of dynamic systems and resource allocation. They are also a common subroutine of many artificial intelligence (AI) algorithms, and account for a substantial portion of the time and energy costs in some cases. The most common method to perform these primitives is LU decomposition, whose time-complexity scales as $O(d^3)$. Many proposals have been made to accelerate such primitives, for example using iterative methods such as the conjugate gradient method.
In the last decade, these primitives have been accelerated by hardware improvements, notably by graphical processing units (GPUs), fueling massive parallelization. However, the scaling of these methods is still a prohibitive factor, and obtaining a good approximate solution to a dense matrix of more than a few tens of thousand dimensions remains challenging.

Exploiting physics to solve mathematical problems is a deep idea, with much focus on solving optimization problems <cit.>. In the context of linear algebra, much attention has been paid to quantum computers <cit.>, since the mathematics of discrete-variable quantum mechanics matches that of linear algebra. A quantum algorithm <cit.> to solve linear systems has been proposed, which for sparse and well-conditioned matrices scales as $\log d$. However, the resource requirements <cit.> for this algorithm are far beyond current hardware capabilities. More generally building large-scale quantum hardware has remained difficult <cit.>, and variational quantum algorithms for linear algebra <cit.> have battled with vanishing gradient issues <cit.>.

Therefore, the search for alternative hardware proposals that can exploit physical dynamics to accelerate linear algebra primitives has been ongoing. Notably, memristor crossbar arrays have been of interest for accelerating matrix-vector multiplications <cit.>. Solving linear systems has also been the subject of analog computing approaches <cit.>.

Recently, we defined a new class of hardware, built from stochastic, analog building blocks, which is ultimately thermodynamic in nature <cit.>. (See also probabilistic-bit computers <cit.> and thermodynamic neural networks <cit.> for alternative approaches to thermodynamic computing <cit.>). AI applications like generative modeling are a natural fit for this thermodynamic hardware, where stochastic fluctuations are exploited to generate novel samples.

In this work, we surprisingly show that the same thermodynamic hardware from Ref. <cit.> can also be used to accelerate key primitives in linear algebra. Thermodynamics is not typically associated with linear algebra, and connecting these two fields is therefore non-trivial.
Here, we exploit the fact that the mathematics of harmonic oscillator systems is inherently affine (i.e., linear), and hence we can map linear algebraic primitives onto such systems. (See also Ref. <cit.> for a discussion of harmonic oscillators in the context of quantum computing speedups.) We show that simply by sampling from the thermal equilibrium distribution of coupled harmonic oscillators, one can solve a variety of linear algebra problems.

Specifically we develop thermodynamic algorithms for the following linear algebraic primitives: (i) solving a linear system $Ax = b$, (ii) estimating a matrix inverse $A^{-1}$, (iii) solving Lyapunov equations <cit.> of the form $A\Sigma+ \Sigma A^{\intercal} = \mathds{1}$ and (iv) estimating the determinant of a symmetric positive definite matrix $A$. We show that if implemented on thermodynamic hardware, these methods scale favorably with problem size compared to digital algorithms. Our numerical simulations corroborate our analytical scaling results and also provide evidence of the fast convergence of these primitives with the wall-clock time, with the speedup relative to digital methods getting more pronounced with increasing dimension and condition number.

We remark that there is a connection between our thermodynamic algorithms and digital Monte-Carlo (MC) algorithms that were developed for linear algebra <cit.>. Namely, our algorithms can be viewed as a continuous-time version of these digital MC algorithms. However, we emphasize that the continuous time (i.e., physics-based rather than physics-inspired) nature of our algorithms is crucial for obtaining our predicted asymptotic speedup. Additionally, thermodynamic algorithms can be run on a single device <cit.> whereas efficient digital MC linear algebra requires extensive parallelization <cit.>.

§ RESULTS

§.§ Algorithmic Scaling

In Table <ref>, we summarize the asymptotic scaling results for our thermodynamic algorithms as compared to the best state-of-the-art (SOTA) digital methods for dense symmetric positive-definite matrices. The derivations of these results can be found in the Supplemental Information, and are based on bounds obtained for physical thermodynamic quantities, including correlation times, equilibration times, and free energy differences. As one can see from Table <ref>, an asymptotic speedup is predicted for our thermodynamic algorithms relative to the digital SOTA algorithms. Specifically, a speedup that is linear in $d$ is expected for each of the linear algebraic primitives (ignoring a possible dependence of $\kappa$ on $d$). We remark that the complexity of analog algorithms is subtle <cit.> and depends, e.g., on assumptions of how the hardware size grows with problem size. The assumptions made to obtain our scaling results are detailed in the Methods section. In what follows, we systematically present our thermodynamic algorithms for various linear algebraic primitives.

Problem     Digital SOTA

Linear System    $O(\min\{d^\omega,d^2 \sqrt{\kappa}\})$     $O( d \kappa^2\varepsilon^{-2})$

Matrix Inverse      $O(d^\omega)$      $O( d^2 \kappa \varepsilon^{-2})$

Lyapunov Equation    $O(d^3)$       $O( d^2 \kappa \varepsilon^{-2})$

Matrix Determinant     $O(d^\omega)$       $O(d \kappa  \ln(\kappa)^3 \varepsilon^{-2})$

Comparison of asymptotic complexities of linear algebra algorithms. Here, $d$ is the matrix dimension, $\kappa$ is the condition number, and $\varepsilon$ is the error. For our thermodynamic algorithms, the complexity depends on the dynamical regime. Here we display the overdamped dynamics which have marginally better complexity than the underdamped equivalents. For the digital SOTA, the complexity of solving symmetric, positive definite linear systems, matrix inverse, Lyapunov equation, and matrix determinant problems are respectively for algorithms based on: conjugate gradient method <cit.>, fast matrix multiplication/inverse <cit.>, Bartels-Stewart algorithm <cit.>, and LUP decomposition <cit.>. $\omega\approx 2.3$ denotes the matrix multiplication constant.

§.§ Solving Linear Systems of Equations

The celebrated linear systems problem is to find $x \in \mathbb{R}^d$ such that

\begin{equation}
\label{eq:linear-system-def}
A x = b,
\end{equation}

given some invertible matrix $A \in \mathbb{R}^{d \times d}$ and nonzero $b \in \mathbb{R}^d$. We may assume without loss of generality that the matrix $A$ in Eq. (<ref>) is symmetric and positive definite (SPD); if $A$ is not SPD, then we may consider the system $A^\intercal A x = A^\intercal b$, whose solution $x=A^{-1} b$ is also the solution of $A x = b$. Note that this will affect the total runtime, but still allows for asymptotic scaling improvements with respect to digital methods, in some cases[Constructing an SPD system from a generic one in this way results in the squaring of the condition number, which influences performance.]. In what follows, we will therefore assume that $A$ is SPD.

Diagram of our thermodynamic algorithm for solving linear systems and inverse estimation. The system of linear equations, or the matrix $A$, is encoded into the thermodynamic hardware, the system is then allowed to evolve until the stationary distribution has been reached, when the trajectory is then integrated to estimate the sample mean or covariance. This gives estimates of the solution of the linear system or the inverse of $A$ respectively.

Now let us connect this problem to thermodynamics. We consider a macroscopic device with $d$ degrees of freedom, described by classical physics. Suppose the device has potential energy function:

\begin{equation}
\label{eq:device-potential}
U(x) = \frac{1}{2}x^\intercal A x - b^\intercal x,
\end{equation}

where $A \in \text{SPD}_{d}(\mathbb{R})$. Note that this is a quadratic potential that can be physically realized with a system of harmonic oscillators, where the coupling between the oscillators is determined by the matrix $A$, and the $b$ vector describes a constant force on each individual oscillator. (We remark that while Figure <ref> depicts mechanical oscillators, from a practical perspective, one can build the device from electrical oscillators such as RLC circuits.)

Suppose that we allow this device to come to thermal equilibrium with its environment, whose inverse temperature is $\beta = 1/k_B T$. At thermal equilibrium, the Boltzmann distribution describes the probability for the oscillators to have a given spatial coordinate: $f(x) \propto \exp (- \beta U(x))$. Because $U(x)$ is a quadratic form, $f(x)$ corresponds to a multivariate Gaussian distribution. Thus at thermal equilibrium, the spatial coordinate $x$ is a Gaussian random variable

\begin{equation}
\label{eq:equilibrium-dist}
x\sim \mathcal{N}[A^{-1}b, \beta^{-1} A^{-1}].
\end{equation}

The key observation is that the unique minimum of $U(x)$ occurs where $Ax-b = 0$, which also corresponds to the unique maximum of $f(x)$. For a Gaussian distribution, the maximum of $f(x)$ is also the first moment $\braket{x}$. Thus, we have that, at thermal equilibrium, the first moment is the solution to the linear system of equations:
\begin{equation}
\label{eq:firstmoment}
\braket{x} = A^{-1}b.
\end{equation}

From this analysis, we can construct a simple thermodynamic protocol for solving linear systems, which is depicted in Figure~\ref{fig:LinearSystems}. Namely, the protocol involves realizing the potential in Eq.~\eqref{eq:device-potential}, waiting for the system to come to equilibrium, and then sampling $x$ to estimate the mean $⟨x|$⟩ of the distribution. This mean can be approximated using a time-average, defined as

\begin{equation}
\label{eq:x-bar-def}
\braket{x}\approx \bar{x}(\tau) = \frac{1}{\tau} \int_{t_0}^{t_0+\tau} dt' x(t'),
\end{equation}

where $t_0$ must be sufficiently large to allow for equilibration and $\tau$ must be sufficiently large for the average to converge to a desired degree of precision. The eventual convergence of this time average to the mean is the content of the ergodic hypothesis <cit.>, which is often assumed for quite generic thermodynamic systems. It should be mentioned that the mean could also be approximated as the average of a sequence of samples; however the integration approach has the advantage that it can conveniently be implemented in a completely analog way (for example, using an integrator electrical circuit), which obviates the need for transferring data from the physical device until the end of the protocol.

Figure <ref> shows the equilibration process for both a single trajectory (left panel) and the overall distribution (right panel). One can see the ergodic principle illustrated in this figure, since the time dynamics of a single trajectory at thermal equilibrium are representative of the overall ensemble.

Equilibration of the thermodynamic system. The process of equilibration is depicted on the single-trajectory level (left) and on the distribution level (right). The trajectory dynamics are described by the overdamped Langevin equation and the distributional dynamics by the Fokker-Planck equation <cit.> The system displays ergodicity, as the time average of a single trajectory (blue curve, left) approaches the ensemble average (dots, right) in the long-time limit. Time and the coordinate vector $(x_1,x_2)$ are in arbitrary units.

The overall protocol can be summarized as follows.

[title=Linear System Protocol, breakable]

* Given a linear system $A x = b$, set the potential of the device to

\begin{equation}
U(x) = \frac{1}{2} x^\intercal A x - b^\intercal x
\end{equation}

at time $t=0$.

* Choose equilibration tolerance parameters $\varepsilon_{\mu0}, \varepsilon_{\Sigma 0} \in \mathbb{R}^+$, and choose the equilibration time

\begin{equation}
\label{eq:UDL-eq-time}
t_0 \geq \widehat{t}_0,
\end{equation}

where $\widehat{t}_0$ is computed from the system's physical properties or using heuristic methods based on Eqs. (<ref>), (<ref>). Allow the system to evolve under its dynamics until $t=t_0$, which ensures that $\left\|\braket{x} - A^{-1} b\right\|/\|A^{-1} b\| \leq \varepsilon_{\mu 0}$ and $\left\|\Sigma - \beta^{-1} A^{-1} \right\|/\|\beta^{-1}A^{-1}\| \leq \varepsilon_{\Sigma0}$.

* Choose error tolerance parameter $\varepsilon_x$ and success probability $P_\varepsilon$, and choose the integration time

\begin{equation}
\tau \geq \widehat{\tau},
\end{equation}

where $\widehat{\tau}$ is computed from the system's physical properties, Eq. (<ref>) or (<ref>). Use an analog integrator to measure the time average

\begin{equation}
\bar{x} = \frac{1}{\tau} \int_{t_0}^{t_0 + \tau} dt\, x(t),
\end{equation}

which satisfies $\left\|A \bar{x} - b\right\|/\|b\| \leq \varepsilon_x$ with probability at least $P_\delta$.

In order to implement the protocol above, the necessary values of $\widehat{t}_0$ and $\widehat{\tau}$ must be identified, which requires a more quantitative description of equilibration and ergodicity. To obtain such a description, a model of the system's microscopic dynamics may be introduced. Given that the system under consideration is composed of harmonic oscillators in contact with a heat bath, it is natural to allow for damping (i.e., energy loss to the bath) and stochastic thermal noise, which always accompanies damping due to the fluctuation-dissipation theorem <cit.>. The Langevin equation accounts for these effects, and specifically we consider two common formulations, the overdamped Langevin (ODL) equation  and the underdamped Langevin (UDL) equations. In the Methods section, we provide additional details on ODL and UDL dynamics, and we provide explicit formulas for $\widehat{t}_0$ and $\widehat{\tau}$ for the overdamped and underdamped regimes.

§.§ Estimating the Inverse of a Matrix

The results of the previous section rely on estimating the mean of $x$, but make no use of the fluctuations in $x$ at equilibrium. By using the second moments of the equilibrium distribution, we can go beyond solving linear systems. For example it is possible to find the inverse of an SPD matrix $A$. As mentioned, the stationary distribution of $x$ is $ \mathcal{N}[A^{-1} b, \beta^{-1} A^{-1}]$, meaning the inverse of $A$ can be obtained by evaluating the covariance matrix of $x$. This can be accomplished in an entirely analog way, using a combination of analog multipliers and integrators. By setting $b=0$ for this protocol, we ensure that $\braket{x} = 0$, so the stationary covariance matrix is, by definition

\begin{equation}
\Sigma_\text{s} = \lim_{t \to \infty}\braket{x(t) x^\intercal(t)}.
\end{equation}

In order to estimate this, we again perform time averages after allowing the system to come to equilibrium

\begin{equation}
\Sigma_\text{s} \approx \overline{x x^\intercal} = \frac{1}{\tau}
\int_{t_0}^{t_0 + \tau} dt\,  x(t) x^\intercal(t).
\end{equation}

It is therefore necessary to have an analog component which evaluates the product $x_i(t) x_j(t)$ for each pair $(i,j)$, resulting in $d^2$ analog multiplier components. Each of these products is then fed into an analog integrator component, which computes one element of the time-averaged covariance matrix

\begin{equation}
\Sigma_{\text{s},ij} \approx \frac{1}{\tau} \int_{t_0}^{t_0 + \tau}dt \, x_i(t) x_j(t).
\end{equation}

While the equilibration time is the same as for the linear system protocol, the integration time is different, because in general the covariance matrix is slower to converge than the mean. We now give a detailed description of the inverse estimation protocol, assuming ODL dynamics (the corresponding results for underdamped dynamics can be found in the Supplemental Information). In the Methods section,  we provide explicit formulas for $\widehat{t}_0$ and $\widehat{\tau}$ for the Inverse Estimation Protocol. We remark that our matrix inversion algorithm is a special case of our general algorithm for solving Lyapunov equations; the latter is presented in the Supplemental Information.
[title=Inverse Estimation Protocol,breakable]

* Given a positive definite matrix $A$, set the potential of the device to

\begin{equation}
U(x) = \frac{1}{2} x^\intercal A x
\end{equation}

at time $t=0$.

* Choose equilibration tolerance parameter $\varepsilon_{\Sigma 0} \in \mathbb{R}^+$, and choose the equilibration time

\begin{equation}
\label{eq:UDL-eq-time}
t_0 \geq \widehat{t}_0,
\end{equation}

where $\widehat{t}_0$ is computed from the system's physical properties, Eq. (<ref>) or (<ref>). Allow the system to evolve under its dynamics until $t=t_0$, which ensures that $\left\|\Sigma - \beta^{-1} A^{-1} b\right\|/\|\beta^{-1}A^{-1}\| \leq \varepsilon_\Sigma$.

* Choose error tolerance parameter $\varepsilon_\Sigma$ and success probability $P_\varepsilon$, and choose the integration time

\begin{equation}
\tau \geq \widehat{\tau},
\end{equation}

where $\widehat{\tau}$ is computed from the system's physical properties, Eq. (<ref>) or (<ref>). Use analog multipliers and integrators to measure the the time averages

\begin{equation}
\overline{x_i x_j} = \frac{1}{\tau^2} \int_{t_0}^\tau dt\, x_i(t) x_j(t),
\end{equation}

which satisfies $\|\overline{x x^\intercal} - \beta^{-1}A^{-1}\|_F/\|\beta^{-1}A^{-1}\|_F \leq \varepsilon_A$ with probability at least $P_\varepsilon$.

§.§ Estimating the Determinant of a Matrix

The determinant of the covariance matrix appears in the normalization factor of a multivariate normal distribution, whose density function is

\begin{equation}
f_{\mu;\Sigma}(x) = (2 \pi)^{-d/2}\left| \Sigma\right|^{-1/2}
\exp \left(-\frac{1}{2} x^\intercal \Sigma^{-1} x \right),
\end{equation}

and it is therefore natural to wonder whether hardware which is capable of preparing a Gaussian distribution may be used to somehow estimate the determinant of a matrix. This can in fact be done, as the problem is equivalent to the estimation of free energy differences, an important application of stochastic thermodynamics. Recall that the difference in free energy between equilibrium states of potentials $U_1$ and $U_2$ is <cit.>

\begin{equation}
\Delta F = F_2 - F_1 = -\beta^{-1} \ln \left(\frac{\int dx \, e^{-\beta U_2(x)}}{\int dx \, e^{-\beta U_1(x)}}\right).
\end{equation}

Suppose the potentials are quadratic, with $U_1(x) = x^\intercal A_1 x$ and $U_2(x) = x^\intercal A_2 x$. Then each integral simplifies to the inverse of a Gaussian normalization factor,

\begin{equation}
\int dx\,  e^{-\beta V_j(x)} = (2\pi)^{d/2}\sqrt{\beta^{-1}\left|A_j^{-1} \right|},
\end{equation}

\begin{equation}
\Delta F = -\beta^{-1} \ln \left(\sqrt{\frac{\left| A_2^{-1}\right|}{\left| A_1^{-1}\right|}}\right) = -\beta^{-1} \ln \left(\sqrt{\frac{\left| A_1\right|}{\left| A_2\right|}}\right).
\end{equation}

This suggests that the determinant of a matrix $A_1$ can found by comparing the free energies of the equilibrium states with potentials $U_1$ and $U_2$ (where $A_2$ has known determinant), and then computing

\begin{equation}
\left| A_1\right| = e^{-2 \beta \Delta F} \left| A_2\right|.
\end{equation}

Fortunately, the free energy difference $\Delta F$ can be found, assuming we have the ability to measure the work which is done on the system as the potential $U(x)$ is changed from $U_1$ to $U_2$. According to the Jarzynski equality <cit.>, the free energy difference between the (equilibrium) states in the initial and final potential is

\begin{equation}
e^{-\beta \Delta F} = \braket{e^{-\beta W}},
\end{equation}

where $\braket{\cdot}$ denotes an average over all possible trajectories of the system between time $t=0$ and time $t=τ$, weighed by their respective probabilities. This may be approximated by an average over $N$ repeated trials,
\begin{equation}
e^{-\beta \Delta F}\approx \overline{e^{-\beta W}}  \equiv \frac{1}{N}\sum_{j=1}^N e^{-\beta W_j}.
\end{equation}
However, while Jarzynski's relation may be applied directly to estimate the free energy difference, this estimator has large bias and is slow to converge. Far more well-behaved estimators have been found based on work measurements. For simplicity, we here provide the expression based on Jarzynski's estimator, while in the Methods section and the Supplemental Information we refer to more suitable estimators. In summary, the determinant of $A_1$ is approximated by
\begin{equation}
\left| A_1\right| \approx \left(\overline{e^{-\beta W}}\right)^2 \left| A_2\right|.
\end{equation}
In practice we will generally be interested in the log determinant to avoid computational overflow. This is
\begin{equation}
\ln\left(\left| A_1\right|\right) \approx 2 \ln\left(\overline{e^{-\beta W}}\right)  + \ln\left( \left| A_2\right|\right).
\end{equation}
It is shown in the Supplemental Information that to estimate the log determinant to within (absolute) error $δ_LD$ with probability at least $P_δ$, the total amount of time required is roughly
\begin{align}
\tau \approx \frac{ d\,\ln(\kappa)^2}{\delta^2_\text{LD}(1-P_\delta)} \ln\left(\chi^2 \kappa^{3/2}\varepsilon_{\Sigma0}^{-1}\left[\frac{1}{4\zeta_\text{max}^{2}}+1\right]\right)
\tau_\text{r(UD)} = O(d \, \ln(\kappa)^3).
\end{align}
We also present numerical simulations of a protocol for determinant estimation that does not include directly measuring the work in the Supplemental Information.

\subsection{Convergence and comparison to digital algorithms}

\subsubsection{Convergence}

\begin{figure}[t]
\centering
\includegraphics[width=0.98\linewidth]{figures/TLA_times.pdf}
\caption{\textbf{Error of our thermodynamic algorithms as a function of the analog integration time for different dimensions.} Matrices $A$ are drawn from a Wishart distribution with $2d$ degrees of freedom. Vertical dashed lines are the times $t_C$ at which error goes below a threshold (horizontal dashed line). Inset: Crossing time $t_C$ as a function of dimension $d$. (A) For the linear systems algorithm, a linear relationship between dimension and the analog dynamics runtime is observed.
(B) For the matrix inversion algorithm, a quadratic relationship between dimension and the analog dynamics runtime is observed.}
\label{fig:ErrorVsAnalogTime}
\end{figure}

We now present several numerical experiments to corroborate our analytical results. Figure~\ref{fig:ErrorVsAnalogTime}(A) displays the convergence of the absolute error, $||x̅ - A^-1b||$ where $||.||$ denotes the 2-norm, as a function of time for our thermodynamic linear systems algorithm. This plot shows that the expected convergence time to reach a given error is linearly proportional to the dimension of the system, which is in agreement with the analytical bounds that we presented above.

Similarly, let us examine the performance of the inverse estimation protocol. We employ the absolute error on the inverse, $ Ã^-1 - A^-1_F$ where $·_F$ denotes the Frobenius norm. Figure~\ref{fig:ErrorVsAnalogTime}(B) shows the convergence of the error as a function of the analog dynamics time for our thermodynamic inverse estimation algorithm. We see that the expected convergence time to reach a given error is quadratic ($∝d^2$) in the dimension, in agreement with the analytical bounds presented above.

\begin{figure}[t]
\centering
\includegraphics[width=\linewidth]{figures/error_vs_time_linearsys.pdf}
\caption{\textbf{Comparison of the error $||\bar{x} - A^{-1}b||$ of the thermodynamic algorithm (TA) to solve linear systems with the conjugate gradient method and Cholesky decomposition as a function of total runtime.} Panels (a)-(c): the TA is shown for different values of $k_BT$ (units of 1/$\gamma$) for each dimension in $\{100, 1000, 5000\}$. Random matrices are drawn from the Wishart distribution and then mixed with the identity such that their condition numbers are respectively 120, 1189, 5995. Panels (d)-(f): same quantities with a fixed condition number $\kappa$, respectively 199, 1190, and 7880 for fixed dimension $d = 1000$. Calculations were performed on an Nvidia RTX 6000 GPU.}
\label{fig:linearsys_comparison_d}
\end{figure}

\subsubsection{Comparison to digital algorithms}

Another question of key importance is how the thermodynamic algorithm is expected to perform in practical scenarios, i.e., when being run on real thermodynamic hardware. Due to the hardware being analog in nature, this involves additional digital-to-analog compilation steps. To investigate this question, we consider a timing model for the thermodynamic algorithm, based on the hardware proposal described  Ref.~\cite{coles2023thermodynamic} (See the Supplemental Information for a brief summary of this hardware, whose dynamics correspond to the overdamped regime as in Eq.~\eqref{eq:ODLdynamics_main}). This model includes all the digital, digital-to-analog and analog operations needed to solve the problem, starting with a matrix $A$ stored on a digital device, and sending back the solution $x$ from the thermodynamic system to the digital device. Note that this includes a compilation step that scales as $O(d^2)$, which is absent for the digital methods\footnote{Cholesky and conjugate gradients are run on a digital computer, and the initial matrix is stored on that same computer, hence there is no transfer cost, unlike for the thermodynamic algorithm.}. Assumptions about this model are detailed in the Methods section. Note that analog imprecision is not taken into account in these experiments, and is the subject of further investigations~\cite{aifer2024error}.

Figure~\ref{fig:linearsys_comparison_d} plots the absolute error for solving linear systems as a function of time for the thermodynamic algorithm (TA), the conjugate gradient (CG) method, and the Cholesky decomposition (which is exact). In panels (a) - (c) we explore how the methods converge with varying $κ$ and $d$. While at low dimensions our method performs poorly with respect to the Cholesky decomposition and only slightly better than CG, it  becomes very competitive for dimensions $d = 1000$ and $d=5000$. Panels (d) - (f) show the error as a function of time for different condition numbers, at fixed dimension. One can see that as $κ$ grows (as conditioning is worse) our method becomes more competitive with CG. This suggests that, even in practical scenarios where we account for realistic computational overhead issues, our thermodynamic linear systems algorithm can outperform SOTA digital methods, especially for large $d$ and large $κ$.

Figure~\ref{fig:linearsys_comparison_d} also shows that the thermodynamic algorithm performs significantly better than the CG method at early times, although the CG method ultimately achieves a higher quality result for later times. This suggests that the thermodynamic algorithm is ideally suited to providing an approximate solution in a short amount of time. Nevertheless, we note that the effective temperature of the thermodynamic hardware is an important parameter, and one can lower this temperature to achieve higher precision solutions from the thermodynamic hardware, as can be seen from the curves in Fig.~\ref{fig:linearsys_comparison_d}.

Using a timing model similar to that employed for the linear systems protocol, we performed a runtime comparison to Cholesky decomposition for the task of matrix inversion. The results are shown in Fig.~\ref{fig:inverse_estim_comparison_d}, where the error is plotted as a function of physical time for dimensions $100, 1000$, and $5000$. The dashed lines represent the corresponding times for Cholesky decomposition, for given dimensions. We see that as the dimension grows, the advantage with respect to the Cholesky decomposition also grows, thus highlighting a practical thermodynamic advantage. Our method for the inverse estimation therefore has the advantage of having well-defined convergence properties as a function of dimension and condition number (compared to other approximate methods for inverting dense matrices, which do not have well defined convergence properties), as well as leading to reasonable error values in practical settings.

Overall, these numerical experiments highlight the potential utility of thermodynamic hardware by showing the opportunity for speedup over SOTA digital methods, based on a simulated timing model of the thermodynamic device.

\begin{figure}[t]
\centering
\includegraphics[width=0.6\linewidth]{figures/error_vs_time_inverse.pdf}
\caption{\textbf{Comparison of the error of the thermodynamic algorithm (TA) to invert matrices with the Cholesky decomposition as a function of total runtime.} Dimensions $d = 100, 1000, 5000$, respectively in light green, light blue, and purple, are shown for the thermodynamic algorithm (solid lines) and the Cholesky decomposition (dashed lines). Here the condition numbers are respectively $\{120, 1189, 5995\}$. Calculations were performed on an Nvidia RTX A600 GPU.}
\label{fig:inverse_estim_comparison_d}
\end{figure}

\section{Discussion}

Various types of physics-based computers have been devised, which are supposed to expedite calculations by using physical processes to evaluate expensive functions \cite{haensch2018next, small1993general,feynman1982simulating,preskill2018quantum, nielsen2000quantum}. These devices (which include quantum computers and a number of distinct analog architectures) have been shown to offer theoretical advantages for solving certain problems, including linear systems of equations, but they have not found common use commercially. A key obstacle to harnessing the power of physical computing is that fluctuations in the system's state tend to cause errors that compound over time, and which cannot be corrected in a straightforward way \cite{huang2017analog} (as can be done for digital computers).

For this reason, we have considered \emph{thermodynamic} algorithms, which treat the naturally-present fluctuations as a resource, or at the very least are indifferent to them. In fact, we have introduced three distinct classes of thermodynamic algorithms: first-moment based, second-moment based, and all-moment based algorithms. Other thermodynamic algorithms will likely be discovered making use of third and higher moments, implying that such methods form a hierarchy. In some sense, using higher moments allows us to solve ``harder'' problems, for example inverting a matrix (which uses the second moments) is harder than solving a linear system of equations (which uses the first moment). Whether a precise relationship can be found between computational hardness and the hierarchy of thermodynamic algorithms is currently an open question.

Another open question concerns the optimality of these new thermodynamic algorithms. Our analysis implies that, while the time and energy costs of linear-algebraic primitives are negotiable, the product of time and energy necessary for a computation is fundamentally constrained (see Methods). It is therefore of interest to search for thermodynamic algorithms which achieve lower values of the energy-time product for these computations, and also to see whether such constraints may apply to other problems as well. We anticipate that non-equilibrium thermodynamics will be a crucial tool in exploring such resource tradeoffs for computation. For example, we have used the fact that a thermodynamic distance may be defined between equilibrium configurations of a system, and this distance determines the minimal amount of dissipated energy necessary to transition from one configuration to another in a finite time \cite{crooks2007measuring, cafaro2022thermodynamic, quevedo2007geometrothermodynamics,andresen1996finite,chen2021extrapolating}; the shorter the time of transition, the more dissipation must occur. Perhaps, then, the search for algorithms which have minimal energy-time product may be framed as a variational problem of minimizing length on the thermodynamic manifold. Although proofs of optimal algorithmic performance are notoriously hard to find in digital computing paradigms \cite{saptharishi2015survey}, unavoidable resource tradeoffs are  relatively mundane in thermodynamic analyses \cite{chiribella2022nonequilibrium, riechers2018transforming}, suggesting that computational cost may be fruitfully studied within the thermodynamic computing paradigm.

Aside from the theoretical questions mentioned, clearly the task of actually implementing our algorithms remains an important one. Recently, we created an electrical thermodynamic computing device~\cite{melanson2023thermodynamic} on a printed circuit board and used it to demonstrate the thermodynamic matrix inversion algorithm presented here, inverting 8 by 8 matrices using 8 coupled electrical oscillators. A potential next step would be to experimentally verify our predicted scaling of integration time with dimension (e.g., linear scaling for linear systems and quadratic scaling for matrix inversion), thus confirming our predicted speedup over digital methods. We anticipate that other researchers may independently seek to verify our results experimentally, leading to a rapid development of thermodynamic hardware. As a result, we predict that these methods will become appealing alternatives to digital algorithms, particularly in settings where it is desirable to trade some accuracy for better time and energy scaling.

In addition to our work's direct impact, the broader impact is laying the theoretical, mathematical foundations for the emerging paradigm of thermodynamic computing~\cite{conte2019thermodynamic}. Our work provides the first mathematical analysis, as well as the first numerical benchmarks, of potential speedups for thermodynamic hardware. Thus we have taken the somewhat vague notion of thermodynamic computing and made it concrete and precise, with a clear set of applications. Moving forward, we expect new applications to be discovered, beyond linear algebra, since one can simply modify the potential energy function $U(x)$ to solve, e.g., non-linear algebraic problems. There is also the exciting prospect of running multiple applications, such as the linear algebra ones here and the probabilistic AI ones discussed in Ref.~\cite{coles2023thermodynamic}, on the same thermodynamic hardware, providing the user with a flexible programming experience. One can envision that much of the amazing technological developments (compilers, simulators, programming languages, etc.) that have happened in quantum computing will likely happen for thermodynamic computing in the near future.

\section{Methods}

\subsection{Timing parameters for overdamped and underdamped regimes}

\subsubsection{Linear Systems}

For our linear systems algorithm, the overdamped Langevin (ODL) equation takes the form:
\begin{equation}
\label{eq:ODLdynamics_main}
dx = - \frac{1}{\gamma}(A x - b) dt +\mathcal{N}\left[0,2\gamma^{-1}\beta^{-1}\, dt\right],
\end{equation}
where $γ>0$ is called the damping constant and $β= 1/k_B T$ is the inverse temperature of the environment. The system has a physical timescale (which is clear from dimensional analysis) that we call the \emph{relaxation time} $τ_r = γ/A$. The condition number of $A$ is $κ= α_max/α_min$, where $α_1…α_d$ are the eigenvalues of $A$. With these definitions, we arrive (see Supplemental Information for derivation) at the following formulas for $t_0$ and $τ$ in the overdamped case, which can be used in the linear systems algorithm:
\begin{align}
\label{eq:t0-tau-linear-systems-odl}
\widehat{t}_0 =
\max \left\{\kappa \tau_\text{r}\ln \left(\kappa \varepsilon^{-1}_{\mu 0}\right),\frac{1}{2} \kappa \tau_\text{r}\ln\left(2\kappa\varepsilon_{\Sigma 0}^{-1}\right)\right\},\quad
\widehat{\tau} = \frac{2 \kappa^2 d\, \|A\|  }{\beta \|b\|^2 \varepsilon_x^2 (1-P_\varepsilon)}\tau_\text{r}.
\end{align}

The underdamped model is instead described by the UDL equations,
\begin{equation}
\label{eq:UDL-x}
dx  = \frac{1}{M} p \, dt, \qquad  dp = -(A x - b) \, dt- \frac{\gamma}{M} p \, dt + \mathcal{N}[0,2 \gamma \beta^{-1}  \mathbb{I} dt].
\end{equation}
We define $ξ= γ/2M$, $ω_j = √(α_j/M)$, and $ζ_j = ξ/ω_j$. Moreover, a timescale $τ_r(UD)$ can be identified for the underdamped system which is analogous to the quantity $τ_r$ associated with the overdamped system. In particular, we define $τ_r(UD) = ξ^-1$. We introduce a dimensionless quantity $χ$ as well, which is $χ= (1+ ξ/ω_min)^1/2(1- ξ/ω_min)^-1/2$.
With these definitions, we arrive (see Supplemental Information for derivation) at the following formulas for the timing parameters in the underdamped case:
\begin{align}
\label{eq:t0-tau-linear-systems-udl}
\widehat{t}_0 =  \max \left\{ \tau_{\text{r(UD)}} \ln \left( \kappa^{1/2}\chi \varepsilon_{\mu 0}^{-1}\right),\frac{1}{2} \tau_\text{r(UD)}\ln\left(\chi^2 \kappa^{3/2}\varepsilon_{\Sigma0}^{-1}\left[\frac{1}{4\zeta_\text{max}^{2}}+1\right]\right)\right\},\quad
\widehat{\tau} = \frac{2\sqrt{\kappa} \chi d\|A\| }{ \beta \|b\|^2\varepsilon_x^2(1-P_\varepsilon)} \tau_{\text{r(UD)}}.
\end{align}

An important distinction between the ODL and UDL regimes is that the random variable $x$ undergoes a Markov stochastic process in the ODL case, but is non-Markovian in the UDL case \cite{doerries2021correlation,h1989colored}. The simple interpretation of this non-Markovianity is that the underdamped system exhibits inertia, which is a form of memory-dependence. This inertia has a non-trivial (and sometimes beneficial) effect on the algorithm's overall performance, which is apparent from the scaling results in Table~\ref{tab:runtime_scaling}.

\subsubsection{Matrix Inversion}

The timing parameters for the inverse estimation protocol (as derived in the Supplemental Information) are, for the overdamped case,
\begin{align}
\label{eq:t0-tau-inverse-odl}
\widehat{t}_0 =
\frac{1}{2} \kappa \tau_\text{r}\ln\left(2\kappa\varepsilon_{\Sigma 0}^{-1}\right),\quad
\widehat{\tau} = \frac{4 \kappa d(d+1)}{(1-P_\varepsilon)\varepsilon_\Sigma^2}\tau_\text{r},
\end{align}
and for the underdamped case
\begin{align}
\label{eq:t0-tau-inverse-udl}
\widehat{t}_0 = \frac{1}{2} \tau_\text{r(UD)}\ln\left(\chi^2 \kappa^{3/2}\varepsilon_{\Sigma0}^{-1}\left[\frac{1}{4\zeta_\text{max}^{2}}+1\right]\right),\quad
\widehat{\tau} =\frac{4 \kappa d(d+1)}{(1-P_\varepsilon)\varepsilon_\Sigma^2}\tau_\text{r(UD)}.
\end{align}

\subsection{Energy-Time Tradeoff}

Note the appearance of the ratio $A/b^2$ in the time required to solve a linear system given by Eqs.~\eqref{eq:t0-tau-linear-systems-odl} and \eqref{eq:t0-tau-linear-systems-udl}. It is tempting to imagine that one might solve the system faster simply by multiplying $b$ by some constant $c$. Then, the time required to solve the system is apparently reduced by a factor of $c^2$, and the solution to original problem is obtained (up to a factor of $c$). A similar approach would be to multiply $A$ by a small number; of course, in practice it is not possible to solve linear systems of equations in vanishingly short periods of time, which is reflected in an energy-time tradeoff. As we explain in the Supplemental Information, there is an energy cost associated with initializing the system proportional to $b^⊺A^-1 b$, and this results in a re-formulation of Eqs.~\eqref{eq:t0-tau-linear-systems-odl} and \eqref{eq:t0-tau-linear-systems-udl} as lower bounds on the product of energy and time. If $ℰ$ is the energy required to solve the system $A x = b$, and $τ$ is the necessary time, then we have, in the overdamped case:
\begin{align}
\label{eq:time-energy-odl}
\mathcal{E}\tau &\geq \frac{2 \kappa^2 d  }{ \varepsilon_x^2 (1-P_\varepsilon)} \beta^{-1}\tau_\text{r}
\end{align}
and in the underdamped case:
\begin{align}
\label{time-energy-udl}
\mathcal{E}\tau \geq  \frac{2\sqrt{\kappa} \chi  d}{  \varepsilon_x^2(1-P_\varepsilon)} \beta^{-1}\tau_{\text{r(UD)}}.
\end{align}
This fundamental energy-time tradeoff appears naturally within this computational model. While digital computations can often be accelerated by investing more energy (for example, via parallelization), it is generally less obvious what exact form the relationship between time and energy cost takes, suggesting that thermodynamic algorithms may offer a new and useful perspective on algorithmic complexity.

\subsection{Detailed Algorithmic Scaling}

In the main text, we presented a simplified version of the detailed table shown in Table~\ref{tab:runtime_scaling_methods}. Table~\ref{tab:runtime_scaling_methods} breaks down the scaling into the overdamped and underdamped regimes, whereas Table~\ref{tab:runtime_scaling} just takes the best scalings of our algorithms. The latter is essentially the scaling associated with the underdamped regime. However, in practice, there can be some engineering advantages to working in the overdamped regime, and hence it is useful to see the complexities of both regimes.
\begin{table}[t]
\centering
\renewcommand{\arraystretch}{2}
\begin{tabular}{c|c|c|c}
\textbf{Problem} & \textbf{Digital SOTA}
& \textbf{This work (Overdamped)} & \textbf{This work (Underdamped)} \\
\hline
Linear System &$O(\min\{d^\omega,d^2 \sqrt{\kappa}\})$& $ O(d \kappa^2 \varepsilon^{-2})$ & $O(\gamma d \kappa^{3/2}  \varepsilon^{-2})$ \\
Matrix Inverse &  $O(d^\omega)$ & $O(d^2 \kappa \varepsilon^{-2})$ & $ O(\gamma  d^2 \kappa^2 \varepsilon^{-2})$\\
Lyapunov Equation &$O(d^3)$ & $O( d^2 \kappa \varepsilon^{-2})$&  $ O(\gamma d^2 \kappa^2\varepsilon^{-2})$ \\
Matrix Determinant &  $O(d^\omega)$& $O(d  \kappa \ln(\kappa)^3 \varepsilon^{-2})$ & $ O(\gamma d \kappa  \ln(\kappa)^3 \varepsilon^{-2})$\\
\end{tabular}
\caption{\textbf{Asymptotic complexities of linear algebra algorithms, including overdamped and underdamped regimes.} For our thermodynamic algorithms, the complexity depends on the dynamical regime, i.e., whether the dynamics are overdamped and underdamped, as shown in this table.}
\label{tab:runtime_scaling_methods}
\end{table}

\subsection{Assumptions}

A number of assumptions are made in the analytical derivations of the findings presented in the Results section. Certain aspects of the problem have been idealized in order to reveal the fundamental performance characteristics of the thermodynamic algorithms. The main assumptions are the following
\begin{itemize}
\item The dynamics of the system may be described by the ODL equation or the UDL equations.
\item The potential function $U(x)$, and in particular the matrix $A$ and vector $b$, can be switched instantaneously between different values.
\item The potential energy function $U(x)$ can be implemented to arbitrary accuracy.
\item Before a protocol begins, the system may be taken to be in an equilibrium distribution $\mathcal{N}[0, \beta^{-1}\|A\|^{-1}\mathbb{I}]$.
\end{itemize}

\subsection{Numerical simulations}

We outline the methods used for our numerical simulations of the thermodynamic algorithms. In general, we simulated the overdamped system dynamics because the performance is similar to the underdamped case, and the overdamped system is more numerically stable. The ODL equation, $dx = - 𝒜x + 𝒩[0, ℬ dt]$, is often written using an Itô integral \cite{gardiner1985handbook}:
\begin{equation}
\label{eq:ou-process-ito}
x(t) = e^{-\mathcal{A} t}x_0 + \int_{0}^t dt' e^{-\mathcal{A}(t-t')} Ld W_t,
\end{equation}
where $L L^⊺=ℬ$, and in this form it is apparent that the deterministic and stochastic parts of the evolution (the first and second terms above) can be evaluated separately. As this is a Gaussian process, the corresponding Fokker-Planck dynamics are fully captured by the behavior of the first and second moments, which can be evaluated directly using the well-known solution to the Ornstein-Uhlenbeck (OU) process. The JAX~\cite{frostig2018compiling} library was used to simulate the system at high dimensions, leveraging efficient implementations of matrix exponentials, diagonalization, and convolution to evaluate the various terms in the solution of the OU process. In what follows, we describe the timing model employed for benchmarking our thermodynamic algorithm, assuming an implementation using electrical circuits.

\subsubsection{Timing model}

To obtain the comparisons to other digital methods, we considered the following procedure to run the TA on electrical hardware. For more details on our model for the hardware implementation we refer the reader to the Supplemental Information. We take the $RC=1/γ= 1 μs$, which sets the characteristic timescale of the thermodynamic device. The determinant estimation procedure is excluded here for clarity, as it involves directly measuring work, which may involve a more complicated hardware proposal.
\begin{enumerate}
\item Compute the values of the resistors $\{R_{ij}, R'_i\}$ entering the matrix $\mathcal{J}$ that encodes the $A$ matrix.
\item Digital-to-analog (DAC) conversion of the $\mathcal{J}$ matrix and the $b$ vector with a given bit-precision.
\item Let the dynamics run for $t_0$ (the equilibration time). Note that for simulations this time was chosen heuristically by exploring convergence in the solutions of the problem of interest.
\item Switch on the integrators (and multipliers for the inverse estimation) and let the system evolve for time $\tau$.
\item Analog-to-digital (ADC) conversion of the solution outputted from the integrators sent back to the digital device.
\end{enumerate}
For step 1, we measured the time for the digital operation to be performed, and for the other steps we estimated the time, based on the following assumptions:
\begin{itemize}
\item 16 bit-precision
\item 5000 ADC/DAC channels with a sampling rate: 250 Msamples/s.
\item $R = 10^3 \,\Omega$, $C= 1 \,\text{nF}$, which means $RC = 1 \mu \text{s}$ is the characteristic timescale of the system.
\end{itemize}
Finally, note that in all cases that were investigated, the dominant contribution to the total runtime was the digital compilation step. This step includes $O(d^2)$ operations and involves conversion of the matrix $A$ to $𝒥$, which is detailed in the Supplemental Information. Hence some assumptions about the DAC/ADC may be relaxed and the total thermodynamic runtime would be similar. The $RC$ time constant may also be reduced to make the algorithm faster.

\bibliography{thermo_new.bib}

\pagebreak

\begin{appendix}

\begin{center}
\Large Supplemental Information for ``Thermodynamic Linear Algebra''
\end{center}

In this Supplemental Information we provide a derivation of the stationary distributions and equilibration times for systems governed by overdamped and underdamped Langevin equations. We also formulate rigorous statements of ergodicity (the time average of an observable converging to its ensemble average in the long-time limit). Ergodicity is proved for the first and second moments of the stationary distribution, and the integration times necessary for ergodicity of these observables are identified. In addition, we provide further details on our algorithms for determinant estimation and solving the Lyapunov equation.

\section{Thermodynamic Framework}

The algorithms detailed in the main text were developed using simple thermodynamic arguments, and the analysis of their performance makes use of ideas from the field of stochastic thermodynamics. We provide a brief summary of the relevant concepts here.

\subsection{Equilibrium}
Suppose that a system's state can be completely described (on a microscopic level) by a vector of generalized coordinates $x = (x_1 \dots x_d)^\intercal$ and canonically conjugate momenta $p = (p_1 \dots p_d)^\intercal$, and that the energy of the system is given by a Hamiltonian function $H(x,p)$
\begin{equation}
E = H(x,p).
\end{equation}
A statistical ensemble is an imaginary collection of copies of this system, each of which has its own coordinates $(x,p)$ and energy $E$. The relative density of copies in different parts of the coordinate space is described by a non-negative density function $f(x,p)$, satisfying a normalization constraint $\int dr f(x,p) = 1$. Ensembles describe macroscopic states (macrostates) of complex systems, where the microscopic state (microstate) is not precisely known. An observable quantity which depends on the microstate $O(r)$ has an \emph{ensemble average} $\braket{O}$, defined as
\begin{equation}
\braket{O} = \int dx \, dp \, f(x,p) O(x,p).
\end{equation}
One such ensemble, called the \emph{canonical ensemble}, has the density function \cite{callen1998thermodynamics}
\begin{equation}
f_\beta(r) = \frac{1}{Z} \exp\left(-\beta H(x,p)\right),
\end{equation}
where $\beta = 1/k_B T$ is the inverse temperature and $Z$ is partition function
\begin{equation}
Z = \int dx\, dp \, \exp\left( -\beta H(x,p)\right).
\end{equation}
The canonical ensemble is the macrostate of a system in \emph{thermal equilibrium} with an environment at temperature $T$, which is large and does not change in response to changes in the system, a heat reservoir.

In nature, we do not encounter anywhere the infinite collection of copies implied by the statistical ensemble; instead we have access to one particular instance of a system at different moments in time. The usage of the ensemble concept is justified by the property of \emph{ergodicity} observed empirically in many systems. A system is ergodic insofar as its observable quantities may be averaged over time to yield the appropriate ensemble averages, that is \cite{zwanzig2001nonequilibrium}
\begin{equation}
\overline{O} =   \frac{1}{\tau} \int_{t_0}^\tau dt \, O(t) = \braket{O}
\end{equation}
for any $t_0$ and sufficiently long $\tau$. This characterization of ergodicity should not be taken as a definition, but more of a template; in practice, we will only claim (and prove) that some particular observable $O$ is ergodic when averaged over a sufficiently long length of time, which depends on the observable.

Of interest in this work are Hamiltonians of the form
\begin{equation}
H(x,p) = U(x) + \frac{1}{2M} p^\intercal p,
\end{equation}
where $M\in \mathbb{R}^+$ and the potential is of the form
\begin{equation}
\label{eq:full-potential}
U(x) = \frac{1}{2} x^\intercal A x - b^\intercal x,
\end{equation}
for a symmetric, positive-definite $A\in \mathbb{R}^{d\times d}$ and $b\in \mathbb{R}^d$. The above can be written
\begin{equation}
\label{eq:full-hamiltonian}
H(x,p) = \frac{1}{2}(x - A^{-1}b)^\intercal A(x - A^{-1} b) - \frac{1}{2}b^\intercal A^{-1} b + \frac{1}{2M}p^\intercal p.
\end{equation}
As a rule, the Hamiltonian may be increased or decreased by a constant without any observable effect, and equivalently any term appearing in the Hamiltonian that does not depend on $x$ or $p$ may be absorbed into the partition function. In this case, we may take
\begin{equation}
\label{eq:full-hamiltonian-shifted}
H(x,p) = \frac{1}{2}(x - A^{-1} b)^\intercal A(x - A^{-1} b)+\frac{1}{2M}p^\intercal p,
\end{equation}
so the density function of the canonical ensemble is
\begin{equation}
f(x,p) = \frac{1}{Z}\exp \left( -\frac{\beta}{2}(x - A^{-1} b)^\intercal A(x - A^{-1} b) + \frac{1}{2M} p^\intercal p \right),
\end{equation}
which we often write as
\begin{equation}
\label{eq:full-hamiltonian-gibbs-dist}
x \sim \mathcal{N}[A^{-1}b, \beta^{-1}A^{-1}], \: \: p \sim  \mathcal{N}\left[0, \beta^{-1} M \mathbb{I}\right].
\end{equation}

\subsection{Stochasticity}
The property of ergodicity outlined previously is very strict, and is generally not valid on small timescales, where the system may display rich behavior which is not captured by the long time averages. This behavior is the province of stochastic thermodynamics. Due to the largeness of the heat reservoir, its interactions with the system result in effectively random displacements of the coordinate vector, which we model by Brownian motion. In the absence of any other dynamics, (anisotropic) Brownian motion is given by the equation \cite{doob1942brownian}
\begin{equation}
\label{eq:pure-brownian-sde}
dx = \mathcal{N}[0, B dt],
\end{equation}
where $B\in \mathbb{R}^{d\times d}$ is symmetric and positive definite, and we abuse notation by writing $\mathcal{N}[\mu, \Sigma]$ to denote a vector drawn from this distribution. Equation \eqref{eq:pure-brownian-sde} is self-consistent in that subdividing the time step $dt$ into two equal parts gives
\begin{equation}
x(t+dt) = x(t) + \mathcal{N}[0,B\,dt/2] +\mathcal{N}[0,
\end{equation}
and because adding two random normal vectors can be accomplished by summing their means and covariance matrices, we recover Eq. \eqref{eq:pure-brownian-sde}. If $x(0)$ is a random normal vector with mean $\mu_0$ and covariance matrix $\Sigma_0$ then, under Eq. \eqref{eq:pure-brownian-sde}, the distribution at a later time $t$ is
\begin{equation}
x(t) \sim \mathcal{N}\left[\mu_0, \Sigma_0 + t B\right].
\end{equation}
Note that this implies that the limit $\lim_{t\to \infty} \Sigma(t)$ does not exist, meaning there is no stationary distribution for the pure Brownian process. The overdamped Langevin (ODL) equation is a useful generalization of Brownian motion which includes a drift term determined by the potential \cite{lemons1997paul, hatano2001steady}:
\begin{equation}
\label{eq:general-odl}
dx = -\frac{1}{\gamma} \nabla U(x) dt +\frac{1}{\gamma}\mathcal{N}[0, B dt],
\end{equation}
where $\gamma \in \mathbb{R}^+$ is called the damping constant. When $V$ is given by Eq. \eqref{eq:full-potential}, the ODL equation can be put in a particularly simple form by making a change of variables $y = x - A^{-1} b$, in terms of which the potential is $V(y) = \frac{1}{2} y^\intercal A y$. We also define $\mathcal{A} = \gamma^{-1}A$, as well as $\mathcal{B} = \gamma^{-2}B$, and now Eq. \eqref{eq:general-odl} reads
\begin{equation}
\label{eq:odl-ou-process}
dy = - \mathcal{A}y + \mathcal{N}[0, \mathcal{B} dt],
\end{equation}
which is the defining equation of the Ornstein-Uhlenbeck (OU) process (with the requirement that $A$ is positive definite). The OU process is a Gaussian process \cite{mackay1998introduction}, meaning that if $y$ is initially a Gaussian random variable then it will be Gaussian at all times. Suppose this is the case, and $x_0 \sim \mathcal{N}[\mu_0,\Sigma_0]$. The mean $\mu(t) = \braket{y(t)}$ is given by
\begin{equation}
\label{eq:ou-solution-mu}
\mu(t) = e^{-\mathcal{A} t}\mu_0.
\end{equation}
The second moments of $y(t)$ can also be found explicitly. Because the process is indexed both by the components of $y$ and by the time, the covariances are collected in a matrix-valued function of two time parameters,
\begin{equation}
G(t,s) = \braket{y(t)y^\intercal(s)},
\end{equation}
which is called the correlation function. For the OU process,
\begin{equation}
\label{eq:ou-solution-G}
G(t,s) = e^{-\mathcal{A}t}\Sigma_0 e^{-\mathcal{A}^\intercal t} + \int_{0}^{\text{min}(t,s)} dt' \, e^{-\mathcal{A} (t-t')} \mathcal{B} e^{- \mathcal{A} (s - t')},
\end{equation}
and upon setting $s=t$, we find that the covariance matrix is given by
\begin{equation}
\Sigma(t) = e^{-\mathcal{A} t}\Sigma_0 e^{-\mathcal{A}^\intercal t} + \int_{0}^{t} dt' \, e^{-\mathcal{A} (t-t')} \mathcal{B} e^{- \mathcal{A} (t - t')}.
\end{equation}
The OU process has a unique stationary distribution, which is a Gaussian $\mathcal{N}[0,\Sigma_\text{s}]$, where $\Sigma_\text{s}$ is the unique solution to the equation
\begin{equation}
\label{eq:ou-stationary-sigma}
\mathcal{A} \Sigma_\text{s} + \Sigma_\text{s} \mathcal{A}^\intercal= \mathcal{B}.
\end{equation}
Notice that if $\mathcal{B} = \gamma^{-2}B = 2\gamma^{-1}\beta^{-1}\mathbb{I}$, then Eq. \eqref{eq:ou-stationary-sigma} implies that the stationary solution to the ODL equation is
\begin{equation}
x \sim \mathcal{N}[A^{-1} b, \beta^{-1}A^{-1}],
\end{equation}
agreeing with Eq. \eqref{eq:full-hamiltonian-gibbs-dist}. When the system is in its stationary distribution, the correlation function can be simplified
\begin{equation}
G_\text{s}(t, s) = e^{-\mathcal{A} (s-t)}\Sigma_\text{s},
\end{equation}
when $s \geq t$, and $G_\text{s}(t,s) = G_\text{s}(s,t)^\intercal$. These results are quite useful due to the wide variety of cases where the OU process is applicable. For example, Eq. \eqref{eq:general-odl} may be modified to account for inertial effects, resulting in the underdamped Langevin (UDL) equations, which are
\begin{equation}
\label{eq:udl-x-general}
dx  = \frac{1}{M} p \, dt,
\end{equation}
\begin{equation}
\label{eq:udl-p-general}
dp = -\nabla V \, dt- \frac{\gamma}{M} p \, dt + \mathcal{N}[0,B dt],
\end{equation}
where $M, \gamma \in \mathbb{R}^+$. Taking the potential in \eqref{eq:full-potential}, Eq. \eqref{eq:udl-p-general} reads
\begin{equation}
\label{eq:udl-p-device}
dp = -(A x - b)\, dt- \frac{\gamma}{M} p \, dt + \mathcal{N}[0,B dt].
\end{equation}
In order to put the UDL equations into the form of an OU process, we first define dimensionless coordinates
\begin{equation}
\tilde{y} =\sqrt{\frac{\gamma^2 \beta}{M}}y, \: \: \: \tilde{p} = \sqrt{\frac{\beta}{M}}p,
\end{equation}
and define the vector $r$ as the concatenation of the dimensionless position and momentum vectors
\begin{equation}
\label{r-def-app}
r = \left(\tilde{y}_1, \dots \tilde{y}_d, \dots \tilde{p}_1, \dots \tilde{p}_d\right)^\intercal.
\end{equation}
Then we have
\begin{equation}
\label{UDL-vec}
dr = -\mathcal{A} r\, dt + \mathcal{N}[0, \mathcal{B} dt],
\end{equation}
where the matrices $\mathcal{A}$ and $\mathcal{B}$ may be written as block matrices, each with four $d \times d$ blocks
\begin{equation}
\mathcal{A} = \begin{pmatrix}
0 & -\gamma M^{-1}\mathbb{I} \\
\gamma^{-1}A & \gamma M^{-1} \mathbb{I}
\end{pmatrix}, \: \: \: \:
\mathcal{B} = \begin{pmatrix}
0 & 0  \\
0 & \beta M^{-1} B
\end{pmatrix}.
\end{equation}
If we take $B = 2 \gamma \beta^{-1} \mathbb{I}$, then Eq. \eqref{eq:ou-stationary-sigma} gives
\begin{equation}
x \sim \mathcal{N}[A^{-1}b, \beta^{-1}A^{-1}], \: \: p \sim  \mathcal{N}\left[0, \beta^{-1}M \mathbb{I}\right],
\end{equation}
which again agrees with Eq. \eqref{eq:full-hamiltonian-gibbs-dist}.

\subsection{Thermodynamics}
Thermodynamics is concerned with the relationships between macroscopically observable quantities which are well-characterized by their equilibrium ensemble averages. For example, the average energy is defined as
\begin{equation}
\braket{E} = \int dx H(x,p) f_\beta(x,p).\\
\end{equation}
For the Hamiltonian given in Eq. \eqref{eq:full-hamiltonian}, a Gaussian integral yields
\begin{equation}
\label{eq:average-energy}
\braket{E} =  d \, \beta^{-1}-\frac{1}{2}b^\intercal A^{-1} b.
\end{equation}
According to the results of the last section, the system asymptotically approaches the canonical equilibrium distribution (in both the overdamped and underdamped models) if $B = 2\gamma \beta^{-1}\mathbb{I}$, which can be used to eliminate temperature from Eq. \eqref{eq:average-energy}, giving
\begin{equation}
\label{eq:average-energy-2}
\braket{E} = \frac{d\left|B\right|}{2 \gamma} -\frac{1}{2}b^\intercal A^{-1} b.
\end{equation}
This relationship between temperature, damping, and noise variance is a statement of the fluctuation-dissipation relation which occurs frequently in stochastic thermodynamics. It is often explained with reference to the fact that both damping and noise arise from interactions with the bath, and stronger coupling between the system and bath tends to increase both noise and damping. Equation \eqref{eq:average-energy-2} provides an intuitive view of the fluctuation-dissipation relation, namely that energy increases with more noise and decreases with more damping. Having evaluated the average energy, we next consider changes in energy, which are separated into work and heat. The work is defined as the part of the change in energy which is due to the variation in time of one or more control parameters. If $\lambda_1 \dots \lambda_k$ are the complete set of control parameters, then the work done on the system in a time interval $dt$ is \cite{seifert2008stochastic} given by:
\begin{equation}
dW = \sum_{i=1}^k \frac{\partial H}{\partial \lambda_i} \frac{d \lambda_i}{d t} dt = -\sum_{i=1}^k X_i \dot{\lambda}_i,
\end{equation}
where we have defined the forces conjugate to the control parameters as
\begin{equation}
X_i = - \frac{\partial H}{\partial \lambda_i}.
\end{equation}
The heat is defined as the remaining part of the energy change, associated with the change in the system's coordinates
\begin{equation}
dQ = \sum_{i=1}^d \frac{\partial H}{\partial x_i} \dot{x}_i dt + \frac{\partial H}{\partial p_i} \dot{p}_i dt ,
\end{equation}
and in terms of these we can state the first law of thermodynamics,
\begin{equation}
\frac{d \braket{E}}{dt} = \frac{\braket{dW}}{dt} + \frac{\braket{dQ}}{dt}.
\end{equation}
The second law can be expressed in terms of the Helmholtz free energy, which is defined as \cite{callen1998thermodynamics}
\begin{equation}
F = -\beta^{-1} \ln(Z),
\end{equation}
and can be interpreted as the amount of energy in the system that can be converted to work at constant temperature $\beta^{-1}$. For the Hamiltonian given in Eq. \eqref{eq:full-hamiltonian}, we get
\begin{equation}
F = -\frac{1}{2} b^\intercal A^{-1} b+ \frac{1}{2\beta} \ln \left(\left| A\right| \right) -\frac{d}{\beta} \ln\left(\frac{\beta}{2\pi \sqrt{M}}\right),
\end{equation}
and so two states with the same temperature have free energy difference
\begin{equation}
\Delta F = F_2 - F_1 = \frac{1}{2}(b_1^\intercal A_1^{-1}b_1-b_2^\intercal A_2^{-1}b_2) + \frac{1}{2\beta} \ln \left( \frac{\left|A_2\right|}{\left|A_1\right|}\right).
\end{equation}
The second law reads
\begin{equation}
\label{eq:second-law-helmholtz}
\braket{W} \geq \Delta F.
\end{equation}
Whereas Eq. \eqref{eq:second-law-helmholtz} is an inequality, Jarzynski identifies the following equality \cite{jarzynski1997nonequilibrium}
\begin{equation}
\label{eq:jarzynski-relation}
e^{-\beta \Delta F} = \Braket{e^{-\beta W}},
\end{equation}
and it follows from Jensen's inequality \cite{jarzynski2011equalities} that the latter implies Eq. \eqref{eq:second-law-helmholtz}. In the limit of an infinitely slow (quasistatic) process we would have $\braket{W} = W =  \Delta F$, but in general there will be some excess (dissipated) work,
\begin{equation}
\braket{W_\text{ex}} = \braket{W} - \Delta F.
\end{equation}
Intuitively, there is more dissipation when the system is driven between very different states in a short period of time, as friction-like effects are more significant when the system changes quickly. This notion is made rigorous by introducing the thermodynamic metric tensor $g$ \cite{crooks2007measuring, shenfeld2009minimizing}
\begin{equation}
\label{eq:friction-tensor-def}
g_{ij}(\lambda) = \Braket{\left(X_i - \Braket{X_i}\right)\left(X_j - \Braket{X_j}\right)},
\end{equation}
where the average is over the canonical ensemble with inverse temperature $\beta$ and control parameters $\lambda$. The thermodynamic length $\mathcal{L}$ of a trajectory in the control parameter space $\lambda(t)$ is then given by
\begin{equation}
\mathcal{L} = \int_0^\tau dt \sqrt{\dot{\lambda}^\intercal g \dot{\lambda}}.
\end{equation}
Suppose that a system is driven from one state to another, and then back again by the reverse process, and that it is allowed to come to equilibrium $N$ times during both the forward and reverse process. It has been shown \cite{crooks2007measuring} that the excess work during the combined forward and reverse process (called the hysteresis) is lower bounded as
\begin{equation}
W_\text{ex} \geq \frac{\mathcal{L}^2}{N}.
\end{equation}
Interestingly, the thermodynamic length is also equivalent to the Fisher-Rao distance, which is \cite{amari2000methods}
\begin{equation}
\mathcal{L} = \beta \int_{0}^\tau dt \sqrt{ \int dx \frac{\dot{f}(x)^2}{f(x)}}.
\end{equation}
The above relations imply that the more distinguishable the initial and final distributions, and the less time allowed to evolve between them, the more work will be dissipated.

Equation \eqref{eq:jarzynski-relation} suggests that the free energy difference can be estimated by repeatedly transforming one potential $U_1$ into another $U_2$ and measuring the work done in the process ($n_J$ times, say), then  taking an average. The estimate of the free energy obtained in this way is called the Jarzynski estimator $\Delta \hat{F}_J$
\begin{equation}
\Delta F \approx \Delta \hat{F}_J = \frac{1}{n_J}\sum_{j=1}^{n_J} e^{-\beta W_j}.
\end{equation}
The Jarzynski estimator is theoretically sound, but in practice it can be very slow to converge as it is dominated by terms with large negative work, and there are very few such terms in the average because the occurrence of large negative work is unlikely \cite{jarzynski2006rare}. This problem is mitigated somewhat by using an estimator which makes use of both the forward and reverse process. The Bennett Acceptance Ratio (BAR) estimator $\Delta \hat{F}_\text{BAR}$ is defined implicitly as the solution to \cite{shirts2003equilibrium, bennett1976efficient}
\begin{equation}
\sum_{i=1}^{N} \frac{1}{1 +  e^{ \beta W_{\text{(f)}i} - \beta \Delta \hat{F}_\text{BAR}}}  =   \sum_{i=1}^{N} \frac{1}{1 +  e^{ -\beta W_{\text{(r)}i} +\beta \Delta \hat{F}_\text{BAR}}},
\end{equation}
where it is assumed that there are $N$ realizations of both the forward and reverse process, resulting in work measurements of $W_{\text{(f)}i}$ and $W_{\text{(r)}i}$. This estimator can be generalized to the Multistate Bennett Acceptance Ratio \cite{shirts2008statistically} (MBAR), where intermediate equilibrium states are used between the initial and final states \cite{shirts2008statistically}. Interestingly, the variance of any unbiased estimator of the free energy difference is lower bounded in terms of the thermodynamic length between the terminal states and the total number of observations $n$, which corresponds to the number of times that equilibrium must be reached during the protocol \cite{shenfeld2009minimizing}
\begin{equation}
\text{Var}(\Delta \hat{F}) \geq \frac{\beta^{-2} \mathcal{L}^2}{n}.
\end{equation}

\section{Analysis of the Overdamped Regime}

\subsection{Stationary Distribution}
The overdamped Langevin equation reads
\begin{equation}
\label{eq:ODL-app}
dx = - \frac{1}{\gamma}(A x - b) dt + \frac{1}{\gamma}\mathcal{N}\left[0,B\, dt\right].
\end{equation}
We set $B = 2 \gamma \beta^{-1} \mathbb{I}$ (see Methods), where $\beta =(k_B T)^{-1}$ is the inverse temperature of the environment. We then change variables to $y = x - A^{-1}b$, and have the transformed equation
\begin{equation}
\label{ODL-transformed-app}
dy = - \frac{1}{\gamma}A y dt + \mathcal{N}\left[0, \frac{2}{\beta \gamma} \mathbb{I}\, dt\right],
\end{equation}
which can be written as a vector Ornstein-Uhlenbeck (OU) process \cite{gardiner1985handbook}
\begin{equation}
\label{ODL-transformed-app}
dy = - \mathcal{A} y dt + \mathcal{N}[0, \mathcal{B}\, dt],
\end{equation}
where $\mathcal{A} = \gamma^{-1}A$ and $\mathcal{B} = 2 \gamma^{-1}\beta^{-1}\mathbb{I}$. The stationary distribution for $y$ then has mean zero (which follows from Eq. \eqref{eq:ou-solution-mu}) and its variance $\Sigma_\text{s}$ satisfies Eq. \eqref{eq:ou-stationary-sigma},
\begin{equation}
\mathcal{A} \Sigma_\text{s} + \Sigma_\text{s} \mathcal{A} = \mathcal{B}.
\end{equation}
In this case, we have
\begin{equation}
\frac{1}{\gamma}A \Sigma_\text{s} + \frac{1}{\gamma}\Sigma_\text{s} A = \frac{2}{\beta \gamma}\mathbb{I},
\end{equation}
which is satisfied for the choice
\begin{equation}
\label{ODL-stationary-sigma}
\Sigma_\text{s} = \beta^{-1} A^{-1}.
\end{equation}
The uniqueness of this solution is guaranteed, because Eq. \eqref{ODL-stationary-sigma} always has a unique solution when $A$ is positive definite. Transforming back to the original coordinates, we see that at equilibrium $x$ is distributed as
\begin{equation}
\label{eq:stationary-x-linear-system}
x \sim \mathcal{N}[A^{-1}b, \beta^{-1} A^{-1}].
\end{equation}
Alternatively we may take $B = 2\gamma\beta^{-1}R$ (or $\mathcal{B} = 2\gamma^{-1} \beta^{-1}R$) for some symmetric positive-definite matrix $R$. Then we would have stationary covariance matrix satisfying
\begin{equation}
A \Sigma_\text{s} +\Sigma_\text{s} A =2\beta^{-1}R,
\end{equation}
which is a Lyapunov equation. In the latter case, the stationary mean of $x$ is still $A^{-1} b$.
\subsection{Equilibration and Correlation Time}
We assume that the initial distribution is Gaussian. After some time has passed, the system will be in the equilibrium distribution described by Eq. \eqref{eq:stationary-x-linear-system}. Here we quantify the amount of time one should wait (called the equilibration time) to ensure that the distribution is arbitrarily close to the equilibrium distribution. Note that because $A$ is symmetric positive definite, the spectral norm of $e^{-A t/\gamma}$ is
\begin{equation}
\label{exp-A-norm-bound-app}
\left\|e^{-A t/\gamma}\right\| = e^{-\alpha_\text{min} t/\gamma}.
\end{equation}
As Eq. \eqref{ODL-transformed-app} describes an OU process, the mean evolves in time according to
\begin{equation}
\braket{y(t)} = e^{-A t/\gamma}\braket{y(0)}.
\end{equation}
Let $\alpha_1 \dots \alpha_d$ be the eigenvalues of $A$. The norm of $\braket{y}$ is exponentially bounded,
\begin{equation}
\label{odl-mean-exp-bound-app}
\left\|\braket{y(t)}\right\| =\left \|e^{- A t/\gamma} \braket{ y(0)}\right\| \leq e^{-\alpha_\text{min} t/\gamma} \left\| \braket{y(0)}\right\|.
\end{equation}
Given $\delta_{\mu 0}>0$, in order to have $\left\|\braket{x(t_0)} -A^{-1} b\right \| \leq \delta_{\mu 0}$, it is sufficient to require
\begin{equation}
t_0 \geq \frac{\gamma}{\alpha_\text{min}} \ln \left( \frac{\left\|\braket{x(0)} - A^{-1} b\right\|}{\delta_{\mu0}}\right).
\end{equation}
We assume the initial distribution has mean zero, $\braket{x(0)}=0$. Defining the condition number $\kappa = \alpha_\text{max}/\alpha_\text{min}$, and the relative error tolerance $\varepsilon_{\mu 0} = \delta_{\mu 0 }/\|A^{-1} b\|$ it is sufficient to require
\begin{equation}
\label{eq:odl-mean-equilibration-time-app}
t_0 \geq \kappa \gamma \|A\|^{-1} \ln \left(\kappa \varepsilon^{-1}_{\mu 0}\right).
\end{equation}
The correlation matrix is defined by
\begin{equation}
G(t,s) = \braket{[x(t) - \braket{x(t)}][x(s) - \braket{x(s)}]^\intercal}
\end{equation}
For $t=s$, the correlation matrix reduces to the covariance matrix
\begin{equation}
G(t,t) = \braket{[x(t) - \braket{x(t)}][x(t) - \braket{x(t)}]^\intercal} = \Sigma(t)
\end{equation}
For an OU process, the dynamics of the correlation matrix (and the covariance matrix in particular) can be expressed in terms of the following propagator function
\begin{equation}
\mathcal{P}_t(X) = e^{-A t/\gamma} X e^{-A^\intercal t/\gamma}.
\end{equation}
The covariance matrix for $y$ is then given by
\begin{equation}
\label{sigma-time-dep-app}
\Sigma(t) = \mathcal{P}_t(\Sigma_0) + \int_0^t dt' \mathcal{P}_{t'} (\mathcal{B}),
\end{equation}
and for any $t > t_0$ we have
\begin{equation}
\label{sigma-error-time-dep-app}
\Sigma(t) - \Sigma(t_0)=\mathcal{P}_{t}(\Sigma_0) -\mathcal{P}_{t_0}(\Sigma_0) + \int_{t_0}^{t} dt' \mathcal{P}_{t'} (\mathcal{B}),
\end{equation}
Taking the limit as $t \to \infty$, we evaluate
\begin{align}
\left\|\Sigma_\text{s} - \Sigma(t_0) \right\| &= \left\| \Sigma_\text{s} -\mathcal{P}_{t_0}(\Sigma_0) + \int_{t_0}^{\infty} dt' \mathcal{P}_{t'} (\mathcal{B})\right\| \\
& =  \left\|\mathcal{P}_{t_0}( \Sigma_\text{s} - \Sigma_0) + \int_{t_0}^{\infty} dt' \mathcal{P}_{t'} (\mathcal{B})\right\| \\
& \leq
\left\| \mathcal{P}_{t_0} (\Sigma_\text{s} - \Sigma_0) \right\|+ \int_{t_0}^\infty dt' \left\| \mathcal{P}_{t'}(\mathcal{B})\right\|\\
& \leq
\left\| e^{-A t/\gamma}\right\|^2 \cdot \left\| \Sigma_0 - \Sigma_\text{s}\right\|+ \left\|\mathcal{B}\right\|\int_{t_0}^\infty dt' \left\| e^{-A t'/\gamma}\right\|^2 \\
& \leq     e^{- 2\alpha_\text{min} t_0/\gamma }\left\| \Sigma_0 - \Sigma_\text{s}\right\|+ \left\|\mathcal{B}\right\|\int_{t_0}^\infty dt' e^{-2 \alpha_\text{min} t'/\gamma} \\
& =  e^{- 2\alpha_\text{min} t_0 }\left\| \Sigma_0 - \Sigma_\text{s}\right\|+\frac{\gamma}{2\alpha_\text{min}}e^{-2 \alpha_\text{min}t_0/\gamma} \left\|\mathcal{B}\right\| \\
& =  e^{- 2\alpha_\text{min} t_0/\gamma }\left(\left\| \Sigma_0 - \Sigma_\text{s}\right\|+\frac{\gamma}{2  \alpha_\text{min}}\|\mathcal{B}\|\right).
\end{align}
So in order to have $\|\Sigma(t_0) - \Sigma_\text{s}\| \leq \delta_
\Sigma$, we take
\begin{equation}
t_0   \geq  \frac{\gamma}{2 \alpha_\text{min}}\ln\left(\frac{\|\Sigma_0 - \Sigma_\text{s}\| + (\gamma / 2\alpha_\text{min})\|\mathcal{B}\|}{\delta_{\Sigma 0}}\right).
\end{equation}
We use the fact that $\Sigma_\text{s} = \beta^{-1} A^{-1}$, and $\mathcal{B} = 2 \gamma^{-1} \beta^{-1} \mathbb{I}$. Also, we assume the system is initially at equilibrium in a potential proportional to identity $A_0 = \|A\|\mathbb{I}$, so $\Sigma_0 = \beta^{-1} \|A\|^{-1}\mathbb{I}$. We have
\begin{align}
t_0   &\geq  \frac{\kappa\gamma}{2 \|A\|}\ln\left(\frac{\beta^{-1}\|(\|A^{-1}\|\mathbb{I}) - A^{-1}\| + \kappa \beta^{-1}/\|A\|}{\delta_{\Sigma 0}}\right)
\end{align}
We can give a looser requirement
\begin{align}
t_0   &\geq  \frac{\kappa\gamma}{2\|A\| }\ln\left(\frac{2\kappa\beta^{-1}\|A^{-1}\|}{\delta_\Sigma }\right).
\end{align}
Defining the variance relative error tolerance $\varepsilon_{\Sigma 0} = \delta_{\Sigma 0}/\beta^{-1} \|A^{-1}\|$, we have
\begin{align}
\label{eq:odl-variance-equilibration-time-app}
t_0   &\geq  \frac{\kappa\gamma}{2 \|A\|}\ln\left(\frac{2\kappa}{\varepsilon_{\Sigma 0}}\right).
\end{align}
Combining Eqs.~\eqref{eq:odl-mean-equilibration-time-app} and \eqref{eq:odl-variance-equilibration-time-app} gives
\begin{equation}
t_0 \geq \max \left\{\frac{\kappa \gamma}{\|A\|} \ln \left(\kappa \varepsilon^{-1}_{\mu 0}\right),\frac{\kappa\gamma}{2\|A\| }\ln\left(2\kappa\varepsilon_{\Sigma 0}^{-1}\right)\right\},
\end{equation}
Finally, we define the relaxation time $\tau_\text{r} = \gamma/\|A\|$, and we have
\begin{equation}
t_0 \geq \max \left\{\kappa \tau_\text{r}\ln \left(\kappa \varepsilon^{-1}_{\mu 0}\right),\frac{1}{2} \kappa \tau_\text{r}\ln\left(2\kappa\varepsilon_{\Sigma 0}^{-1}\right)\right\}.
\end{equation}

More generally we consider the case where $\mathcal{B} = 2 \gamma^{-1} \beta^{-1} R$, and $\|R\|=1$. Now the variance equilibration time is
\begin{align}
t_0  & \geq  \frac{\gamma}{2 \alpha_\text{min}}\ln\left(\frac{\|\Sigma_0 - \Sigma_\text{s}\| + (\gamma \kappa / 2\|A\|)\|\mathcal{B}\|}{\delta_{\Sigma 0}}\right)\\
&= \frac{1}{2}\kappa\tau_\text{r}\ln\left(\frac{\|\Sigma_0 - \Sigma_\text{s}\| + \kappa \beta^{-1} \|A\|^{-1}
}{\delta_{\Sigma 0}}\right)
\end{align}
However, without knowing the solution to the Lyapunov equation the above cannot be sharply bounded, so in this case we are left with the form
\begin{equation}
t_0 \geq \max \left\{\kappa \tau_\text{r}\ln \left(\kappa \varepsilon^{-1}_{\mu 0}\right),\frac{1}{2}\kappa\tau_\text{r}\ln\left(\frac{\|\Sigma_0 - \Sigma_\text{s}\| + \kappa \beta^{-1} \|A\|^{-1}
}{\delta_{\Sigma 0}}\right)
\right\}.
\end{equation}
Once the equilibration time has passed, we may consider the correlation matrix of the stationary system
\begin{equation}
G_\text{s}(\tau) = \lim_{t \to \infty} G(t, t+\tau),
\end{equation}
which is given by
\begin{equation}
G_\text{s}(\tau) = e^{-A \tau/\gamma} \Sigma_\text{s}.
\end{equation}
\begin{equation}
\|G_\text{s}(\tau)\| \leq \|e^{-A t/\gamma}\| \cdot \|\Sigma_\text{s}\| \leq e^{- \alpha_\text{min} \tau/\gamma}  \left\|\Sigma_\text{s}\right\|.
\end{equation}
Suppose we sample once per interval $\tau_c$. To obtain samples with bounded correlation $\left\| G_\text{s} \right \| \leq \delta_G $, we need
\begin{equation}
\tau_c \geq \frac{\gamma}{\alpha_\text{min}} \ln \left(\frac{\|\Sigma_\text{s}\|}{\delta_G}\right).
\end{equation}
Defining $\varepsilon_G = \delta_G/\|\Sigma_\text{s}\|$, we get
\begin{equation}
\tau_c \geq\kappa \tau_\text{r} \ln \left(\kappa \varepsilon_G^{-1}\right).
\end{equation}
Interestingly, the correlation time is dimension-independent, ignoring any incidental dependence of condition number on dimension.

\subsection{Ergodicity of Mean}

Define the time-average of $y$
\begin{equation}
\bar{y} = \frac{1}{\tau} \int_{t_0}^{t_0 + \tau} dt y(t).
\end{equation}
We assume that $t_0$ is large enough that, to a good approximation, the system has reached equilibrium by time $t_0$. Therefore
\begin{equation}
\braket{\bar{y}} = \frac{1}{\tau} \int_{t_0}^{t_0 + \tau} dt \braket{y(t)} = 0.
\end{equation}
We now bound the covariance matrix of $\bar{y}$
\begin{align}
\left\|\braket{\bar{y}\bar{y}^\intercal}\right\| &=\left \| \frac{1}{\tau^2}\int_{t_0}^{t_0+\tau}\int_{t_0}^{t_0+\tau} dt''dt' \braket{y(t')y(t'')^\intercal}\right \| \\
&=  \left\|\frac{1}{\tau^2}\int_{t_0}^{t_0+\tau}\int_{t_0}^{t_0+\tau}  dt''dt' G_\text{s}(t''-t') \right\|\\
&=  \left\|\frac{2}{\tau^2}\int_{t_0}^{t_0 + \tau}dt''\int_{t_0}^{t''} dt' G_\text{s}(t''-t') \right\|\\
\frac{2}{\tau^2}\int_{t_0}^{t_0 + \tau}dt''\int_{t_0}^{t''} dt' \left\|G_\text{s}(t''-t') \right\|\\
\frac{2}{\tau^2}\int_{t_0}^{t_0 + \tau}dt''\int_{t_0}^{t''} dt' e^{-\alpha_\text{min}(t''-t')/\gamma} \left\| \Sigma_\text{s}\right\|\\
\frac{2\gamma}{\alpha_\text{min}\tau^2} \left\| \Sigma_\text{s}\right\|\int_{t_0}^{t_0 + \tau}dt'' 1-e^{-\alpha_\text{min}(t'' - t_0)/\gamma} \\
& \leq
\frac{2\gamma}{\alpha_\text{min}\tau^2} \left\| \Sigma_\text{s}\right\|\int_{t_0}^{t_0 + \tau}dt''1 \\
\frac{2\gamma}{\alpha_\text{min}\tau} \left\| \Sigma_\text{s}\right\|
\\ & = \frac{2 \gamma \kappa^2}{\beta\|A\|^2\tau}\\
& =  \frac{2  \kappa^2 \tau_\text{r}}{\beta\|A\|\tau}.
\end{align}
According to Chebyshev's inequality
\begin{equation}
P(y^\intercal \braket{\bar{y} \bar{y}^\intercal}^{-1} y > k^2 ) \leq \frac{d}{k^2}.
\end{equation}
Now note that the eigenvalues of $\braket{\bar{y} \bar{y}^\intercal}^{-1}$ are the inverses of the eigenvalues of $\braket{\bar{y} \bar{y}^\intercal}$, which are the same as the singular values as the covariance matrix is positive definite. Therefore because the eigenvalues of $\braket{\bar{y} \bar{y}^\intercal}$ are at most $\left\|\braket{\bar{y}\bar{y}^\intercal}\right\|$, the eigenvalues of $\braket{\bar{y} \bar{y}^\intercal}^{-1}$ are at least $ \left\|\braket{\bar{y}\bar{y}^\intercal}\right\| ^{-1}$. This implies
\begin{equation}
y^\intercal \braket{\bar{y} \bar{y}^\intercal}^{-1} y \geq \|y\|^2\frac{\beta\|A\|^2\tau}{2 \gamma \kappa^2 },
\end{equation}
and so we have the following proposition
\begin{equation}
\|y\|^2 > \frac{2 \gamma\kappa^2}{\beta\|A\|^2\tau}k^2 \Rightarrow
y^\intercal \braket{\bar{y} \bar{y}^\intercal}^{-1}y > k^2.
\end{equation}
Moreover, $\|Ay\|^2 \leq \|A\|^2 \|y\|^2$ so the following also holds
\begin{equation}
\|A y\|^2 > \frac{2 \gamma\kappa^2}{\beta\tau}k^2 \Rightarrow
y^\intercal \braket{\bar{y} \bar{y}^\intercal}^{-1}y > k^2.
\end{equation}
This means that
\begin{equation}
P\left(  \frac{\|A y\|^2}{\|b\|^2} > \frac{2 \gamma\kappa^2}{\|b\|^2\beta\tau}k^2\right) \leq P( y^\intercal \braket{\bar{y} \bar{y}^\intercal}^{-1}y \geq k^2) \leq \frac{d}{k^2}.
\end{equation}
Now let
\begin{equation}
\varepsilon_y = \sqrt{\frac{2 \gamma \kappa^2}{\beta\|b\|^2\tau}} k,
\end{equation}
so $k = \sqrt{\beta \|b\|^2 \tau/2\gamma \kappa^2}\delta_y$. Then we have
\begin{equation}
P\left(  \frac{\|A y\|}{\|b\|} \geq \varepsilon_y\right) \leq \frac{d}{k^2} = \frac{2 \gamma \kappa^2d}{\beta \|b\|^2 \tau \varepsilon_y^2}.
\end{equation}
We change back to the original coordinates. If we would like to have
\begin{equation}
\left \|A \bar{x} - b \right\| \leq \varepsilon_x \|b\|
\end{equation}
with probability at least $P_\varepsilon$, then we can require the integration time is at least
\begin{align}
\tau &\geq \frac{2 \gamma\kappa^2 d  }{\beta \|b\|^2 \varepsilon_x^2 (1-P_\varepsilon)} \\
&=\frac{2 \kappa^2 d\, \|A\| \tau_\text{r} }{\beta \|b\|^2 \varepsilon_x^2 (1-P_\varepsilon)} .
\end{align}
Here, we note that the energy cost $\mathcal{E}$ can be estimated as $\mathcal{E} \geq b^\intercal A^{-1} b \geq \|A\|^{-1}\|b\|^2$, which is the depth of the bottom of the potential well relative to the point $x=0$. Therefore $\mathcal{E}^{-1} \leq \|A\|\|b\|^{-2}$, allowing us to write
\begin{align}
\mathcal{E}\tau &\geq \frac{2 \kappa^2 d  }{ \varepsilon_x^2 (1-P_\varepsilon)} \beta^{-1}\tau_\text{r}.
\end{align}

\subsection{Ergodicity of Covariance Matrix}
The covariance matrix can also be estimated by a finite time integral
\begin{equation}
\overline{yy^\intercal} =\frac{1}{\tau} \int_{t_0}^{t_0 + \tau} dt' y(t') y^\intercal(t').
\end{equation}
\begin{equation}
\text{var}(\overline{y_iy_j}) = \braket{\overline{y_i y_j}^2} - \Sigma^2_{\text{s}, ij}
\end{equation}

\begin{equation}
\braket{\overline{y_i y_j}^2}  = \frac{1}{\tau^2}\int_{t_0}^{t_0 + \tau} dt'' \int_{t_0}^{t_0 + \tau} dt' \braket{y_i(t') y_j(t') y_i(t'') y_j(t'')}.
\end{equation}
Using Isserlis's theorem \cite{koopmans1995spectral},
\begin{align}
\braket{y_i(t') y_j(t') y_i(t'') y_j(t'')} &= \braket{y_i(t')y_j(t')}\braket{y_i(t'')y_j(t'')} +\braket{y_i(t')y_i(t'')}\braket{y_j(t')y_j(t'')}\\ & +\braket{y_i(t')y_j(t'')}\braket{y_j(t')y_i(t'')}\\
& = (\Sigma_{\text{s},ij})^2 + G_{\text{s},ii}(t''-t')G_{\text{s},jj}(t''-t') + G_{\text{s},ij}(t''-t')^2,
\end{align}
\begin{align}
\text{var}(\overline{y_i y_j})
&=\braket{\overline{y_i y_j}^2} -\Sigma_{\text{s},ij}^2  \\
&=  \frac{1}{\tau^2}\int_{t_0}^{t_0 + \tau} dt'' \int_{t_0}^{t_0 + \tau} dt' G_{\text{s},ii}(t''-t')G_{\text{s}.jj}(t''-t') + G_{\text{s},ij}(t''-t')^2.
\end{align}
Summing over $i$ and $j$,
\begin{align}
\sum_{ij}\text{var}(\overline{y_i y_j})
&=  \frac{1}{\tau^2}\int_{t_0}^{t_0 + \tau} dt'' \int_{t_0}^{t_0 + \tau} dt'\sum_{ij} G_{\text{s},ii}(t''-t')G_{\text{s}.jj}(t''-t') + G_{\text{s},ij}(t''-t')^2 \\
&\leq  \frac{2}{\tau^2}\int_{t_0}^{t_0 + \tau} dt'' \int_{t_0}^{t_0 + \tau} dt'\,(d+1)\|G_\text{s}(t'' - t')\|_\text{F}^2
\end{align}
where $\|\cdot\|_\text{F}$ is the Frobenius norm. As the Frobenius norm is bounded by $\sqrt{d}$ times the spectral norm, we have
\begin{align}
\sum_{ij}\text{var}(\overline{y_i y_j})
&\leq \frac{2d(d+1)}{\tau^2}\int_{t_0}^{t_0 + \tau} dt'' \int_{t_0}^{t_0 + \tau} dt'\|G_\text{s}(t'' - t')\|^2\\
& =  \frac{4d(d+1)}{\tau^2}\int_{t_0}^{t_0 + \tau} dt'' \int_{t_0}^{t''} dt'\|G_\text{s}(t'' - t')\|^2 \\
& \leq   \frac{4d(d+1)}{\tau^2}\int_{t_0}^{t_0 + \tau} dt'' \int_{t_0}^{t''} dt'e^{-2\alpha_\text{min} (t'' - t')/\gamma} \left\| \Sigma_\text{s}\right\|^2 \\
\frac{4d(d+1)}{\tau^2} \left\| \Sigma_\text{s}\right\|^2\int_{t_0}^{t_0 + \tau}dt'' \frac{\gamma-\gamma e^{-2\alpha_\text{min}(t'' - t_0)/\gamma}}{2\alpha_\text{min}} \\
& \leq
\frac{4d(d+1)}{\tau^2} \left\| \Sigma_\text{s}\right\|^2\int_{t_0}^{t_0 + \tau}dt'' \frac{\gamma}{2\alpha_\text{min}} \\
& = \frac{4 d (d+1) \,\gamma\left\| \Sigma_\text{s}\right\|^2}{\tau \alpha_\text{min}}. \\
\end{align}
Recall the generalization of Chebyshev's inequality, for an arbitrary norm $\|\cdot\|_\nu$,
\begin{equation}
P\left( \|X-\mu\|_\nu \geq k \sigma_\nu\right) \leq \frac{1}{k^2},
\end{equation}
\begin{equation}
\sigma_\nu^2 = \text{E}[\|X - \mu\|_\nu^2].
\end{equation}
We use the Frobenius norm for $X = \overline{y y^\intercal}$, resulting in
\begin{align}
\sigma_\text{F}^2 = \text{E}\left[ \left\|
\overline{yy^\intercal} - \braket{\overline{y y^\intercal}}
\right\|^2_\text{F}\right] &=
\text{E}\left[ \left\|
\overline{yy^\intercal} - \Sigma_\text{s}
\right\|_\text{F}^2\right]\\
\text{E}\left[ \sum_{ij} (\overline{y_i y_j} - \braket{y_i y_j})^2\right]\\
& = \sum_{ij} \text{var}(\overline{y_i y_j})\\
&\leq \frac{4 d (d+1) \,\gamma\left\| \Sigma_\text{s}\right\|^2}{\tau \alpha_\text{min}}
\end{align}
If we set
\begin{equation}
\delta_\Sigma^2 = k^2\frac{4 d (d+1) \,\gamma\left\| \Sigma_\text{s}\right\|^2}{\tau \alpha_\text{min}},
\end{equation}
\begin{align}
P\left( \|\overline{y y^\intercal} - \Sigma_\text{s}\|_\text{F} \geq \delta_\Sigma \right)
&\leq \frac{1}{k^2}\\
& = \frac{4 d (d+1) \gamma \|\Sigma_\text{s}\|^2}{\tau \alpha_\text{min} \delta_\Sigma^2}.
\end{align}
We define the relative error $\varepsilon_\Sigma = \delta_\Sigma/\|\Sigma_\text{s}\|$, and have
\begin{align}
P\left( \frac{\|\overline{y y^\intercal} - \Sigma_\text{s}\|_\text{F}}{\|\Sigma_\text{s}\|} \geq \varepsilon_\Sigma \right)
& \leq \frac{4 d (d+1) \gamma \|\Sigma_\text{s}\|^2}{\tau \alpha_\text{min} \varepsilon_\Sigma^2 \|\Sigma_\text{s}\|^2}\\
& =  \frac{4 d (d+1) \kappa \tau_\text{r}}{  \varepsilon_\Sigma^2 \tau }\\
\end{align}
So if we would like to have  (Frobenius) relative error of at most $\varepsilon_\Sigma$ with probability $P_\varepsilon$ then it is sufficient to allow integration time of at least
\begin{equation}
\tau \geq \frac{4 \kappa d(d+1)}{(1-P_\varepsilon)\varepsilon_\Sigma^2}\tau_\text{r}.
\end{equation}

\subsection{Hardware implementation}\label{app:hw}
Here we describe an electronic device comprised of $d$ coupled $RC$ cells that maps to the overdamped Langevin process described previously. The equation of motion for voltages $\mathbf{v} = (v_1, v_2, \ldots, v_d)$ across the capacitor of each cell is given by
\begin{equation}d\mathbf{v} = \mathbf{C}^{-1}(-\mathbf{J}\mathbf{v}dt+\mathbf{R}^{-1}d\mathbf{w})\end{equation}
where $\mathbf{C} = \mathrm{diag}(C_1, C_2, \ldots, C_d)$, $\mathbf{R} = \mathrm{diag}(R_1, R_2, \ldots, R_d)$, $\mathbf{w}$ is uncorrelated Brownian motion and the elements of $\mathbf{J}$ are given by
\begin{equation}
\begin{cases}
J_{ij} = -\frac{1}{R_{ij}} \text{ if } i\neq j \\
J_{ij} = \frac{1}{R_{i}} + \frac{1}{R'_{i}} + \frac{1}{R_{ij}} \text{ if } i= j.
\end{cases}
\end{equation}
Hence there are two in-cell resistors, allowing for freedom on the diagonal elements of the $\mathbf{J}$ matrix independently of the $\mathbf{R}$ matrix, and one resistor coupling each cell. One can set $\mathcal{J} = \mathbf{J} \mathbf{R}$, and $\gamma = 1/RC$ as $\mathbf{R},\mathbf{C}$ are diagonal matrices. By denoting $x = \mathbf{v}$ and choosing $\mathcal{J} = A$, this leads to:

\begin{equation}
dx =  \frac{1}{\gamma}(-Axdt + d\mathbf{w})
\end{equation}
which reduces to overdamped dynamics of the linear systems solver provided the mean of the noise is $\gamma b$, and its variance is $2\gamma/\beta$. For more details on this implementation, see ref.~\cite{coles2023thermodynamic}.

Note that to implement this in hardware the values of the in-cell and out-of-cell resistors must be computed based on the $A$ matrix, therefore this incurs a cost of $O(d^2)$ operations at the initialization. For electrical hardware, the effective temperature is related to Johnson-Nyquist noise as $\overline{v_n^2} = \sqrt{4k_BTR\Delta f}$, with $\Delta f$ the bandwidth of the system. For a bandwidth of $1 \text{MHz}$, we obtain $\Bar{v_n^2} = 0.1 \mu V$ for $T = 300K$.

\section{Analysis of the Underdamped Regime}
\label{app:UDL}

\subsection{Stationary Distribution }

Consider the following underdamped Langevin equations
\begin{equation}
\label{UDL-x-app}
dx  = \frac{1}{M} p \, dt,
\end{equation}
\begin{equation}
\label{UDL-p-app}
dp = -(A x - b) \, dt- \frac{\gamma}{M} p \, dt + \mathcal{N}[0,2 \gamma \beta^{-1}  \mathbb{I} dt],
\end{equation}
where $M, \gamma, \beta \in \mathbb{R}^+$. As before, define $y = x - A^{-1} b$, resulting in the transformed equations
\begin{equation}
\label{UDL-x-app-transformed}
dy  = \frac{1}{M} p \, dt,
\end{equation}
\begin{equation}
\label{UDL-p-app-transformed}
dp = -A y \, dt- \frac{\gamma}{M}p \, dt + \mathcal{N}[0,2 \gamma \beta^{-1} \mathbb{I} dt].
\end{equation}
Equations \eqref{UDL-x-app-transformed} and \eqref{UDL-p-app-transformed} can be made dimensionless by defining $\tilde{x}$ and $\tilde{p}$ as follows
\begin{equation}
\tilde{y} =\sqrt{\frac{\gamma^2 \beta}{M}}y, \: \: \: \tilde{p} = \sqrt{\frac{\beta}{M}}p.
\end{equation}
and note that these choices imply that $\tilde{p} = M \gamma^{-1} \dot{\tilde{x}}$. The underdamped Langevin equations can now be formulated as a single dimensionless equation by concatenating the vectors $\tilde{y}$ and $\tilde{p}$ into a single vector $r$,
\begin{equation}
\label{r-def-app}
r = \left(\tilde{y}_1, \dots \tilde{y}_d, \dots \tilde{p}_1, \dots \tilde{p}_d\right)^\intercal,
\end{equation}
As in the Methods section, we write the UDL equations as a single matrix equation
\begin{equation}
dr = -\mathcal{A} r\, dt + \mathcal{N}[0, \mathcal{B} dt],
\end{equation}
\begin{equation}
\mathcal{A} = \begin{pmatrix}
0 & -\gamma M^{-1}\mathbb{I} \\
\gamma^{-1}A & \gamma M^{-1} \mathbb{I}
\end{pmatrix}, \: \: \: \:
\mathcal{B} = \begin{pmatrix}
0 & 0  \\
0 & \beta M^{-1} B
\end{pmatrix}.
\end{equation}
Following the analysis from the Methods section further, we find that the stationary distribution, in the original coordinates, is
\begin{equation}
x \sim \mathcal{N}[A^{-1}b, \beta^{-1}A^{-1}], \: \: p \sim  \mathcal{N}\left[0, \beta^{-1}M \mathbb{I}\right].
\end{equation}
If instead we had $B = 2\gamma\beta^{-1}R$ for an SPD matrix $R$, the stationary covariance matrix for $x$ would be the solution to the Lyapunov equation
\begin{equation}
A \Sigma_\text{s} +\Sigma_\text{s} A =2\beta^{-1}R,
\end{equation}
and the mean would still be $A^{-1}b$.

\subsection{Transient Behavior}
We review some results about the transient behavior of a system of coupled damped harmonic oscillators. Consider Eqs. \eqref{UDL-x-app-transformed} and \eqref{UDL-p-app-transformed} with the noise term removed. These can be combined into a single second order equation
\begin{equation}
M   \ddot{y} + \gamma \dot{y} + A y =0.
\end{equation}
The system may be decoupled into $d$ independent equations by expanding $y$ in the (orthonormal) eigenbasis of $A$,
\begin{equation}
y = \sum_j c_j a_j,
\end{equation}
where $a_j$ are the eigenvectors of $A$ with eigenvalues $\alpha_j$. Now we have the set of uncoupled equations
\begin{equation}
\label{decoupled-DHO-app}
\ddot{c}_j + \frac{\gamma}{M} \dot{c}_j+ \frac{\alpha_j}{M} c_j =0.
\end{equation}
We define $\omega_j = \sqrt{\alpha_j/M}$, $\xi = \gamma/2M$, and $\zeta_j = \xi/\omega_j $. The general solution to Eq. \eqref{decoupled-DHO-app} is
\begin{equation}
\label{DHO-solution-app}
c_j(t) = c_{j+} e^{\lambda_{j+} t} + c_{j-}e^{\lambda_{j-}t},
\end{equation}
\begin{equation}
\label{DHO-eigenvalues-app}
\lambda_{j\pm}=  - \xi \pm \sqrt{\xi^2 - \omega_j^2}.
\end{equation}
If $\xi < \omega_j$, the square root in Eq. \eqref{DHO-eigenvalues-app} is imaginary, in which case the system is said to be underdamped, and when $\xi > \omega_j$ the system is overdamped. For $\xi = \omega_j$ Eq. \eqref{DHO-solution-app} is invalid, which is the critically damped case. In what follows, we assume that $\xi < \omega_j$. In this case, the solution is
\begin{equation}
c_j(t) = c_{j0}e^{-\xi t}\sin(\omega_j t + \phi_j),
\end{equation}
We define the energy functions
\begin{equation}
E_j(t) = \frac{1}{2} \alpha_j c_j^2 + \frac{1}{2} M \dot{c}_j^2.
\end{equation}
It can be shown that the energy is
\begin{equation}
E_j(t) =  \frac{1}{2} M c_{j0}^2\omega_j^2e^{-2 \xi t} \left(1 + \frac{\xi  }{\omega_j}\cos(2 \omega_j t + \phi_j) \right),
\end{equation}
which can be bounded as
\begin{equation}
\frac{1}{2} M c_{j0}^2\omega_j^2e^{-2 \xi t} \left(1 - \frac{\xi  }{\omega_j} \right) \leq  E_j(t) \leq  \frac{1}{2} M c_{j0}^2\omega_j^2e^{-2 \xi t} \left(1 + \frac{\xi  }{\omega_j} \right).
\end{equation}
This implies that the energy decays exponentially,
\begin{equation}
E_j(t) \leq \frac{1+\xi/\omega_j}{1-\xi/\omega_j}e^{-2 \xi t} E_j(0).
\end{equation}
Therefore the total energy, which is the sum of the mode energies $E_j$, decays at least as fast as the mode whose energy decays slowest,
\begin{equation}
E(t) = \frac{1}{2} y^\intercal A y +\frac{1}{2M} p^2 \leq  \frac{1+\xi/\omega_\text{min}}{1-\xi/\omega_\text{min}}e^{-2 \xi t} E(0)
\end{equation}
When regarded as a function of the state vector $r\in \mathbb{R}^{2d}$, the square root of the energy is a norm, which we denote by $\|r\|_E = \sqrt{\beta E(r)}$.
\begin{align}
\|r\|_E^2 &= \frac{\beta}{2}y^\intercal A y + \frac{\beta}{2M}p^\intercal p\\
&= \frac{M}{2\gamma^2}\tilde{y}^\intercal A \tilde{y} + \frac{1}{2 } \tilde{p}^\intercal \tilde{p}
\end{align}
Note that $y^\intercal A y \geq \alpha_\text{min} \|y\|^2$, so the energy norm gives the following upper bound on the Euclidean norm of the vector $y$
\begin{equation}
\|y\| \leq  \sqrt{\frac{y^\intercal A y}{\alpha_\text{min}}} \leq \sqrt{\frac{2}{\beta \alpha_\text{min}}} \left \|r \right\|_E.
\end{equation}
We may also consider the matrix norm induced by the energy norm, which we also denote $\| \cdot\|_E$
\begin{equation}
\|O\|_E = \max_{r\in \mathbb{R}^{2d}} \sqrt{\frac{\beta E(Or)}{\beta E(r)}}.
\end{equation}
It is instructive to look at the upper left $d \times d$ submatrix of $O$, which we call $O_{y}$. In particular, observe that
\begin{align}
\|O\|_E &= \max_{r\in \mathbb{R}^{2d}} \sqrt{\frac{\beta E(Or)}{\beta E(r)}} \\
& \geq \max_{y \in \mathbb{R}^{d}}\sqrt{\frac{(O_{y}y)^\intercal A (O_{y}y)}{y^\intercal A y}} \\
& \geq
\max_{y \in \mathbb{R}^{d}}\sqrt{\frac{\alpha_\text{min}(O_{y}y)^\intercal (O_{y}y)}{\alpha_\text{max}y^\intercal  y}} \\
& = \kappa^{-1/2} \|O_{y}\|,
\end{align}
\begin{equation}
\|O_y\| \leq \sqrt{\kappa}\|O\|_E,
\end{equation}
a fact which will be used later. The earlier derivation shows that in the underdamped limit,
\begin{equation}
E\left(e^{-\mathcal{A} t} r \right) \leq \frac{1+\xi/\omega_\text{min}}{1-\xi/\omega_\text{min}} e^{-2 \xi t} E(r),
\end{equation}
and so the operator exponential $e^{-\mathcal{A} t}$ is bounded in the energy norm
\begin{align}
\left\|e^{-\mathcal{A} t}\right\|_E\leq  \sqrt{\frac{1+\xi/\omega_\text{min}}{1-\xi/\omega_\text{min}}}e^{- \xi t} \equiv \chi e^{-\xi t},
\end{align}
where we have defined
\begin{equation}
\chi = \sqrt{\frac{1+\xi/\omega_\text{min}}{1-\xi/\omega_\text{min}}}.
\end{equation}

\subsection{Equilibration and Correlation Time}
We proceed in much the same way as in the overdamped case. We now have a bound on the energy norm of $e^{- \mathcal{A} t}$,
\begin{align}
\left\|e^{-\mathcal{A} t}\right\|_E\leq  \chi e^{- \xi t}.
\end{align}
The mean of the distribution is governed by
\begin{equation}
\braket{r(t)} = e^{-\mathcal{A} t}\braket{r(0)}.
\end{equation}
We get a bound on the Euclidean norm of $\braket{y}$
\begin{align}
\left\|\braket{y(t)}\right\| & \leq \sqrt{\frac{2}{\beta \alpha_\text{min}}} \left\|\braket{r(t)}\right\|_E \\
& \leq \sqrt{\frac{2}{\beta\alpha_\text{min}}}
\chi e^{- \xi t} \left\|\braket{r(0)} \right\|_E \\
& =
\sqrt{\frac{2E(\braket{r(0)})}{\alpha_\text{min}}}\chi e^{-\xi t}\\
& = \sqrt{\frac{2\kappa E(\braket{r(0)})}{\|A\|}}\chi e^{-\xi t},
\end{align}
where $E_0$ is the energy of the state $\braket{r_0}$. Define the relative error tolerance as $\varepsilon_{\mu 0} > \|A x - b\|/\|b\| = \|A y\|/\|b\|$. Also assume that $\braket{x(0)} = \braket{p(0)} = 0$, so $E(\braket{r(0)}) = \frac{1}{2} b^\intercal A^{-1} b$. We then see that
\begin{align}
\frac{\|Ay\|}{\|b\|} & \leq\frac{ \|A\|}{\|b\|} \sqrt{\frac{2\kappa E(\braket{r(0)})}{\|A\|}}\chi e^{-\xi t}\\
&=  \sqrt{\frac{2\|A\|\kappa E(\braket{r(0)})}{\|b\|^2}}\chi e^{-\xi t}\\
& \leq  \sqrt{\kappa}\chi e^{-\xi t}\\
\end{align}
In order to ensure that the relative error is less than $\varepsilon_{\mu 0}$ , we require
\begin{equation}
t_0 \geq \frac{1}{\xi} \ln \left( \frac{\sqrt{\kappa} \chi}{\varepsilon_{\mu 0}}\right).
\end{equation}
Defining the underdamped relaxation time $\tau_\text{r(UD)}=1/\xi$, we have
which may be written
\begin{equation}
t_0 \geq \tau_{\text{r(UD)}} \ln \left( \sqrt{\kappa}\chi \varepsilon_{\mu 0}^{-1}\right).
\end{equation}
As in the overdamped case, we proceed to determine the convergence rate of the (dimensionless) covariance matrix. The covariance matrix for $r$ obeys
\begin{align}
\left\|\Sigma_\text{s} - \Sigma(t_0) \right\|_E
& \leq
\left\| e^{-\mathcal{A} t}\right\|^2 \cdot \left\| \Sigma_0 - \Sigma_\text{s}\right\|_E+ \left\|\mathcal{B}\right\|_E\int_{t_0}^\infty dt' \left\| e^{-\mathcal{A} t'}\right\|_E^2 \\
& \leq    \chi^2 e^{- 2\xi t_0 }\left\| \Sigma_0 - \Sigma_\text{s}\right\|_E+ \left\|\mathcal{B}\right\|_E\int_{t_0}^\infty dt'\chi^2 e^{-2\xi t'} \\
& =  \chi^2 e^{- 2\xi t_0 }\left\| \Sigma_0 - \Sigma_\text{s}\right\|_E+\frac{\chi^2}{2\xi}e^{-2 \xi t_0} \left\|\mathcal{B}\right\|_E \\
& = \chi^2 e^{- 2\xi t_0 }\left(\left\| \Sigma_0 - \Sigma_\text{s}\right\|_E+\frac{\gamma }{M \xi}\right)\\
& =  \chi^2 e^{- 2\xi t_0 }\left(\left\| \Sigma_0 - \Sigma_\text{s}\right\|_E+2\right)
\end{align}
Now we assume the initial covariance matrix in the dimensionless coordinates is $\Sigma_{0,r} = \gamma^2 M^{-1} \|A^{-1}\| \mathbb{I}\oplus \mathbb{I}$, which corresponds to a potential $V_0(x) = \|A\| x^\intercal x$. Therefore, $\Sigma_{\text{s},r} - \Sigma_{0,r} = \gamma^2 M^{-1}( A^{-1} - \|A^{-1}\|\mathbb{I}) \oplus 0$, which implies
\begin{align}
\|\Sigma_{\text{s},r} - \Sigma_{0,r}\|_E &\leq \|\gamma^2 M^{-1} A^{-1}\oplus 0\|_E \\
& =\gamma^2 M^{-1} \max_{v} \sqrt{\frac{v^\intercal A^{-1}AA^{-1} v}{v^\intercal A v}} \\
& \leq \gamma^2 M^{-1} \|A\|^{-1} \kappa.
\end{align}
This implies that the spectral norm error of the covariance matrix for $\tilde{y}$ is bounded by
\begin{align}
\left\|\Sigma_{\text{s},\tilde{y}} - \Sigma_{\tilde{y}}(t_0)\right\| &\leq
\sqrt{\kappa}\chi^2 e^{- 2\xi t_0 }\left(\left\| \Sigma_0 - \Sigma_\text{s}\right\|_E+2\right)\\
& \leq \sqrt{\kappa}\chi^2 e^{- 2\xi t_0 }\left(\gamma^2 M^{-1} \|A\|^{-1} \kappa+2\right)
\end{align}
which gives the corresponding result for the covariance matrix of $y$
\begin{equation}
\left\|\Sigma_{\text{s},y} - \Sigma_{y}(t_0)\right\| \leq \beta^{-1}
\sqrt{\kappa}\chi^2 e^{- 2\xi t_0 }\left(\kappa\|A\|^{-1}+M \gamma^{-2}\right).
\end{equation}
We want the relative error to be less than tolerance $\varepsilon_{\Sigma0}$
\begin{equation}
\frac{\left\|\Sigma_{\text{s},y} - \Sigma_{y}(t_0)\right\|}{\beta^{-1} \|A^{-1}\|}\leq \varepsilon_{\Sigma0},
\end{equation}
which is satisfied if we require
\begin{equation}
t_0 \geq \frac{1}{2} \tau_\text{r(UD)}\ln\left(\chi^2 \varepsilon_{\Sigma0}^{-1}\left[\kappa^{3/2}+\kappa^{1/2}M \gamma^{-2} \|A\|\right]\right),
\end{equation}
\begin{equation}
t_0 \geq \frac{1}{2} \tau_\text{r(UD)}\ln\left(\chi^2 \kappa^{3/2}\varepsilon_{\Sigma0}^{-1}\left[\frac{1}{4\zeta_\text{max}^{2}}+1\right]\right).
\end{equation}
Combining the equilibration times for mean and variance gives
\begin{equation}
t_0 \geq \max \left\{ \tau_{\text{r(UD)}} \ln \left( \kappa^{1/2}\chi \varepsilon_{\mu 0}^{-1}\right),\frac{1}{2} \tau_\text{r(UD)}\ln\left(\chi^2 \kappa^{3/2}\varepsilon_{\Sigma0}^{-1}\left[\frac{1}{4\zeta_\text{max}^{2}}+1\right]\right)\right\}.
\end{equation}

\subsection{Ergodicity of Mean}
Once equilibrium has been reached, we may consider the correlation matrix of the stationary system
\begin{equation}
G_\text{s}(\tau) = \lim_{t\to\infty}G(t,t+\tau),
\end{equation}
which is given by
\begin{equation}
G_\text{s}(\tau) = e^{-\mathcal{A} \tau}\Sigma_\text{s}.
\end{equation}
We define the $A$-norm of a vector as $\|v\|_A ^2 = v^\intercal A v$, which induces an $A$-norm on operators. It clearly follows from the definitions that for all vectors $r = (u, v)^\intercal$, $\|u\|_A \leq \|r\|_E$, and similarly if $X_{y}$ is the upper left $d \times 2$ submatrix of a $2  d \times 2d$ matrix $X$, we have $\|X_{y}\|_A \leq \|X\|_E$.
\begin{align}
\|G_{\text{s},\tilde{y}}(\tau)\|_A &\leq \|G_\text{s}(\tau)\|_E\\
& \leq \|e^{-\mathcal{A} t}\|_E \cdot \| \Sigma_\text{s}\|_E \\
& \leq  \chi e^{-\xi \tau} \|\Sigma_\text{s}\|_E.
\end{align}
The quantity $\|\Sigma_\text{s}\|_E$ can be evaluated by the following arguments: after changing to the dimensionless coordinates, the covariance matrix is
\begin{equation}
\Sigma_\text{s} =
\begin{pmatrix}
\gamma^2 M^{-1} A^{-1} & 0\\
0 & \mathbb{I}
\end{pmatrix}.
\end{equation}
By definition,
\begin{align}
\|\Sigma_\text{s}\|_E^2 &= \max_{y,p} \frac{ M \gamma^{-2}y^\intercal \Sigma_{\text{s},\tilde{y}} A \Sigma_{\text{s},\tilde{y}} y  + p^\intercal p}{M \gamma^{-2} y^\intercal A y + \ p^\intercal p} \\
& =
\max_{y,p} \frac{ y^\intercal  A^{-1}y  + p^\intercal p}{M^2 \gamma^{-4} y^\intercal A y + \ p^\intercal p} \\
& =
\max\left\{1, M^{-2}\gamma^{4}\|A^{-1}\|_A^2\right\}.
\end{align}
We then see that
\begin{equation}
\|\Sigma_\text{s}\|_E = \max\{1, \gamma^{2} M^{-1} \alpha_\text{min}^{-1}\}.
\end{equation}
By the assumption that the system is in the underdamped regime, $2M \gamma^{-2}\alpha_\text{min}^{-1} >1$, so we finally obtain
\begin{equation}
\|\Sigma_\text{s}\|_E = \gamma^{2}M^{-1}\alpha_\text{min}^{-1} =\gamma^{2}M^{-1}\kappa \|A\|^{-1}.
\end{equation}
By the same manipulations as in the overdamped case, we arrive at
\begin{align}
\left\|\braket{\bar{y}\bar{y}^\intercal}\right\| &\leq
\frac{2 M\sqrt{\kappa}}{\tau^2\gamma^2 \beta}\int_{t_0}^{t_0 + \tau}dt''\int_{t_0}^{t''} dt' \left\|G_{\text{s}}(t''-t') \right\|_E\\
\frac{2M\sqrt{\kappa}}{\tau^2\gamma^2 \beta}\chi  \left\|\Sigma_\text{s}\right\|_E\int_{t_0}^{t_0 + \tau}dt''\int_{t_0}^{t''} dt' e^{-\xi(t''-t')}\\
\frac{2\sqrt{\kappa}}{\xi\tau^2\beta\|A\|} \chi  \int_{t_0}^{t_0 + \tau}dt'' 1-e^{-\xi(t'' - t_0)} \\
& \leq
\frac{2 \sqrt{\kappa}}{\xi\tau^2 \beta \|A\|} \chi  \int_{t_0}^{t_0 + \tau}dt''1 \\
\frac{2\sqrt{\kappa} \chi \tau_{\text{r(UD)}}}{ \beta \|A\|\tau}
\end{align}
The exact same reasoning employed in the overdamped case now gives the following statement: In order to have
\begin{equation}
\left \|A \bar{x} - b \right\| \leq \varepsilon_x \|b\|
\end{equation}
with probability at least $P_\varepsilon$, it is sufficient to require that
\begin{align}
\tau \geq  \frac{2\sqrt{\kappa} \chi \|A\| d}{ \beta \|b\|^2\varepsilon_x^2(1-P_\varepsilon)} \tau_{\text{r(UD)}}.
\end{align}
We note that by the same reasoning as was applied in the overdamped case, we may state the above as an energy-time tradeoff
\begin{align}
\mathcal{E}\tau \geq  \frac{2\sqrt{\kappa} \chi  d}{  \varepsilon_x^2(1-P_\varepsilon)} \beta^{-1}\tau_{\text{r(UD)}}.
\end{align}

\subsection{Ergodicity of Covariance Matrix}
Repeating the steps from before, we evaluate
\begin{align}
\sum_{ij}\text{var}(\overline{y_i y_j})
&\leq \frac{2d(d+1)}{\tau^2}\int_{t_0}^{t_0 + \tau} dt'' \int_{t_0}^{t_0 + \tau} dt'\|G_{\text{s},y}(t'' - t')\|^2\\
& =  \frac{4d(d+1)M^2\kappa}{\gamma^4 \beta^2 \tau^2}\int_{t_0}^{t_0 + \tau} dt'' \int_{t_0}^{t''} dt'\|G_\text{s}(t'' - t')\|^2 \\
& \leq \frac{4d(d+1)M^2\kappa}{\gamma^4 \beta^2 \tau^2} \int_{t_0}^{t_0 + \tau} dt'' \int_{t_0}^{t''} dt'e^{-2\xi(t''-t')} \left\| \Sigma_\text{s}\right\|_E^2 \\
\frac{4d(d+1)\kappa^3}{\beta^2\|A\|^2 \tau^2 } \int_{t_0}^{t_0 + \tau} dt'' \int_{t_0}^{t''} dt'e^{-2\xi(t''-t')}  \\
\frac{4d(d+1)\kappa^3}{\beta^2  \|A\|^2\tau^2 \xi} \int_{t_0}^{t_0 + \tau}dt'' 1 - e^{-2\xi(t''-t_0)}\\
& \leq \frac{4d(d+1)\kappa^3}{\beta^2  \|A\|^2\tau \xi} \\
& = \frac{4d(d+1)\kappa \|\Sigma_{\text{s},y}\|^2\tau_{\text{r(UD)}}}{\tau }
\end{align}

Then repeating the same steps as in the overdamped, we find that in order to ensure that $\|\overline{y y^\intercal} - \Sigma_\text{s}\|_F/\|\Sigma_\text{s}\|_F \leq \varepsilon_\Sigma$, it is sufficient to require that
\begin{equation}
\tau \geq \frac{4 \kappa d(d+1)}{(1-P_\varepsilon)\varepsilon_\Sigma^2}\tau_\text{r(UD)}.
\end{equation}

\subsection{Asymptotic scaling with condition  number}
The results derived for an underdamped system are only valid insofar as the system is really in the underdamped regime, i.e. when $\xi < \omega_\text{min}$, where $\xi = \gamma/ 2M$. However, if the parameters of the system are fixed,  then as $\kappa$ is increased (and the smallest eigenvalue $\alpha_\text{min}$ decreases), the system is eventually not in the underdamped regime anymore, and would begin to behave as an overdamped system. In this sense, the above results for an underdamped system cannot be interpreted as precisely ``asymptotic" in $\kappa$, as they are invalidated when $\kappa$ tends to infinity. Using the definitions $\omega_j = \sqrt{\alpha_j/M}$, and $\kappa = \|A\|/\alpha_\text{min}$, we see that the underdamped treatment is valid when
\begin{equation}
\frac{\gamma}{2M}< \sqrt{\frac{\|A\|}{\kappa M}},
\end{equation}
\begin{equation}
M > \frac{\gamma^2 \kappa}{4\|A\|}
\end{equation}
Then, because $\tau_\text{r(UD)} = 2 M / \gamma$, we have
\begin{equation}
\tau_\text{r(UD)}>   \frac{\gamma \kappa}{2\|A\|}.
\end{equation}
Therefore, as the averaging time $\tau_\text{r(UD)}$ is proportional to $\gamma \kappa$, we include a factor of $\kappa \gamma$ in the expressions for asymptotic time-complexity in the underdamped case.
\section{Determinant Estimation}
Suppose we are given a series of work measurements $W_1 \dots W_N$ resulting from changing a potential from $U_1$ to $U_2$. Importantly, the system must be allowed to come to equilibrium once per measurement (so $N$ times in total). It has been shown, via the Cramér-Rao bound \cite{cramer1999mathematical}, that if $\Delta \hat{F}$ is \emph{any} asymptotically unbiased estimator of the difference in equilibrium free energies between the two potentials, then its variance is bounded below \cite{shenfeld2009minimizing}
\begin{equation}
\text{var}\left(\Delta \hat{F} \right)\geq \frac{\mathcal{L}^2}{N},
\end{equation}
where $\mathcal{L}$ is the Fisher-Rao length of the path in parameter space between the two equilibrium states. Explicitly, the Fisher-Rao length of a path through the space of probability distributions is
\begin{equation}
\mathcal{L} = \int_0^\tau dt \sqrt{\int dx \frac{\dot{f}_t(x)^2}{f_t(x)}}.
\end{equation}
For Gaussian distributions sharing the same mean, an explicit formula has been found for the minimal Fisher-Rao length between two distributions. In particular, if one of the Gaussian distributions has covariance matrix $\beta^{-1}A^{-1}$ and the other has covariance matrix $\beta^{-1} \|A^{-1}\|\mathbb{I}$, and both share the same mean, then the minimal Fisher-Rao distance is given by \cite{nielsen2023simple}
\begin{equation}
\mathcal{L}_{\text{min}}  = \sqrt{\frac{1}{2} \sum_{j=1} \ln( \lambda_j(A\|A\|^{-1}))^2},
\end{equation}
which is bounded by
\begin{align}
\mathcal{L}_{\text{min}} &\leq  \sqrt{\frac{1}{2} d \left|\ln( \kappa^{-1})\right|^2} \\
&\leq  \sqrt{\frac{1}{2} d \ln( \kappa)^2} \\
\end{align}
Therefore the minimal variance of an unbiased estimator is
\begin{equation}
\text{var}\left(\Delta \hat{F} \right)\geq \frac{d\,\ln(\kappa)^2}{2N \beta^2},
\end{equation}
Moreover, this lower bound can be attained in the limit of large $N$ by the Bennett Acceptance Ratio (BAR) estimator or its multistate generalization (MBAR). The free energy difference between equilibrium states of the device is
\begin{equation}
\Delta F = F_2 - F_1 = \frac{1}{2}(b_1^\intercal A_1^{-1}b_1-b_2^\intercal A_2^{-1}b_2) + \frac{1}{2\beta} \ln \left( \frac{\left|A_2\right|}{\left|A_1\right|}\right).
\end{equation}
Apparently, if $b_1 = b_2 = 0$ and $A_1 = a_1 \mathbb{I} $ then
\begin{equation}
\ln\left(\left|A_2\right|\right)= 2\beta \Delta F +  d\, \ln(a_1).
\end{equation}
Therefore the variance of an unbiased estimator of the log determinant $\ln(\left|A_2\right|)$ is lower bounded as follows
\begin{equation}
\text{var}(\widehat{\text{LD}}(A_2)) \geq \frac{2 d\,\ln(\kappa)^2}{N}
\end{equation}
where we write $\text{LD}$ for the log determinant and $\widehat{\text{LD}}$ for the estimator of the log determinant. To achieve absolute error of at most $\delta_\text{LD}$ with probability $P_\delta$, according to Chebyshev's inequality we need a number of samples given by
\begin{equation}
\text{var}(\widehat{\text{LD}}(A_2)) \leq \delta_\text{LD}^2 (1-P_\delta),
\end{equation}
which can be guaranteed (assuming the estimator is optimal) by requiring
\begin{equation}
N> \frac{2 d\, \ln(\kappa)^2}{\delta^2_\text{LD}(1-P_\delta)}.
\end{equation}
As the equilibration does not require a change of the mean of the distribution but just the covariance matrix, the equilibration time is, in the overdamped case
\begin{equation}
\frac{1}{2} \kappa \tau_\text{r}\ln\left(2\kappa\varepsilon_{\Sigma 0}^{-1}\right),
\end{equation}
for some sufficiently small $\varepsilon_{\Sigma 0}$. The total time of the optimal protocol would then be approximated by
\begin{equation}
\tau \approx \frac{ d\, \kappa \ln(\kappa)^2\ln(2 \kappa\varepsilon^{-1}_{\Sigma0})}{\delta^2_\text{LD}(1-P_\delta)} \tau_\text{r} = O(d \, \kappa \ln(\kappa)^3).
\end{equation}
In the underdamped case, the equilibration time would be
\begin{equation}
\frac{1}{2} \tau_\text{r(UD)}\ln\left(\chi^2 \kappa^{3/2}\varepsilon_{\Sigma0}^{-1}\left[\frac{1}{4\zeta_\text{max}^{2}}+1\right]\right),
\end{equation}
for a total time of
\begin{equation}
\tau \approx \frac{ d\,\ln(\kappa)^2}{\delta^2_\text{LD}(1-P_\delta)} \ln\left(\chi^2 \kappa^{3/2}\varepsilon_{\Sigma0}^{-1}\left[\frac{1}{4\zeta_\text{max}^{2}}+1\right]\right)
\tau_\text{r(UD)} = O(d \, \ln(\kappa)^3).
\end{equation}

\begin{figure}
\centering
\includegraphics[width=0.5\linewidth]{figures/rel_det_error_w_samples_F.pdf}
\caption{\textbf{Convergence of determinant estimation using the Jarzynski estimator.} This shows relative error of the determinant estimate with increasing numbers of samples used in the Jarzynski estimator. Solid lines show the median and the error bars show the interquartile range over 10 different $A$ matrix realizations.}
\label{fig:det}
\end{figure}

\subsection{Numerical Results}\label{app:det_numerics}

To explore the performance of our thermodynamic algorithm for estimating the log determinant of a matrix we numerically explore using a Jarzynski estimator of the free energy difference in the simplest case. Here the potential is initialized using an identity matrix and this is instantaneously switched to the matrix for which we desire to evaluate the log determinant. We explore the performance of the estimated log determinant when using a Jarzynski estimator with different numbers of samples. Recall the free energy can be estimated using
\begin{equation}
e^{-\beta \Delta F}\approx \overline{e^{-\beta W}}  \equiv \frac{1}{N}\sum_{j=1}^N e^{-\beta W_j}.
\end{equation}
This in turn is then used to estimate $\ln\left( \left| A_2\right|\right)$ with the relation
\begin{equation}
\ln\left(\left| A_1\right|\right) \approx 2 \ln\left(\overline{e^{-\beta W}}\right)  + \ln\left( \left| A_2\right|\right).
\end{equation}
In Fig.~\ref{fig:det} we show the relative error $\epsilon_{\log(|A|), \widetilde{\log(|A|)}}$ of the log determinant obtained for several dimensions, estimated with increasing numbers of samples. This shows estimating the determinant in this fashion is a valid protocol.

\section{Thermodynamic Algorithm for Solving the Lyapunov Equation}

In this section we present our general algorithm for solving the Lyapunov Equation. Note that our matrix inversion algorithm from the main text is a special case of our general Lyapunov Equation algorithm.

In what follows, we assume we have access to a device with a controllable noise source such that the covariance matrix of the noise term may be chosen to be an arbitrary symmetric positive definite matrix. We do not include the linear $b^\intercal x$ term in the potential, and therefore obtain the following overdamped Langevin equation
\begin{equation}
\label{eq:odl-lyapunov}
dx = -\frac{1}{\gamma}A x dt + \mathcal{N}\left[0, \frac{2}{\gamma \beta } R\, dt\right],
\end{equation}
where $R$ is symmetric and positive definite.
In this case, the stationary distribution has mean zero and covariance matrix $\Sigma_\text{s}$, which is a solution to the Lyapunov equation
\begin{equation}\label{eq:Lyapunov}
A\Sigma_\text{s}+ \Sigma_\text{s} A^{\intercal} = 2\beta^{-1}R.
\end{equation}
We propose the following protocol for solving the Lyapunov equation.
\medskip
\begin{tcolorbox}[title={Lyapunov Equation Protocol},breakable]
\begin{enumerate}
\item Given two symmetric positive definite matrices $A$ and $R$, set the potential of the device to
\begin{equation}
U(x) = \frac{1}{2} x^\intercal A x,
\end{equation}
and the noise term in the overdamped Langevin equation to $\mathcal{N}\left[0,2\gamma^{-1}\beta^{-1} R dt\right]$ at time $t=0$. That is, the system evolves under the dynamics of Eq. \eqref{eq:odl-lyapunov}.

\item Choose equilibration tolerance parameter $\varepsilon_{\Sigma 0} \in \mathbb{R}^+$, and choose the equilibration time
\begin{equation}
\label{eq:UDL-eq-time}
t_0 \geq \widehat{t}_0,
\end{equation}
where $\widehat{t}_0$ is computed from the system's physical properties, Eq.~\eqref{eq:t0-tau-lyapunov-odl}. Allow the system to evolve under its dynamics until $t=t_0$, which ensures that $\left\|\Sigma - \beta^{-1} A^{-1} b\right\|/\|\beta^{-1}A^{-1}\| \leq \varepsilon_\Sigma$.
\item  Choose error tolerance parameter $\delta_\Sigma$ and success probability $P_\delta$, and choose the integration time
\begin{equation}
\tau \geq \widehat{\tau},
\end{equation}
where $\widehat{\tau}$ is computed from the system's physical properties, Eq.~\eqref{eq:t0-tau-lyapunov-odl}. Use analog multipliers and integrators to measure the the time averages
\begin{equation}
\overline{x_i x_j} = \frac{1}{\tau^2} \int_{t_0}^\tau dt\, x_i(t) x_j(t),
\end{equation}
which satisfies $\|\overline{x x^\intercal} - \Sigma_\text{s}\|_F \leq \delta_\Sigma$ with probability at least $P_\delta$.
\end{enumerate}
\end{tcolorbox}
\medskip

The timing parameters for the Lyapunov equation protocol are, for the overdamped case,
\begin{align}
\label{eq:t0-tau-lyapunov-odl}
\widehat{t}_0 =
\frac{1}{2}\kappa\tau_\text{r}\ln\left(\frac{\|\Sigma_0 - \Sigma_\text{s}\| + \kappa \beta^{-1} \|A\|^{-1}
}{\delta_{\Sigma 0}}\right),\quad
\widehat{\tau} = \frac{4 \kappa d(d+1)}{(1-P_\varepsilon)\varepsilon_\Sigma^2}\tau_\text{r},
\end{align}
For the underdamped case, $\widehat{t_0}$ would be somewhat different, but $\widehat{\tau}$ would be the same, because the behavior of the equilibrium correlation function does not depend on the noise, so the same result derived for the matrix inverse protocol is applicable. Note that Eq. \eqref{eq:t0-tau-lyapunov-odl} only vaguely determines the equilibration time, as the target covariance matrix $\Sigma_\text{s}$ is not known beforehand. The corresponding equilibration time $\widehat{t_0}$ for an underdamped system could also be evaluated in principle; however, this would only result in a similarly vague expression, which is anyway not necessary to determine the asymptotic time-complexity scaling of the algorithm, so it is not pursued here. Moreover, the relative error cannot be bounded as straightforwardly as was done for the other protocols given that there is no explicit formula for the target covariance matrix. For this reason, we have used absolute error as the error tolerance in the above protocol.

\end{appendix}

\end{document}