Thermodynamic Linear Algebra

TL;DR

This work proposes thermodynamic computing hardware that leverages the equilibrium distribution of coupled harmonic oscillators to accelerate core linear-algebra primitives. By mapping problems such as solving linear systems, matrix inversion, Lyapunov equations, and determinant estimation to sampling from Gaussian distributions generated by overdamped or underdamped Langevin dynamics, the authors derive asymptotic speedups that scale linearly with dimension $d$ under plausible assumptions. They provide explicit algorithmic protocols, convergence analyses, and a timing model that compares favorably to digital methods like conjugate gradient and Cholesky for large-scale, ill-conditioned dense matrices, while highlighting practical energy-time tradeoffs and hardware considerations. The work also outlines a pathway to experimental demonstrations and broader implications for thermodynamic computing, including non-linear extensions and multi-application hardware reuse. Overall, the paper establishes a concrete mathematical foundation and numerical benchmarks for using thermodynamic hardware to accelerate linear algebra in near-term devices.

Abstract

Linear algebraic primitives are at the core of many modern algorithms in engineering, science, and machine learning. Hence, accelerating these primitives with novel computing hardware would have tremendous economic impact. Quantum computing has been proposed for this purpose, although the resource requirements are far beyond current technological capabilities, so this approach remains long-term in timescale. Here 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. In this work, 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 (1) solving linear systems of equations, (2) computing matrix inverses, (3) computing matrix determinants, and (4) solving Lyapunov equations. 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 computing hardware.

Figures (6)

• Figure 1: 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.
• Figure 2: 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 fokker1914mittlere 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.
• Figure 3: 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.
• Figure 4: 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.
• Figure 5: 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.