KR20250109682A - A thermodynamic computational system for sampling high-dimensional probability distributions - Google Patents

Abstract

A thermodynamic computing platform samples a multivariate distribution, such as a multivariate Gaussian distribution, using a network of coupled analog unit cells controlled by a digital processor. The analog unit cells are coupled to each other capacitively, resistively, and/or inductively. Each analog unit cell includes a tunable resistor, capacitor, inductor, current source, and/or voltage source. The digital processor sets initial values for these components. These values, or input parameters, represent the multivariate distribution. The coupled analog unit cell network is then evolved over a series of time steps until it reaches thermal equilibrium. While the network evolves, the digital processor samples voltages at nodes in the network at each time step. These voltages represent samples of the multivariate distribution.

Description

A thermodynamic computational system for sampling high-dimensional probability distributions

This application claims the benefit of U.S. patent application Ser. No. 63/518,545, filed Aug. 9, 2023; U.S. patent application Ser. No. 63/492,832, filed Mar. 29, 2023; and U.S. patent application Ser. No. 63/379,561, filed Oct. 14, 2022, under 35 U.S.C. 119(e). Each of these applications is incorporated herein by reference in its entirety for all purposes.

Traditional artificial intelligence (AI) and machine learning (ML) work best when uncertainty is not a concern. However, uncertainty is often rife. Noisy data adds uncertainty to AI and ML predictions, making them unreliable. This reduces the AI’s ability to judge when to follow human input. This is especially true in the context of high-stakes applications for AI, such as criminal justice, healthcare, energy, finance, defense, manufacturing, and robotics. For example, consider an assistant model that diagnoses melanoma. Ideally, such a model should defer to a medical professional when uncertain. Similarly, an AI in a self-driving car should hand over control to a human driver when it is uncertain whether a pedestrian is in the roadway. However, if uncertainty impairs the AI’s ability to predict, it may not be able to effectively follow human judgment.

Uncertainty can be divided into two main categories. First, there is ambiguity, which is usually caused by noisy data or missing features (“aleatoric ambiguity”). Aleatoric ambiguity persists regardless of how much data an AI or ML system is exposed to. For example, a perception system for a self-driving car may contain noisy sensors and imaging. Second, there is uncertainty due to finite training data (“epistemic uncertainty”). Model inputs may look different from the training data, or there may be few examples in the training data that support the inputs. Epistemic uncertainty decreases as the model is trained with more data.

Traditional AI addresses uncertainty by assigning a level of confidence to its decisions (e.g., 90% confidence that an image shows a particular person or object). By assigning a confidence level to a decision, AI quantifies the uncertainty associated with that decision. Unfortunately, traditional AI and ML generally lack meaningful ways to assign precisely quantified uncertainty.

Figure 1 illustrates how uncertainty affects the ability of a mainstream deep neural network to assign a meaningful level of confidence to its predictions in real-world situations. It shows an idealized and actual scatterplot of the uncertainty (confidence) of the deep neural network's predictions, with darker shading indicating higher confidence. The deep neural network was trained on training data (dots) distributed along the semicircles shown in the two scatterplots in Figure 1. The trained deep neural network was then trained on the y 1.5 to x 0 to x The “out-of-distribution” input data (dots) are presented, clustered up to 1.5. Ideally, the deep neural network should recognize that its predictions for the input data are very uncertain or have low confidence (light shading), as shown in the scatterplot on the left. However, in practice, the deep neural network assigns high confidence (dark shading) to its predictions, even when the input data is far from the training data, as shown in the scatterplot on the right. Unfortunately, common AI and ML models tend to be overconfident in their predictions, especially when dealing with unfamiliar input data.

Fortunately, probabilistic AI and ML can quantify and account for uncertainty in predictions. Probabilistic AI and ML combine AI and ML with probabilistic reasoning, a mathematical process that determines uncertainty as a function of prior beliefs, biases, and data. For any level of stake, probabilistic reasoning should be optimal as long as uncertainty exists. Probabilistic AI and probabilistic reasoning naturally entail continuous distributions, such as the Gaussian distribution, to reflect the user’s state of knowledge. Unfortunately, when performing probabilistic reasoning using these continuous distributions on standard digital computers, the computational difficulty becomes severely worse as the model complexity increases.

Nature performs probabilistic reasoning through the physical principles of thermodynamic equilibrium. Probabilistic AI models are typically learned and inferred through probabilistic reasoning that operates in terms of continuous random variables. A thermodynamic computing platform can perform probabilistic reasoning for probabilistic AI and machine learning using thermodynamic equilibrium. The thermodynamic equilibrium of a system describes the dynamics that evolve a physical system until it can be described by a stable probability distribution over possible states; the technology of the present invention uses this thermodynamic equilibrium for computation, and the dynamics are designed and implemented such that the final probability distribution over the states matches the target probability distribution being computed. Furthermore, utilizing a physical system that is inherently continuous naturally aligns with the mathematical framework of continuous distributions used in probabilistic reasoning.

The present invention can be implemented as a thermodynamic system for sampling a multivariate distribution, such as a multivariate Gaussian probability distribution. The system can include a network of coupled harmonic oscillators operably coupled to a digital controller. The network of coupled harmonic oscillators is characterized by a Hamiltonian that maps to the multivariate distribution. Each coupled harmonic oscillator of the network of coupled harmonic oscillators includes a digitally tunable inductor, a digitally tunable capacitor, a digitally tunable voltage source, and a digitally tunable current source. The digital controller initializes the digitally tunable inductor, the digitally tunable capacitor, the digitally tunable voltage source, and the digitally tunable current source. It repeatedly reads voltages and currents at nodes of the network of coupled harmonic oscillators and updates the digitally tunable voltage sources and the digitally tunable current sources based on the voltages and currents. These voltages and currents correspond to positions and momentums of the Hamiltonian. The digital controller also converts the voltages and currents into samples of a multivariate distribution.

A coupled harmonic oscillator network may include one harmonic oscillator for each dimension of the multivariate distribution. The harmonic oscillators may be coupled to each other capacitively, inductively, and/or resistively.

Another embodiment of the present invention is a thermodynamic system for sampling a target probability distribution. The thermodynamic system comprises both an analog dynamical system and a digital controller. The analog dynamical system is configured to evolve over a series of time steps via Hamiltonian equations to provide samples of the target probability distribution. The digital controller is operably coupled to the analog dynamical system and is configured to receive, at each of the time steps, a proposed sample of the target probability distribution and to accept or reject the proposed sample.

Another embodiment of the present technology is a thermodynamic system for sampling a multivariate distribution, such as a multivariate Gaussian probability distribution. The system can include a network of analog electrical circuits (e.g., coupled harmonic oscillators) operably coupled to a digital controller. Each of the analog electrical circuits includes at least one adjustable passive electrical component (e.g., a digitally adjustable inductor, a digitally adjustable capacitor, and/or a digitally adjustable resistor). In operation, the digital controller adjusts the adjustable passive electrical components according to the multivariate distribution and samples voltages at points within the network of analog electrical circuits.

Such thermodynamic systems can also include a voltage-controlled voltage source operably coupled to an analog electrical circuit network to modify the multivariate distribution. The voltage-controlled voltage source can employ an artificial neural network to relate an input voltage to the voltage-controlled voltage source and an output voltage of the voltage-controlled voltage source. In this case, the input voltage is based on voltages sampled by a digital controller, and the voltage-controlled voltage source can be configured to apply an output voltage to the analog electrical circuit network to modify the potential energy distribution of the entire analog electrical circuit network.

The multivariate distribution can have a non-zero cumulant of at least order 3, in which case the system comprises a multiport element operably coupled to at least three analog electrical circuits in the analog electrical circuit network. The multiport element is capable of generating k -field terms in the potential energy distribution of the analog electrical circuit network, where k is an integer greater than or equal to 3. For example, the multiport element can comprise a transistor, a voltage-controlled capacitor, or both.

Each analog electrical circuit in the thermodynamic system may also include stochastic noise sources, such as thermal noise sources, shot noise sources, or digital pseudo-random noise sources. And the analog electrical circuits in the analog electrical circuit network may be coupled to each other through capacitive and/or resistive coupling.

Thermodynamic systems can generate samples from multivariate distributions via the Langevin Monte Carlo algorithm and/or the Metropolis-adjusted Langevin algorithm.

The covariance matrix of a multivariate distribution can be encoded in at least some of the parameters of an analog electrical circuit network.

Another inventive system for sampling a multivariate distribution (e.g., a multivariate Gaussian distribution or a distribution having a non-zero cumulant higher than order 2) comprises a coupled analog unit cell network operably coupled to a digital controller. The coupled analog unit cell network is characterized by a total energy function that maps to the multivariate distribution. Each coupled analog unit cell comprises a digitally adjustable capacitor and either a digitally adjustable voltage source, a digitally adjustable current source, or both a digitally adjustable voltage source and a digitally adjustable current source. In operation, the digital controller initializes the digitally adjustable capacitor, the voltage source, and/or the current source. The digital controller repeatedly reads voltages and currents at nodes of the coupled analog unit cell network and updates the digitally adjustable voltage sources and/or the digitally adjustable current sources. The voltages and currents correspond to positions and momentums of the total energy function, and the digital controller converts the voltages and currents into samples of the multivariate distribution.

The total energy function can be Hamiltonian, in which case the coupled analog unit cell network has dynamics that are used to execute the Hamiltonian Monte Carlo protocol. Alternatively, the coupled analog unit cell network can have dynamics that are used to execute the Langevin Monte Carlo protocol.

A coupled analog unit cell network may include one analog unit cell for each dimension of the multivariate distribution. Each of the analog unit cells may also include a digitally tunable inductor. The analog unit cells may be capacitively, resistively, and/or inductively coupled to one another.

Another inventive system for sampling a multivariate distribution comprises a gradient calculator (e.g., an analog neural network or a network of resistive circuit elements), a momentum integrator device operably coupled to the gradient calculator, and a position integrator device operably coupled to the momentum integrator device.

Another inventive thermodynamic system for sampling a target probability distribution comprises an analog dynamical system operably coupled to a digital controller. The analog dynamical system is configured to evolve through the Langevin equations and/or the Hamiltonian equations of motion. And the digital controller is configured to receive a proposed sample of the target probability distribution from the analog dynamical system at each time step of the Langevin equations and/or the Hamiltonian equations of motion and to accept or reject the proposed sample.

The methods of the present invention include a method of finding an inverse matrix of a matrix. The method includes the steps of uploading elements of the matrix to respective values of tunable circuit elements of an analog unit cell network, allowing the analog unit cell network to reach thermal equilibrium, and computing a covariance matrix based on dynamic variables of the analog unit cell network. The covariance matrix is an inverse matrix of the matrix. The step of computing the covariance matrix may include the steps of extracting samples from a thermal equilibrium distribution of dynamic variables of the analog unit cell network and computing a covariance matrix of the samples. The step of computing the covariance matrix may also include the step of integrating the dynamic variables over time using a plurality of analog integrators.

Another inventive method comprises a method of solving a system of linear equations represented by a matrix and a vector. The method comprises the steps of uploading elements of the matrix to respective values of adjustable circuit elements of the analog unit cell network, uploading elements of the vector to respective current sources or voltage sources of the analog unit cells, allowing the analog unit cell network to reach thermal equilibrium, and computing averages of dynamic variables of the analog unit cell network. The averages represent solutions to the system of linear equations. The step of computing the averages may comprise the step of integrating the dynamic variables over time using a plurality of analog integrators.

All combinations of the concepts described above with additional concepts discussed in more detail below (so long as such concepts are not mutually inconsistent) are considered to be part of the inventive subject matter disclosed in this disclosure. In particular, all combinations of claimed inventive subject matter that appear at the end of this specification are considered to be part of the inventive subject matter disclosed in this disclosure. Any technical terminology explicitly employed in this disclosure that may appear in any disclosure incorporated by reference herein should have a meaning most consistent with the specific concepts disclosed in this disclosure.

Those skilled in the art will appreciate that the drawings are primarily for illustrative purposes and are not intended to limit the scope of the inventive subject matter described in this disclosure. The drawings are not necessarily to scale; in some instances, various aspects of the invention disclosed in this disclosure may be shown exaggerated or enlarged in the drawings to facilitate understanding of various features. Like reference characters in the drawings generally indicate like features (e.g., functionally similar and/or structurally similar components).
Figure 1 shows the uncertainty scatter plot of a deep neural network in ideal and real cases.
Figure 2 illustrates an overview of a thermodynamic computing system. A digital controller performs a Hamiltonian Monte Carlo (HMC) method that samples a target probability distribution by extracting samples from a time-evolving analog system, also known as an analog dynamical system.
Figure 3 is a schematic of a basic LC oscillator that maps to a one-dimensional normal (i.e., Gaussian) probability distribution.
Figure 4 is a circuit diagram of a pair of capacitively coupled oscillators in a capacitively coupled harmonic oscillator network mapped to a multivariate normal probability distribution.
Figure 5 illustrates the U-Turn that appears while running HMC on a 2D Gaussian.
Figure 6 illustrates a system level view of an example UL hardware implementation. The PC downloads setup information to the onboard controller, which in turn reads out sampled data from the cells.
Figure 7 illustrates two coupled LC resonators forming a simple two-cell system. The resonant frequency of each resonator can be tuned independently, and the capacitive coupling can be controlled. The cells are excited by a Gaussian noise current source.
Figure 8 illustrates a typical all-to-all coupling in a UL hardware setup.
Figure 9a illustrates a tunable cell whose resonant frequency can be digitally controlled using a current noise source implemented using a digital controller followed by a LPF.
Figure 9b illustrates an adjustable cell with initial conditions implemented by adding an additional switch that separates the inductor from the capacitor.
Figure 10 illustrates that instead of initializing the capacitor and inductor of each cell, the capacitor can be initialized while disconnected. At a given point in time after connecting the inductor and capacitor, the desired initial condition is achieved. In this way, one component is initialized.
Figure 11 illustrates a coupling unit between two cells. These cells are connected at nodes A and B. The coupling can be enabled or disabled, and the polarity is also indicated using a control signal.
Figure 12 illustrates the following FPGA implementation of the controller: (a) High-level control where download and upload operations occur, in addition to setting up each cell and coupling unit within the FPGA. (b) The cell interface reads the voltage from the ADC and subsequently stores the data. It also controls the switches of the tunable resonator and controls the noise amplitude. (c) The coupling unit interface controls the coupling and its polarity by turning on appropriate switches. (d) The format of the data downloaded to the FPGA.
Figure 13a shows that the current noise source of Figure 7 can be implemented by emulating a current source using a single-bit PN generator and a series resistor. The RC circuit also acts as an LPF.
Figure 13b illustrates an example of a linear feedback shift register (LFSR) used in a single-bit pseudonoise (PN) generator.
Figure 14 visualizes the entire process. The shaded boxes represent procedures that occur on digital devices, while the white boxes represent procedures that occur on analog devices.
Figure 15 illustrates the coupling between HMC unit cells through mutual inductance.
Figure 16 illustrates direct inductive coupling between HMC unit cells.
Figure 17a illustrates a schematic of a possible physical realization of a single unit cell consisting of noisy resistors and capacitors, illustrating direct inductive coupling between HMC unit cells.
Figure 17b illustrates a schematic of a possible means of coupling two unit cells using coupling resistors.
Figure 18 illustrates an example of a switched capacitor bank controlled by digital logic, where each physical capacitor adds an additional bit of precision to the overall capacitance.
Figure 19 illustrates two varicaps in a back-to-back configuration with the third port in the middle acting as a tuning voltage to adjust the capacitance of the system.
Figure 20 is a circuit diagram of a digital-analog method for implementing bipolar coupling.
Figure 21 shows an active inductor based on a gyrator (top) and an equivalent circuit with an inductor (bottom).
Figure 22 illustrates a schematic of a basic symmetric RC unit cell used for over damped Langevin sampling of a one-dimensional normal (i.e., Gaussian) probability distribution.
Figure 23 shows the schematic of a basic coupling cell used to achieve negative coupling between two RC unit cells without an inductor.
Figure 24 shows a schematic of a basic coupling cell used to achieve positive coupling between two RC unit cells without an inductor.
Figure 25 illustrates different strategies for achieving local connectivity, including a square lattice, a hexagonal lattice, and a King's graph, where each box represents a unit cell.
Figure 26 illustrates additional strategies for achieving local connectivity based on a cluster graph containing square clusters and hexagonal clusters, where each box represents a unit cell.
Figure 27 illustrates a semi-local approach to connecting unit cells (represented by five boxes). Here, unit cells within a row are fully connected, while connections between rows are more sparse. For D = 9, the number of nearest neighbors is n = 3 or n = 4 along the row. For D = 16, it is n = 4 or n = 5 along the row. In general, using this connection scheme, n It is scaled accordingly.
Figure 28 illustrates a circuit schematic of a fully connected architecture.
Figure 29 is a graph of the success rate of the accept/reject step as a function of the number of precision bits in the switched capacitor bank capacitors (error bars represent the standard deviation error). The results are obtained by simulating the process and hardware using an 8-dimensional fully connected capacitance matrix.
Figure 30 is a graph of the success rate of the acceptance/rejection stage as the capacitor defects increase (the shaded area represents the standard deviation error). It is the result of the process and hardware simulation using a 200-dimensional fully connected capacitance matrix. Capacitors have various levels of tolerance, but one can choose a capacitor with, for example, a 2.5% tolerance, which according to the graph will result in a success rate of about 80%.
Figure 31 shows plots (error bars represent standard deviation) illustrating the impact of hardware with limited k -banded connectivity on the acceptance/rejection phase acceptance rate. The results are simulated using process and hardware with increasing dimensionality of k -banded capacitance matrices. (A) is the case where the target distribution has a sparse covariance matrix with elements decreasing with distance from the diagonal; (B) is the case where the target distribution covariance matrix is arbitrarily dense.
Fig. 32 shows the dispersion The target distribution f _t is two distributions f ₁ with variances ε and 2ε, respectively. and f ₂ It is a curve approximated by interpolation between f _t , f ₁ , and f ₂ are Gaussian, and f _a is a Gaussian mixture.
Figure 33 is a plot of the distance L ₁ between the approximate probability density function f _a and the target f _t as a function of ε . The target distribution is a univariate normal distribution with mean 0 and variance 20. With error relaxation (green), the first derivative of the error tends to 0 as ε tends to 0, but without error relaxation (blue), the first derivative does not tend to zero as ε tends to 0.
Figure 34 illustrates a Maxwell's Demon, a concept in thermodynamics whereby an intelligent observer can observe a system and reduce the entropy of the system by applying appropriate forces to the system. A typical example is the separation of a gas mixture, as shown here.
Figure 35 illustrates how the potential energy landscape U(x) can be changed by adding a Maxwell demon. The lower panels show schematically how this changes the probability distribution P(x) from which samples are drawn. The left panel shows the harmonic oscillator potential corresponding to sampling a Gaussian probability distribution. The right panel shows that adding a Maxwell demon system produces more complex potentials, allowing for sampling other distributions, such as multimodal distributions.
Figure 36 illustrates how a Maxwell demon device can be viewed as a voltage-controlled voltage source.
Figure 37 illustrates how a digital Maxwell Daemon system (e.g., processed on a CPU or Field Programmable Gate Arrayt (FPGA)) interacts with analog HMC hardware.
Figure 38 illustrates an uncorrelated approach to constructing an analog Maxwell's demon. A primary unit cell (left) is coupled to an auxiliary unit cell (right). Specifically, the voltage across the capacitor of the primary unit cell is applied as a voltage source to the auxiliary unit cell. This voltage source drives a current through a non-linear element (NLE), and the current is read as a voltage across a resistor. The latter voltage difference, after passing through an inverter, acts as a voltage source in the primary unit cell, corresponding to the output of the Maxwell's demon.
Figure 39 illustrates an alternative configuration of an analog Maxwell Demon without a voltage inverter. That is, by making the resistors of the auxiliary unit cell ungrounded at both terminals, the voltage difference can be extracted and used as a voltage source for the primary unit cell after inverting the wires.
Figure 40 illustrates the characteristics of a diode-based Maxwell's demon. Here, the NLE is assumed to be a forward-biased diode. The first panel shows the current-voltage characteristic of the diode, denoted I(x) . The second panel shows the force f(x) output by the Maxwell's demon, which is assumed to have a negative sign on the I(x) curve. The third and fourth panels show the potential energy U(x) and the sampled probability distribution P(x) , respectively.
Figure 41 illustrates the characteristics of Maxwell's demon based on a tunnel diode. Here, we assume that the NLE is a forward-biased tunnel diode. The first panel shows the current-voltage characteristic of the tunnel diode, denoted I(x) . The second panel shows the force f(x) output by Maxwell's demon, which is assumed to negatively affect the I(x) curve (as discussed in the text) and then shift the voltage by a constant negative amount. The third and fourth panels are the potential energy U(x) and the sampled probability distribution P(x) , respectively.
Fig. 42 illustrates an analog Maxwell daemon for generating correlations between unit cells. An analog coupling operation, which may be an affine transform for example, is applied to the voltage vector associated with the capacitors of each unit cell. After the analog coupling operation, the resulting vector is applied component-wise as a voltage source to each auxiliary unit cell containing a nonlinear element (NLE). The current passing through each NLE is read as a voltage, which is then inverted and finally applied as a voltage source to the corresponding primary unit cell.
Fig. 43 illustrates an analog circuit for implementing matrix-vector multiplication A v. This provides a possible implementation for the analog combine operation illustrated in Fig. 42. There are d resistors associated with each row of matrix A. This gives a total of d² resistors for the entire matrix A , assuming that A is dense.
Figure 44 illustrates an analog circuit that implements a function h ( v ) using transistors. Here, the ith component of the circuit that implements h _i ( v ) is shown as a representative example of other components. Here, there are d transistors associated with h _i , so that the overall conversion This gives a total of d² transistors for h ( v ). Each component v _j of v acts as a gate voltage for the transistor associated with the other components v _k , resulting in higher-order coupling. (For simplicity, only three components of v are shown, but there are d components.)
Figure 45 illustrates a device for performing analog HMC when the probability distribution that the user wishes to sample is represented by a neural network.
Figure 46 illustrates a device for performing analog HMC on a Gaussian distribution, where the user specifies a precision matrix Q = Σ ^-1 instead of a covariance matrix Σ.
Figure 47 illustrates three unit cells coupled through three ports of a voltage controlled capacitor.
Figure 48 illustrates four unit cells coupled through a voltage multiplier (i.e., voltage mixer) and a voltage control capacitor. The voltage multiplier can be implemented on a digital device or as an analog device.
Figure 49 illustrates a method for achieving local coupling through voltage-controlled capacitors and three-body coupling via a square grid.
Figure 50 illustrates a method for achieving local connectivity using voltage controlled capacitors (VCC) and three-body coupling via a hexagonal grid.
Figure 51 illustrates an LMC hardware unit cell with a capacitor, an ideal resistor, and a stochastic voltage source.
Figure 52 illustrates two unit cells of LMC hardware coupled with either (A) a resistive bridge or (B) a capacitive bridge.
Figure 53 illustrates two coupled LMC hardware unit cells for Gaussian sampling. The covariance matrix of the Gaussian is encoded into the values of the circuit elements, in particular the resistances of the resistors.
Figure 54 illustrates an alternative analog device based on integrators for implementing LMC. The user interfaces with this device through compiling the slope of the log probability and collecting samples.
Figure 55 illustrates a possible circuit diagram for two coupled unit cells forming the building blocks of the hardware. In each unit cell, current noise is injected in parallel with the LC oscillator circuit.
Figure 56 is a schematic diagram of the thermodynamic process of solving a system of linear equations.

1 System outline

FIG. 2 illustrates a thermodynamic computing system (200), also known as a thermodynamic computing platform, a thermodynamic system, or a thermodynamic platform, which can sample a multivariate Gaussian distribution using the Hamiltonian Monte Carlo (HMC) method, which is a state-of-the-art Markov chain Monte Carlo (MCMC) method. The HMC sampling method is a Markov chain Monte Carlo (MCMC) method that samples from a Gibbs distribution. The HMC sampling method reduces the problem of sampling from a Gibbs distribution of f to simulating a particle trajectory according to the Hamiltonian dynamical laws, where f acts as a potential energy function. The samples generated by the thermodynamic computing system can be used in probabilistic AI, for example, to weight connections between neurons in various layers of a probabilistic neural network.

The thermodynamic computing system (200) of FIG. 2 is a hybrid digital-analog system comprising a digital system (210) interfacing with an analog system (220), also referred to as an analog dynamical system or a time-evolving analog system. The digital system (210) may be implemented in suitably programmed digital hardware (e.g., a processor, an application-specific integrated circuit, a field-programmable gate array, etc.) and executes both a probabilistic model (212), which is an application-level statistical model inferred by the thermodynamical system, and a controller (214) interfacing with the analog system. The controller may be implemented in firmware. In operation, the controller (214) initializes control parameters (211) of the analog dynamical system (220) and accepts or rejects samples (221) of a target multivariate distribution from the analog system (220). Additionally, the digital system (210) interfaces with the user directly or through another digital computing device (e.g., a personal computer (PC)).

The analog system (220) comprises a dynamical physical system (222) that evolves through Hamilton's equations over a series of time steps to provide samples of a target multivariate distribution, such as a multivariate Gaussian distribution, a Bernoulli distribution, a binomial distribution, a Poisson distribution, or other exponential distribution. More generally, the analog system (220) can sample any distribution that can be defined in terms of an energy function that can be realized in (analog) electronics. The target multivariate distribution can be a high-dimensional distribution (e.g., a distribution having dimensions of 10, 100, 1000, 10000, or more) representing a probability, heat, a population, or other suitable quantity.

At each time step of a series of time steps, the analog system (220) proposes a sample (221) of the target probability distribution to the controller (214), and the controller either accepts the proposed sample or rejects it if it is an outlier. The natural temporal evolution of the analog system (220) maps to the Hamiltonian dynamics, which means that the analog system naturally implements the integration step of the HMC sampling method. In other words, the natural physical dynamics of the analog system replaces the numerical integration traditionally performed in digital hardware. As a result, the joint digital-analog thermodynamic system solves these computations via the temporal evolution.

A controller (214) of the digital system (210) sets input parameters (211) of the analog system, which parameters represent a multivariate distribution, for example a normal (Gaussian) distribution or another distribution from which the analog dynamical system (220) samples. The input parameters (211) may include inductance, capacitance, voltage source, and current source of adjustable inductors, capacitors, voltage sources, and current sources, respectively, that make up the analog system (220), as described below. The controller (214) sets these input parameters (211) via a serial bus or other suitable interface, and then allows the analog system (220) to evolve for a predetermined period of time, also known as a time interval, integration period, or sampling period. Once that period of time has elapsed, the controller (214) reads out output parameters of the analog system via the serial bus, which may include voltages, currents, and/or other analog values representing the final state of the analog system. These output parameters represent samples (221) of a multivariate Gaussian probability distribution. The controller (214) rejects incorrect or outlying parameter readings to correct for analog imperfections of the analog system (220).

The analog component (220) of the thermodynamic system (200) simulates physical (Hamiltonian) dynamics to sample a multivariate Gaussian distribution or other target probability distribution. The Hamiltonian (H) of the physical system (222) specifies the total energy of the physical system, which includes both kinetic and potential energy. In the context of sampling, the Hamiltonian represents a target probability distribution with q sample space variables in D dimension. The HMC method traditionally derives the derivatives of the Hamiltonian with respect to the sample space variables and generalized momentums ( p ), and then derives trajectories for time t over a duration ( ) involves integrating over a series of time steps. In this hybrid digital-analog approach, the dynamics of the physical system replaces both the differentiation and the trajectory integration.

One such physical system (222) is an analog electronic system, particularly a network of coupled harmonic oscillators. As described in more detail below, the kinetic energy of the network of coupled harmonic oscillators is mapped to an energy function (Hamiltonian) of a multivariate Gaussian probability distribution up to a normalization constant. The voltages and currents of the oscillators represent position coordinates of the Hamiltonian of the multivariate Gaussian probability distribution, and so they can be sampled to provide samples of that distribution.

2 1st dimension Normal probability distribution sampling

A single lumped circuit element, such as a basic LC oscillator, can be mapped to a one-dimensional (1D) normal (i.e., Gaussian) probability distribution as follows. As a result, the LC oscillator can be used as an analog dynamical system for sampling a 1D normal probability distribution in a thermodynamic computing platform that performs HMC calculations.

2.1 For the 1D normal probability distribution Hamiltonian Monte Carlo (HMC)

Considering a normal distribution with mean μ and variance σ ² , the probability distribution can be written as

The energy function corresponding to this system is:

Let p be the canonical momentum corresponding to x , and choose a kinetic energy function independent of x . Standard choices are as follows:

Here, m represents the mass of the virtual particle with position x . Then, the total energy or Hamiltonian of the system is

After specifying the Hamiltonian, HMC relies on two principal components to sample from the target distribution:

Hamiltonian dynamics

Given the Hamiltonian ( H ) and initial conditions ₍ x0 and p0 ) for x and p , HMC integrates the Hamiltonian equations over time _as follows:

At the end of the integration, the new values of x and p are the next proposals in the MCMC chain.

Normal distribution

Independent of the dynamical problem, given the Hamiltonian ( H ) and the temperature ( T ), there is a corresponding canonical probability distribution for x and p ,

Here, k _B is the Boltzmann constant.

The structure of HMC is influenced by the property that x and p are independent random variables. This property follows from the fact that U ( x ) and K ( p ) do not overlap in the support interval, as follows:

Connecting the two components is the choice of initial conditions for the integration. The initial values ( x ₀ and p ₀ ) set up a surface with a constant probability density over which the Hamiltonian dynamics moves. Up to normalization, the probability density along the trajectory is

The independence of the marginals allows us to choose a new constant value for the joint probability density by resampling p ; then, Hamiltonian dynamics exchanges this momentum resampling for a position resampling. To quote Neal et al ., "MCMC using Hamiltonian dynamics," in: Handbook of Markov Chain Monte Carlo , 2.11 (2011), "[s]ince [ x ] isn't changed, and p is drawn from its correct conditional distribution given [ x ] (the same as its marginal distribution, due to independence), this step obviously leaves the canonical joint distribution invariant."

2.2 Hamiltonian and electrical circuit

The following subsections describe the three steps for deriving the Hamiltonian from a lumped element circuit.

2.2.1 Selecting Coordinates

Assume that we start with an LC oscillator, as in Fig. 3. The LC oscillator has an inductor with inductance ( L ) and a capacitor with capacitance ( C ). We would like to write down the Lagrangian for this circuit, which will enable us to derive the corresponding Hamiltonian.

The method begins by selecting coordinates for all branches and nodes in the circuit. The branch coordinates ( I _L and I _C ) are chosen to be the currents in the inductor and capacitor, respectively. One node is designated as the reference (ground) and the other nodes are designated as the coordinate V , i.e., the voltage difference from that node to ground.

2.2.2 Lagrangian Induction

Next, we can write the total energy of the circuit using these selected coordinates. For a circuit, the total energy is the sum of the energies of the branches. From Kirchhoff's laws, I _{L =} - I _C is established. According to the definition of a capacitor, I _C = C ( dV /dt ). Therefore, the energy of the capacitor and inductor can be written in terms of the voltage V and its time derivative.

Since voltage is treated as a position coordinate, the capacitor energy is potential energy ( U ) and the inductor energy is kinetic energy ( K ).

Now we can write the Lagrangian for our system as:

Given this Lagrangian, the momentum conjugate p with respect to position coordinate V is as follows:

Therefore, the Lagrangian can be rewritten in terms of position and momentum as follows:

2.2.3 Hamiltonian To get Legendre Conversion

Since we are trying to do a Hamiltonian Monte Carlo, the Hamiltonian of the physical system must match the Hamiltonian corresponding to the probability distribution. To do this, we can now derive the Hamiltonian of the LC oscillator.

The Hamiltonian function of a 1D system is defined as follows:

From this it is derived that

The momentum conjugate p can be written in terms of the current through the capacitor as p = LCI _C.

2.3 Gaussian LC matrix On the oscillator connection

Recall the two Hamiltonians derived earlier. One is the following abstract Hamiltonian corresponding to a Gaussian with virtual mass m ,

The other is the Hamiltonian of the electrical system.

Given these two Hamiltonians, we want to find the values of L and C of the physical Hamiltonian using the variance and mass of the Gaussian Hamiltonian.

First, we note that the capacitance C is equal to the inverse variance. The mean can be implemented in hardware using a constant voltage offset, or in software by simply adding the mean to the other zero-mean samples. The mass is equal to LC ² . Thus, after choosing the mean and variance, the inductance L is a free parameter that allows us to choose the mass.

This is an example of how to map a probability distribution to an equivalent electrical circuit.

3 capacity Combined Harmony With the oscillator Multivariate normal probability distribution

3.1 Distribution Definition

An N-dimensional multivariate normal distribution is an N-dimensional mean vector It is defined using the N×N covariance matrix ∑ .

now = 0. Through this, the multivariate normal distribution with covariance matrix ∑ can be written as follows.

Taking the log gives the following "energy function" corresponding to this distribution:

Here Z = is a constant. As with the one-dimensional Gaussian, an electric circuit can be implemented that fits this energy function.

3.2 Capacity Combined Oscillator

3.2. 1 circuit Lagrangian

Consider a set of N harmonic oscillators, each capacitively coupled to all other oscillators. Any pair of harmonic oscillators i , j can be represented by the subcircuits illustrated in Fig. 4.

The total inductive energy of this circuit is:

If we define the induced matrix as follows:

The total induced energy can be written as

The total capacity energy is:

To simplify the subsequent calculus, the Maxwell capacitance matrix C is defined as follows.

Then, we can rewrite the capacity energy as follows:

We assume that the current through the inductor is treated as a position coordinate. To describe the Lagrangian, we rewrite the voltages using the time derivatives of these position coordinates. For each voltage, There is a relationship of . Therefore, the capacitive energy can be written as follows.

Now the Lagrangian of the system can be written as follows:

3.2. 2 circuits Hamiltonian

The Hamiltonian is a set of N positions q _i , N momenta p _i , and a Lagrangian It is defined as follows using

Here, the momentums are defined as follows:

In this case, the momentum is:

Here we use the following line by line: the definition of conjugate momentum; the product rule for derivatives; the reduced and relabeled indices of the Kronecker delta; and the fact that C is a symmetric matrix. This can be summarized in the matrix equation as follows:

This means that

Then the Lagrangian can be rewritten as follows:

This means that the Hamiltonian of the system is:

In terms of circuit parameters, the voltage across the inductor is:

So the momentum is

And the voltage using momentum is as follows.

Then the Hamiltonian can be computed as follows:

3.3 Change coordinates

Existing code bases for HMC and its extensions assume that the target distribution is encoded by potential energy. Therefore, if the distribution can be encoded similarly to potential energy, the design implementation effort can be greatly reduced. This involves changing the coordinates used to describe the analog electrical system. A new coordinate system The original coordinates are such that the resulting stationary paths of the new Hamiltonian K are the same as the stationary paths of the original Hamiltonian. And what changes can be made to the Hamiltonian H ?

3. 3. 1 General Theory

To begin, consider whether the original Hamiltonian system satisfies:

Here, the integrand is the Lagrangian expressed using the Hamiltonian. Considering an arbitrary transformation of the system, the old and new Hamiltonians are related by:

This is Scala is a function of phase space coordinates F with continuous second derivatives. We now show explicitly that the new Hamiltonian K has the same stationary path as H.

Since it is also Hamiltonian, the following holds.

Here we apply the following line by line: the definition of Hamiltonian dynamics for K and substitution using H ; linearity of integration and differentiation; fundamental theorem of line integrals; independence of differences in endpoint values as a function of changes in the integration path. The static paths of both sides of this equation are identical: this means that the location of the static path is a multiplicative constant. Because it does not depend on .

3.3.2 New Coordinate Proof

Based on the results of the previous section, we had a procedure to check whether the transformation was canonical. Coordinates Given a new Hamiltonian K with a variable , it is a scalar By proving a function G that satisfies equation (5), we can obtain another Hamiltonian It can be confirmed that the same dynamics are induced.

For the given system, we consider the original Hamiltonian as follows

If you specify new coordinates as follows,

This gives the following transformed Hamiltonian:

Because we assume L is symmetric.

and Let it be. Then, we get the following

Using this definition, we check whether equation (5) is satisfied as follows:

Here we apply the following line by line: definition of K ; definition of coordinate transformation and definition of F ; elimination of opposite terms; and definition of H and Definition of . Therefore, this transformation is valid.

3.3.3 Final physical and abstract relationships

The definition is summarized as follows. There is an abstract Hamiltonian such that

Here, the coordinates are related to the circuit quantities as follows:

and

3.3.4 Find the inductor current based on the capacitor current.

Consider again the circuit of Fig. 4. The current at node i is:

The current passing through the capacitor is and so it becomes the following

The sign matrix is defined as follows.

Therefore, we can combine the sums by writing as follows:

Converting to matrix form gives:

3.4 Connection with multivariate normal distribution

If the values of the circuit are set to C = Σ, the potential energy of the electric system in the transformed coordinate system is equal to the energy of the target multivariate normal distribution.

Maxwell's capacitance matrix is defined as:

This means that the following relationship holds equally when discrete capacitors are variables:

Applying the relation in reverse, when the covariance matrix Σ is given, the individual capacitances are set by first setting the coupling capacitors as follows.

Then the capacitances for grounding can be set as follows:

3. 5 regular mode

Consider an example of an HMC process running on a 2D Gaussian, shown in Fig. 5. This figure illustrates the importance of choosing the integration time carefully: integrating for too long will result in returning to a location that has already been visited, which means wasted computation.

Typically, this is addressed by upgrading the HMC to a No U-Turn Sampler (NUTS), which adaptively sets the integration time to avoid these U-turns. However, since the Gaussian distribution is so simple, it is conceivable to compute a reasonable integration time directly.

For this purpose, we assume a coupled oscillator system with mass matrix M and spring constant matrix K. In terms of the physical mass-spring system, we obtain:

and

Then the oscillator position vector The equation of motion can be written as follows with as a variable.

Solving the equations of motion yields the unknown eigenvector The generalized eigenvalue equations for the normal angular frequency ω are derived with as a variable:

This can be reorganized as follows:

A = K ^-1 M and If we define , the standard eigenvalue equation becomes:

Each eigenvalue About By setting the regular frequencies, the normal frequencies are obtained.

Given these normal frequencies, we can know how long we need to integrate before we start a U-turn in the direction of slowest motion (lowest frequency). Let f _min be the minimum normal mode frequency. Then, after one cycle, we are back where we started.

is the upper limit of the integration time. So, a smaller value is chosen as the integration time.

for example,

or one-quarter of the longest period.

3.6 Frequency Scaling

To control the frequency of the oscillator, both the capacitance and inductance can be adjusted to reduce the frequency of the oscillator. Because it is determined by . The proof of the Cholesky decomposition process can be modified to extract samples from a scaled version of the desired Gaussian, which gives a scale parameter η that can adjust the oscillator frequency.

3.6.1 General Theory

average and using the covariance matrix ∑ In , we define the target Gaussian probability density g ,

The bijective continuous differentiable transformation φ is defined as follows.

Suppose we define a new function f as follows:

If X is some random variable with density f , then φ ( X ) has density g . So let's see if it is convenient to sample from the density f .

The following holds true:

,

That is, η times the identity matrix is established. Therefore, the following is established.

For composite functions, the following holds:

. then

and If we multiply by , we get:

This is the covariance matrix is the probability density of the centered Gaussian.

3.6. 2 circuits Parameter perspective

Thus, the scaling procedure is established. Consider that the frequency increase of a physical oscillator occurs by decreasing one or both of its inductance or capacitance. Since we have set Σ _physical = C for capacitively coupled oscillators, for a fixed inductance L, Σ _physical < Σ Setting it will result in a higher frequency.

In other words, if η > 1 and the following holds, the frequency becomes higher

According to the above derivation, the samples extracted from the system are will be distributed according to (0, ∑ _physical ); each sample cast If you map it to Samples distributed according to are obtained.

3.6.3 Quantifying Frequency Scaling

How much has the frequency changed? Recall that the eigenvalue equations of the coupled oscillator system are

In terms of the system covariance matrix, this becomes

Substituting the target covariance Σ and the scaling factor η , we get

If ω represents the normal mode frequency of the original target distribution,

then

Therefore, the coefficient η used in the definition of φ is exactly equal to the product of the normal mode frequencies.

4 Physical Implementation example

FIG. 6 illustrates a hardware implementation that primarily performs the sampling portion of the overall computational process. The hardware implementation includes a digital system, implemented herein as an FPGA controller (610), operably coupled to an analog system, illustrated herein as an LC resonator (620), coupled to an analog-to-digital converter (ADC), a digital-to-analog converter (DAC), a noise source, and a switch (630). The samples are voltage levels across the LC resonators (620) (also referred to as unit cells or cells, as described in detail in Section 4.1) when excited by the noise source (630). A user may interact with the hardware implementation via a personal computer (PC) (602) operably coupled to the FPGA controller (610).

The PC (602) downloads three different types of information from the FPGA controller (610): (1) cell resonance frequency, (2) coupling sections and their respective polarities, and (3) noise level. The FPGA controller (610) (described in detail in Section 4.3) sets up the circuits according to the downloaded data. The ADC samples the voltages and stores the data in the on-chip memory of the FPGA. Once the FPGA controller (610) has collected a predetermined amount of data, it transmits the data back to the PC (602).

The interface between the FPGA controller (610) and the PC (602) illustrated in Fig. 6 is a virtual COM port via USB. Other interface protocols (e.g., PCIe or CXL) are also possible.

4.1 Combined Resonators

Figure 7 illustrates coupled resonators in the form of two capacitively coupled parallel LC circuits. The coupling and cell capacitances are determined by the controller. The capacitors can be adjusted analog or digitally, while the inductances can be kept constant.

Coupling a larger number of resonators improves system performance. The nature of the coupling (e.g., all-to-all, linear, etc.) depends on desired processing and hardware constraints. Figure 8 illustrates a conceptual all-to-all coupling between N cells.

The exact implementation examples of cells and their coupling sections will vary depending on the technology. Sections 4.1.1 and 4.1.2 show examples of how to implement cells with technologies compatible with printed circuit boards (PCBs).

4.1.1 Tunable Resonator

Fig. 9a(a) illustrates the structure of a switch-tunable resonator. The capacitors can change their values using shunt switches (e.g., transistors). The minimum capacitance value ( C 1) does not require a switch. The top plates of the capacitors represent the voltages that are sampled and coupled to other cells.

The noise source can be implemented using a digital controller to reduce or minimize the dependence on analog circuits. As a result, a lowpass filter (LPF) filters out unwanted and/or undesirable high frequency components of the output signal. This is described in detail in Section 4.4.

4.1.2 Control binding

The coupling units between the cells can be implemented as illustrated in Fig. 11. A series switch (controlled by a coupling enable signal) can turn the coupling on or off. The polarity of the coupling is controlled by the signal coupling polarity . Negative and positive polarity are achieved using a transformer setup as illustrated in the drawing.

4.2 Analog-to-digital converter ( ADC )

Voltage readings from each cell are made using an ADC. This configuration uses a single ADC channel for each cell. Alternatively, multiple cells can share a single ADC if they are sampled simultaneously. The input impedance of the ADC must be considered relative to the resonant frequency of the cells.

4.3 Controller (FPGA)

The controller within the PCB handles the following functions:

Download and upload interface with PC

Interfacing with ADC (reading and storing data, initialization, etc.)

Control cell switches to tune the resonators

Control to turn the coupling switches on or off with the appropriate polarity

Generating pseudorandom noise (Section 4.4)

Figure 12 shows a block diagram of the FPGA functions along with the protocol of setup data from a PC.

4.4 Noise injection

Each unit cell has an independent noise source. Since the cells are set up as a parallel LC circuit, the noise sources in Fig. 7 are current sources. One way to implement this is to generate a single bit pseudorandom noise (PN) and filter the output using an RC circuit, as shown in Fig. 13a. The presence of a resistor also causes the PRN voltage source to appear as a current source (Thevenin to Norton transformation).

The noise amplitude can be controlled using pulse density modulation (PDM) for a single-bit noise source.

4.4.1 PN Polynomial

A single-bit pseudorandom sequence is implemented as a linear feedback shift register (LFSR) configuration representing a primitive polynomial. Each cell starts with a different initial condition, causing the cells to generate uncorrelated pseudonoise sequences. The length of the shift register must be large enough so that the noise of each unit cell is sufficiently separated from other sequences. A typical shift LFSR is illustrated in Fig. 13b.

5 How it works

5.0. 1 System reset

At this stage, the device is prepared to extract samples from a specific multivariate normal distribution according to the process illustrated in Fig. 14, starting from system initialization (1410):

1. Definition of the probability model (1411) (digital device): Target covariance matrix Σ and target mean vector is a d- dimensional normal distribution Defines .

2. Hyperparameter selection (1412) (digital device): Σ using the Lanczos process, and Calculate the highest and lowest normal mode frequencies of the oscillator block. The frequency is the highest possible frequency, taking into account the characteristics of the adjustable components and the measurement circuits of the oscillator block implementation. _Threshold and to make them almost the same, that is Choose a scaling factor η such that it becomes _{a threshold} .

Why η ? This allows the speed of the device to be chosen independently of the target covariance matrix. Without η , the frequencies are fixed by the normal mode equation (15).

_{The threshold} is set by the fact that there is some trade-off between device speed (set by η ) and device accuracy (set by the measurement circuitry); the exact value is derived by simulation of the entire oscillator system.

Set the target physical covariance matrix as Σ _physical = η ^{- 2} Σ (see Section 3.5 for details). Integration time to be used for the modified HMC sampling process. (For more information on this selection, see Section 3.6).

3. Oscillator Set component values (1413) (Analog devices) : Select the setting L _i of each adjustable inductor F (e.g., all ) These values form a diagonal matrix L. Compute the capacitance matrix, which is equal to C = ∑ _physical . Each coupled capacitor Set C _ij = - C _ij value. Each i -th adjustable oscillator capacitor E Set to .

4. Markov Chain Initialization (1414) (Digital Device) : Initialize the Markov chain S on a digital device. The initial state is a vector of all zeros whose length is d , which is the dimension of the target normal distribution. This chain represents voltage values.

5.0.2 Sample Extraction

In sampling step 1420, the digital device (the digital controller) and the physical device (the coupled network of analog unit cells) cooperate to perform a Markov chain for the voltages. The iterative part of the procedure is an application of the Hamiltonian Monte Carlo (HMC) process as follows: instead of positions, momentum is used as the chain state, and the integration is performed by the physical device through its natural dynamics rather than by the digital device through numerical integration. The accept/reject step of the loop is a novel form of analog error correction because it corrects for deviations of the physical device from the target distribution. The following steps are repeated N times:

1. V selection (1421) (digital device) : The most recent state in the Markov chain S is the next voltage value. (initial voltage to set the coupling spike J between each oscillator).

2. I Selection (1422) (digital device) : d uncorrelated unit normal distribution Extract samples from , resulting in a vector of length d are generated. Currents The initial value of It is calculated as .

3. Initial V, I setting (1423) (Analog device): Open switch B and close switches A and C. Voltage sources D , and the current sources G are set to Set to .

4. Activate System Evolution (1424) (Analog Device) : Turn on the switch so that B is closed and A , C are open in each oscillator block; turn off the voltage and current sources. This allows the oscillators to start evolving freely.

5. Final V, I measurement (1425) (Analog device) After time T , measure the final trajectory values - the final voltages of the coupled spikes J using each voltmeter X. Measure _the final currents through the oscillator inductors F using each ammeter Y. Measure _{the final} .

6. Calculation of measured value (1426) (digital device) : Initial value of energy _Initial = Wow, the final value of energy Final= Calculate the acceptance rate Computes . Sample a value u from a uniform distribution on the interval [0, 1]. If u < α , then -1. 0* Add _{the final} to the Markov chain S ; otherwise, Add .

5.0.3 End Sequence

In the termination sequence step (1430), the physical quantities are converted back into target distribution samples.

1. Termination (1431) (digital device) : Once N samples are extracted, stop adding values to the Markov chain.

2. Momentum calculation (1432) (digital device) : Change variables from electrical quantities (voltages) to algorithmic quantities (momentum). Each of the Markov chain S About conversion This is done by applying .

3. Markov chain regularization (1433) (digital device) : Angular momentum Convert to and returns the result.

6 Inductive coupling

6.1 Inductive for multivariate distribution sampling Combined Harmony Oscillators

6.1. 1 circuit Lagrangian

Consider a set of N harmonic oscillators, each inductively coupled to all other oscillators through mutual inductances L _ij . Any pair of oscillators i , j can be represented by the following subcircuit:

The total induced energy in this circuit is:

The total capacitive energy is:

To simplify the subsequent calculus, the inductance matrix L is defined as follows.

Then, we can rewrite the induced energy as follows:

We assume that the voltages between the nodes of the circuit are treated as position coordinates. To describe the Lagrangian, we rewrite the currents using the time derivatives of these position coordinates. For each current, There is a relationship called .

Therefore, the induced energy can be written as follows:

Now the Lagrangian of the system can be written as

6.1. 2 circuits Hamiltonian

The Hamiltonian is a set of N positions q _i , N momenta p _i , and a Lagrangian is defined as follows using

Here, the momentums are defined as follows:

In this case, the momentum is:

Here we use the following line by line: the definition of conjugate momentum; the product rule for derivatives; the reduced and relabeled indices of the Kronecker delta; and the fact that L is a symmetric matrix. This can be summarized in the matrix equation as follows:

This means that

Then the Lagrangian can be rewritten as follows:

This means that the Hamiltonian of the system is:

6.1.3 Direct inductive coupling

Consider a set of N harmonic oscillators, each of which is directly inductively coupled to all other harmonic oscillators. Any pair of oscillators i, j can be represented by the subcircuits shown in Fig. 16.

The total induced energy in this circuit is:

The current in the coupling branch here is:

This means that the induced energy does not couple neighboring cells.

6.2 Inductive coupling and Capacitance combination of combinations

More generally, a multivariate normal probability distribution is mapped onto a network of coupled harmonic oscillators, whether the coupling is inductive, capacitive, or both. Again, the network contains one harmonic oscillator for each dimension of the multivariate normal probability distribution, and the mapping is based on the Lagrangian and Hamiltonian of the circuit.

For this general case of inductive coupling and/or capacitive coupling, the derivation of the Lagrangian and Hamiltonian is conceptually similar to the derivation presented above. Similarly, the resulting Hamiltonian can be mapped to a multivariate normal probability distribution following the derivation presented above. Therefore, for simplicity, we omit the derivation and note that it is possible to combine inductive and capacitive coupling.

7 Resistance Combination

An alternative way to couple the unit cells is to use resistors. In this case, one can consider a unit cell with a noisy resistor at non-zero temperature. Figure 17(a) shows a typical equivalent noise model for a noisy resistor consisting of a stochastic voltage noise source δv ( t ) in series with an ideal (noiseless) resistor of resistance R. The terminal capacitance C of the resistor is also added to the equivalent resistor model. The voltage at node 1 (labeled here simply as v ( t )) is a state variable, and its dynamics are as follows:

This stochastic differential equation (SDE) contains a drift term proportional to v ( t ) and a diffusion or stochastic term proportional to δv ( t ).

When constructing a system consisting of many unit cells, the unit cells must be combined to express the correlation and geometric constraints. For example, as shown in Fig. 17(b), two unit cells can be combined through a resistor. Then, the combined unit cells represented by the voltages at nodes 1 and 2 are combined through the following drift terms:

Here and

Here, the self-resistance matrix R , the capacitance matrix C , and the conductance matrix J are introduced. Then, a system of a plurality of coupled unit cells can be constructed using the basic building blocks illustrated in Fig. 17.

8. Adjustable configuration of components

In order to be able to represent any Gaussian distribution, hardware components must be changed or adjusted to suit each distribution. The selection of components must be made with a balance between precision and range in mind.

8.1 Adjustable configuration of capacitors

There are two ways in which tunable capacitors are implemented in HMC hardware: (1) using a switched capacitor bank, and (2) using variable capacitor (varactor) diodes.

8.1.1 Switching capacitor bank

Fig. 18 illustrates a switched capacitor bank. The switched capacitor bank consists of a number of capacitors that can be connected in parallel in any combination to obtain a number of different capacitance values. A digital signal is used to actuate the switches to connect or disconnect the capacitors. Thus, the device essentially operates as an adjustable capacitor whose capacitance value is adjusted by a digital voltage signal, as illustrated in Fig. 18.

Since the device can contain any number of capacitors of any capacitance value, it can theoretically have any range and precision at the expense of large physical area, losses, and parasitic capacitors. Furthermore, since the device is constructed from ordinary capacitors, it generally has good linearity and operates well over a wide range of voltage amplitudes.

8.1.2 Varactor Diode

Varicap diodes are semiconductor devices with voltage-dependent capacitance values. This type of diode has a depletion layer that is particularly sensitive to voltage bias. Therefore, the capacitance value can be changed by varying the voltage bias across the diode. These diodes are specially designed to meet the desired tunability range. Varicaps generally have a narrower tunability range than switching capacitor banks, but they have continuous tunability.

Another consideration of these semiconductor devices is their nonlinearity. This means that they behave like linear capacitors for very small voltage fluctuations. For larger AC fields, the nonlinear distortion becomes severe. Finally, since these devices are diodes, they are polarized. That is, they behave as capacitors only in one direction of the DC voltage bias.

Figure 19 illustrates a possible approach to constructing a tunable capacitor based on varicaps. Here, the two varicaps are oriented in opposite directions, and an external source applies a voltage between the two varicaps. This approach allows controlling the overall capacitance using an external voltage, such as a voltage signal from a digital device.

8.2 Bipolar Configuration of coupling capacitors

Any objective covariance has some elements that are negative, indicating negative correlations between variables. Bipolar (positive and negative) capacitive coupling allows the representation of these elements in the hardware capacitance matrix.

8.2.1 Transformer Method

As discussed above with respect to Fig. 11, negative coupling between elements can be implemented using an adjustable capacitor connected to a transformer. The adjustable capacitor sets the magnitude of the coupling while the transformer changes the sign. The transformer is wound so as to have the effect of nullifying the sign of the voltage from one side to the other, thus implementing effective negative capacitive coupling. Additionally, a switch is placed between the capacitor and the transformer to switch between the transformer and the wires going to the neighboring cell to turn the negative coupling on or off.

8.2.2 Hybrid Digital-Analog Method

An alternative way to achieve heteropolar coupling is the hybrid digital-analog method, illustrated in Figure 20. In this method, the hardware has decoupled unit cells, each of which is equipped with a means to measure voltage and is capacitively coupled to an arbitrary voltage source. This voltage source is the voltage of the other unit cells, V = ( v ₁ , v _{2, ...} , v _{N -1} , v _N ) can be controlled to supply voltage f _i ( V ) to cell i . This supplied voltage can be digitally nullified, thus effectively implementing negative coupling.

8.3 Adjustable configuration of inductors

Real inductors have the disadvantage of taking up a lot of space on the chip. Therefore, having a replacement circuit element that acts like an inductor can help reduce the area used on the chip.

Active inductors can be employed. For example, Fig. 21 shows a schematic of an active inductor based on a gyrator, where two Gm cells and a capacitor provide a transfer function that mimics the behavior of inductors within a given frequency range. By modifying the capacitance of the capacitor, the overall inductance can be adjusted.

This active inductor approach has several limitations. These include limited frequency range, low quality factor, increased power consumption, and nonlinearity. Nevertheless, this circuit can be useful for reducing the on-chip area of the inductor.

9 Inductor Removal Strategy

9.1 Introduction

As mentioned above, the analog computational primitives for Gaussian sampling (and others) involve encoding the variables of interest for the application into the degrees of freedom of the circuit unit cells. Correlations between the variables of interest can be encoded via capacitive coupling of the unit cells. However, as discussed in Section 8.2 above, simple capacitive coupling is not suitable for encoding both negative and positive correlations, or bipolar correlations.

One way to solve this problem is to use a transformer in series with the capacitors, wired in such a way that the polarity of the voltage is reversed when the coupling is desired to be opposite to the pure capacitive coupling. However, this solution has the disadvantage that it requires two inductors per coupling element. Also, in the total-to-total connection graph of the unit cells, the number of coupling elements increases proportionally to the square of the number of unit cells , d ² . This is a serious problem because inductors occupy a large physical area when placed on silicon chips , thus limiting the number of unit cells per chip.

To solve this problem, we use a unit cell structure based on RCR circuits, where the degree of freedom is the voltage difference across the capacitors. Then, the heteropolar coupling can be done through the "parallel" and "cross" connection of a pair of coupling capacitors. This device can be used as a sampling processor following the Langevin Monte Carlo protocol (see Section 19 below).

9.2 Unit Cell

A proposal for a unit cell circuit is illustrated in Fig. 22.

Analyzing this circuit using Kirchhoff's current law yields the following equations describing the currents in each branch:

If we sum these equations,

And the difference is as follows

Equation (24) can be expanded using node voltages in the following manner:

Equation (25) can likewise be developed as follows.

To simplify these equations, we introduce the following new set of variables:

Define the time derivatives. Now equation (26) becomes

Equation (27) becomes:

To characterize the asymmetry between the resistors within a unit cell, the asymmetry parameter α _i is defined as follows:

The sum of the unit cell resistances is defined as R _i = R _ia + R _ib . This asymmetry parameter is -1 < α _i It is defined in the range < 1, and in practice can be kept as close to 0 as possible.

Then equation (30) becomes:

Equation (31) becomes:

Finally, substituting equation (34) into equation (33) gives:

The dynamics of this unit cell has one degree of freedom. can be simulated by , and these dynamics are not affected by the asymmetry of the resistors α _i .

9.3 Combined Cells

The capacitive coupling between unit cells of the type described in Section 9.2 is described below. The circuit contains two switches, which operate in only two possible configurations: one for positive coupling and one for negative coupling.

9.3.1 Negative Combination

In the first case, we consider the case where the effective coupling is negative and the capacitive coupling is between similar nodes of the unit cells, for example, from node ia to node ja , and from node ib to node jb . The circuit diagram is shown in Fig. 23. Analyzing this circuit using Kirchhoff's current law yields equations describing the currents in each branch, which are given by for cell i :

And about cell j

Following the analysis of the unit cell, the sum and difference of the current equations are obtained, and these sums and differences are for cell i .

And about cell j

These equations can be expanded for voltages. For cell i , we obtain

For cell j we obtain:

In the above development, the coordinates defined in equation (28) were used for each cell. Also, the resistance asymmetry parameter from equation (32) and the new capacitance asymmetry parameter defined as follows

And the sum of the coupling capacitances is was used.

9.3.2 Symmetric Case

In the symmetric case where α _i = α _j = 0 and β = 0, we obtain the following set of equations for cell i :

And about cell j

9.3.3 Positive Combination

In the second case, we consider the case where the effective coupling is positive and the capacitive coupling is between different nodes of the unit cells, for example, from node ia to node jb , from node ib to node ja . The circuit diagram for this case is shown in Fig. 24. Following the analysis in Section 9.3.1, we obtain the following equations for the sum and difference of current equations in cell i :

For cell j we obtain:

These equations are very similar to equations (44)-(47), with the degrees of freedom of interest being and The sign of the conjunction between the two has changed.

9.3.4 Symmetric Case

In the symmetric case where α _i = α _j = 0 and β = 0, we obtain the following set of equations for the cell:

9.4 Discussion

The above derivation results show that this method actually changes the switch positions in the combined cell, resulting in two degrees of freedom. and It is shown that the capacitive coupling between the two circuit elements is provided. Also, if there is ideal symmetry between the circuit elements, the unwanted degrees of freedom and Complete the disassociation.

10 Hardware Connection Strategies

10.1 Local Connection table Global Connection

As discussed above, the device of the present invention can be used to extract samples from a Gaussian probability distribution. Conveniently, the Gaussian distribution is completely characterized by its first two moments, i.e., the mean and the covariance. As explained in Section 3.6, one can incorporate the mean as a constant offset for each oscillator variable in postprocessing. However, the covariance matrix must be incorporated into the dynamics of the system, which is represented by the capacitance matrix of the physical device. Given a target Gaussian probability distribution of d variables, a d x d capacitance matrix can be used to completely represent the covariance matrix of the target Gaussian probability distribution, which represents the degree of correlation between any two variables. This matrix consists of the d self-capacitances on the diagonal and the two-oscillator coupling capacitances on the off-diagonal.

Building systems of coupled oscillators involves deciding on the degree of connectivity that can be achieved. The degree of connectivity is how much an oscillator in the system is connected to other oscillators. In the previous sections, coupled oscillator systems were assumed to be all-to-all connected. This is the highest degree of connectivity that can be achieved for a given coupling scheme: all oscillators are connected to all other oscillators.

In this coupling scheme, the hardware has about d ² connections between oscillators (the dimension of the problem). This can be a space bottleneck in larger systems. This problem can be solved by reducing the degree of coupling so that only local couplings exist between a set number of oscillators. The degree of coupling can be defined by the maximum number of connections for a given oscillator. For example, a simple coupling scheme is a connection architecture of degree 2, where oscillators are coupled only to their two nearest neighbors.

However, for less than full-to-full connectivity, there is a mismatch problem between the target covariance matrix and the physically implemented capacitance matrix. As an example, consider a system with d = 5 units. The target covariance matrix for this example is:

Using a device with all-to-all connectivity, we obtain the following capacitance matrix:

This fully describes the target covariance matrix. However, when devices with connectivity k = 2 are used, the following capacitance matrix is obtained:

This does not fully describe the target covariance matrix. We define a useful quantity for further discussion of connectivity, namely the diagonal distance. The diagonal distance measures the distance between a matrix element and the main diagonal of the matrix and is defined as follows:

Here, i and j are the column and row indices of the matrix elements, respectively. In the example above with only nearest neighbor joins, k = 1. It is straightforward to extend this concept to arbitrary dimensions and connectivity.

If the hardware has limited connectivity, the impact of covariance mismatch on the behavior of the overall process must be quantified.

10.2 Hardware with local connectivity

10.2.1 Overview

When a specific application is targeted, the properties of the covariance matrix can be directly emulated in hardware. For example, hardware with degree k connectivity can be designed for certain classes of applications that have banded covariance matrices. Also, the accept/reject step of the process can overcome some mismatch between the target covariance connectivity and the hardware connectivity. For example, if the target covariance matrix is not perfectly banded, but the matrix elements decrease exponentially with distance from the diagonal, a process running on banded k -local hardware should still be able to extract high-quality samples.

Figures 25 and 26 illustrate several possible local connection hardware architectures that could be implemented for distributions of appropriate classes. Implementing these strategies in hardware would provide significant savings in on-chip space and complexity, thus opening up the possibility of larger-scale systems while efficiently extracting samples from a particular family of covariance matrices.

10.2.2 Uniform grid

Figure 25 illustrates three different examples of uniform lattices: a square lattice, a hexagonal lattice, and a king graph. For each unit cell, the number of nearest neighbors (i.e., the number of direct connections) is n = 4, n = 6, and n = 8 for these lattices, respectively.

10.2.3 Cluster Graph

Fig. 26 illustrates an alternative method for achieving local connectivity based on the concept of a cluster graph, where unit cells are grouped into clusters of size n _c . For example, for a square cluster, n _c = 4, and for a hexagonal cluster, n _c = 6. There is full connectivity within a cluster, meaning that each unit cell is connected to every other unit cell in the cluster. There is only sparse connectivity between different clusters, as illustrated in Fig. 26 . For square and hexagonal clusters, the number of nearest neighbors is n = 5 and n = 7, respectively. Other types of clusters, such as octagons or other shapes, can also be considered.

Cluster graph connectivity can be applied to, for example, image processing. This is because the covariance matrix for image processing often has a block cyclic structure. The blocks of this block cyclic structure can correspond to clusters of the cluster graph connectivity. By corresponding the fully connected blocks of the covariance matrix to the fully connected clusters of the hardware connectivity, the structure of the hardware connectivity and the covariance matrix can be made to match (e.g., for image processing).

10.3 Semi-local or hardware with global connectivity

The local connectivity described above has the property that the number of nearest neighbors n does not depend on the dimension D. Now consider a connectivity structure in which n increases with D.

10.3.1 Semi-local Connectivity

At one extreme, global hardware connectivity is difficult to achieve in practice. At the other extreme, local connectivity can limit the speedup that can be achieved with analog hardware compared to digital. Intermediate connectivity, neither fully global nor fully local, can be an attractive alternative.

This intermediate connectivity can be called semi-local, meaning that the number of nearest neighbors n increases monotonically as a function of D , but the connectivity is not completely global.

Figure 27 illustrates the quasi-local connectivity between unit cells in a square grid. Each unit cell is fully connected to other unit cells in the same row, and there are sparse connections between rows of unit cells. In this scheme, the total number of connections is increases asymptotically, which is neither linear nor quadratic in D. The number of nearest neighbors, n , is It increases asymptotically.

This strategy essentially means that each cluster is a row and the cluster size is corresponds to the increasing cluster graph connectivity (discussed above). More generally, the cluster size increases as some function of D (e.g., One can also consider other geometries for the clusters, such as rectangular or square clusters (where D increases proportionally or even linearly in D ).

10.3.2 Full Connectivity

Finally, we consider the extreme case of global connectivity, also called full connectivity, where the total number of connections grows asymptotically with D ² , and the number of nearest neighbors n grows with D .

Hardware with full connectivity can be designed for applications with arbitrarily dense covariance matrices. Figure 28 illustrates one possible layout of eight fully connected oscillators following the description in Section 4.

Fully connected designs can provide the greatest expressibility and allow the greatest speedup over digital hardware, but at a high cost in terms of physical space and complexity on the chip.

The method of achieving heteropolar coupling described in Section 8.2.2 can be used to effectively build fully connected hardware with only order D connections.

11 Handling hardware device faults

11.1 Correcting thermodynamic errors

One of the advantages of using thermodynamic hardware to generate samples for HMC is the inherent error resilience of this approach. This is due to the presence of a Metropolis-Hastings (MH) step in which a sample is accepted or rejected based on the energy difference between the proposed sample and the previous sample. Incorporating this step in conjunction with the analog hardware that proposes the sample can be thought of as thermodynamic error correction. The accept/reject step serves to filter out erroneous samples that are proposed due to hardware constraints such as capacitor imprecision or connection constraints, as discussed below. This natural form of error correction is therefore likely to make the thermodynamic approach to HMC resilient to the presence of imperfections.

The MH phase is typically performed in digital hardware in this methodology. This allows any errors accumulated during the analog dynamics to be corrected. Therefore, the combination of the analog-digital interface (during the MH phase) leads to error correction.

In a sense, the MH step serves to keep the analog dynamics on the right track. In this way, while the analog dynamics may deviate from the desired trajectory due to device imperfections (such as imprecision and noise), the MH step serves to bring the position-momentum coordinates back toward the desired trajectory.

One consideration when studying approaches to thermodynamic error correction is the behavior of the acceptance rate (also known as the success rate or acceptance probability) that reflects the quality of the proposed samples. The following subsections present an in-depth examination of this quantity, emphasizing the behavior of the acceptance rate as a function of capacitor imprecision, tolerance and range, and device connectivity constraints.

11.2 Handling capacitor imprecision, tolerances and ranges

As discussed above, the accept/reject phase of the process is a form of error correction for the hardware. If the distribution described by the physical hardware deviates from the target distribution, samples drawn from this hardware distribution will be rejected, and the overall sample quality will remain intact. This leads to a tradeoff between the overall time taken to extract the desired samples and the accuracy or precision of the hardware. This tradeoff can be quantified by computing the average acceptance ratio of the process using imperfect hardware. This error correction means that the hardware can be configured with a more conservative level of tunability while still reliably representing a useful Gaussian distribution.

The device can then be designed with this concept in mind. The number of precision bits and the tuning range of the capacitors and inductors can be chosen by optimizing the acceptance ratio and the constraints of the hardware architecture. For example, Fig. 29 shows the effect of capacitor imprecision on the acceptance ratio based on numerical simulations. Since the acceptance ratio is only a weak imprecision function, it is possible to use capacitors represented by a small number of bits in hardware (e.g., 2 or 3 bits) while maintaining good performance.

An additional consideration during the design phase is the tolerance of the component. This represents the possible deviation of the component's value from the nominal value. These added deviations can be corrected by the acceptance/rejection phase, see Figure 30. However, these tolerances can be taken into account by a full characterization of the hardware before putting it into service.

11.3 Handling Limited Connectivity

Consider the case of limited connectivity, such as the local connection discussed in Section 10.2. The limited connectivity of the device has several consequences when applying the hardware of the present invention to HMC. In particular, as discussed below, the expected improvement over digital methods should be directly related to the connectivity of the device in some applications, such as generating samples from multivariate Gaussians.

As discussed above, the nature of the acceptance/rejection step serves to ensure that the accepted samples represent the desired distribution. One metric to consider is the acceptance probability, and we must be confident that this probability does not decrease significantly.

Another approach to mitigate the possible negative impact of limited hardware connectivity is to potentially mitigate it by classical preprocessing, which is discussed below for the case of HMC for multivariate Gaussian sampling.

11.3.1 Classical preprocessing for local connectivity constraints

To increase the hardware-induced speedup even with limited connectivity, one can develop methods to increase the size of the covariance values within a given bandwidth. Using permutation matrices, one can reorder the combined cells and find an approximately optimal order that maximizes the weights for a given bandwidth.

One way to achieve this permutation is to construct the Laplacian of the weighted graph using the absolute values of the covariance matrix. The second smallest eigenvector of the Laplacian defines the Fiedler vector. The order of the elements of this vector contains an approximate best-fit order for labeling each of the unit cells.

Therefore, in order to best map the covariance matrix values to the hardware, the second smallest eigenvector of the Laplacian constructed from the covariance matrix is computed. The complexity of this operation depends on the characteristics of the matrix. In the worst case, this is ( d ² ) Eigenvalue error in operations can be computed up to . For sparse covariance matrices, this computation should be cheaper.

Another approach to achieve the same goal is to apply a direct cut to sparse a dense matrix, followed by the Cuthill-Mckee process to band the sparse matrix. The direct cut can be applied during the construction of the covariance matrix, setting some elements to zero if they are smaller than some cutoff value. Thus this can be done while computing the matrix itself. Then the Cuthill-Mckee process scales with how many non-zero elements remain in the sparse matrix, i.e. O ( n _nz ), which should be less than O ( d ² ). Thus applying a simple cut followed by Cuthill-Mckee can take a dense covariance matrix to band form in O ( n _nz ) steps. The number of non-zero elements will depend on the cut, which is a parameter to be set by the user.

For a specific example of multivariate Gaussian sampling for a random covariance matrix where elements decrease exponentially with distance from the diagonal, Fig. 31a shows that the acceptance ratio is approximately constant with increasing dimensionality. For the case where the covariance matrix is a dense random matrix, Fig. 31b shows the exponential decay of the acceptance probability with dimensionality when attempting to generate samples from an approximation of the actual covariance matrix with bandwidth k . These two examples highlight the different expected performances of hardware applied to different types of covariance matrices. They also highlight the potential benefit of preprocessing the covariance matrix (as described above) to bring larger elements closer to the diagonal and potentially prevent exponential scaling.

12 Thermodynamic Error Mitigation Methods

Error mitigation can reduce the impact of errors associated with thermodynamic sampling devices. In particular, this error mitigation method targets imprecision errors, which are errors associated with imprecision of analog hardware components. The acronym for this method is THERMIES (THermodynamic ERror Mitigation via Imprecise Ensemble Sampling).

12.1 Univariate Protocol

We begin with a one-dimensional case to illustrate the core concepts, and then consider a more general case (see the supporting documentation: Error Mitigation for Thermodynamic Sampling Hardware). Imagine a Gaussian sampling device that can realize a random variable X whose probability density function is

That is, the device can sample a mean-zero normally distributed random variable whose variance is a multiple of ε . The mean-zero constraint is not restrictive, because bias can be added after sampling with little computational overhead, and the discretization of the variance is a result of the digital encoding of the parameters that tune the device.

In fact, samples can be drawn from distributions with variances that are not multiples of ε , for example, the normal distribution. Suppose we want to sample (here and in the following represents a normal distribution with mean a and (co)variance b ). This distribution is illustrated in Figure 32.

To (roughly) sample this distribution, we perform the following procedure, which is the basis of Thermies' law:

Then, the probability density function of the random variable X is as follows:

Here, the subscript a indicates that this distribution is an approximation of the target distribution f _t . Since both distributions f ₁ and f ₂ have zero mean, f _a also has zero mean. Furthermore, the variance of f _a is

For this example, the f _a distribution is indicated by the dashed line in Figure 32, which fits the target distribution better than either f ₁ or f ₂ .

instead Forming a mixture of Then, the variance of f _a is as follows:

This is arbitrary This means that we can choose w and m to obtain , so , which can be guaranteed, and is called distributed matching. Obtaining an approximate distribution with is problematic, and this problem is discussed again below.

In the univariate case, the imprecision problem can be avoided by rescaling the random variables to have a feasible variance, i.e., some m About {1, 2, ...} A new probability variable that makes it so can be defined, so that error mitigation is not necessary. This can be difficult or impossible in the multivariate case, which is why error mitigation methods are useful, and we now develop a multivariate generalization of the Thermise method.

12.2 Inaccuracy Dependency

It is desirable that the quality of the samples should not vary with the precision of the hardware implementation. Fortunately, the thermise protocol removes the first-order dependence of the approximate distribution on ε as ε goes to 0. Figure 33 shows the L ₁ distance between f _a and f _t with and without thermise error mitigation. For the mitigated case, the slope vanishes at ε = 0, whereas for the unmitigated case, the slope does not vanish at ε = 0. These numerical results provide evidence that the thermise method (or similar methods) reduces or eliminates the sensitivity of the sample quality to hardware imprecision.

13 Non-Gaussian Extension to distribution

13.1 General Strategy

So far we have described analog hardware for sampling from Gaussian distributions. However, there are many other distributions that may be of interest. Therefore, we will now consider strategies for constructing analog hardware that can generate samples from non-Gaussian distributions.

The framework introduced in Section 2.1 applies to any continuous physical system with coordinates { x , p }, a potential energy function U ( x ), and a kinetic energy function K ( p ). For the Gaussian HMC sampler, the second-order potential can be defined as a function of the position coordinate x . However, other potentials (potentials beyond the second-order potential) can be realized, and this is the conceptual basis for extending beyond Gaussian sampling.

The stochastic dynamics of the HMC process guarantees that in the long time limit, the positions and momentum are distributed according to the Boltzmann distribution given in equation (2), and this distribution is true for a well-defined choice of the potential function U ( x ). Furthermore, the distributions for x and p will be independent, and the marginal distribution for x will be is proportional to . Therefore, we can see that by appropriately selecting the potential U ( x ), we can design and implement a probability distribution from which to extract samples.

For the second-order potential, since the Boltzmann distribution over the position coordinates is Gaussian, we sample the position x from Gaussian~exp( -U ( x )). As mentioned above, this second-order potential is a natural physical potential associated with the LC oscillator system. However, if we can modify the shape of the original second-order potential energy function U ( x ) provided by the LC oscillator circuit, we can extract samples from distributions associated with non-second-order potentials.

In the following, we present two approaches to modifying the shape of the potential function.

An approach based on the thermodynamic concepts of Maxwell's Demon; and

An approach based on multibody coupling of LC oscillators.

13.2 Overview of the Maxwell Demon Approach

Section 14 details an approach to sampling from non-Gaussian distributions based on the concept of Maxwell's demon.

Consider a coupled LC oscillator system. A Maxwell demon is a digital or analog device that applies a state-dependent physical force to a harmonic oscillator system. This force is state-dependent in the sense that it depends on the state of the oscillator, for example, it may depend on the voltage across the capacitor of each unit cell. In practice, this force corresponds to a voltage source (or set of voltage sources) that is physically applied to each unit cell.

In simple terms, the voltage source output by the Maxwell demon device artificially generates a new potential energy function U ( x ) for the original position coordinate x of the LC oscillator circuit. This allows the LC oscillator system to evolve according to a new Hamiltonian. As a result, it is possible to program the probability distribution being sampled by appropriately choosing U ( x ).

13.3 Multi-body Overview of the Combination Approach

Section 16 details an approach to design and implement multi-body combinations (e.g., two-body combinations) of unit cells to sample from non-Gaussian distributions.

The multivariate Gaussian distribution can be described by the first two moments, i.e., the mean vector and the covariance matrix. Therefore, a distribution with non-trivial higher-order moments (e.g., third-order moments) is a non-Gaussian distribution. Therefore, it is possible to shift the Gaussian distribution by focusing on generating higher-order moments in the distribution.

Within the HMC framework, this corresponds to generating higher-order coupling terms in the Hamiltonian. Two-body couplings are very common in physical systems, but higher-order couplings (e.g., cubic) are rare. Therefore, these couplings must be generated carefully in physical systems.

At a high level, the reason for achieving two-body coupling in the design of an LC oscillator system is to use two-port elements (e.g. capacitors) to couple the unit cells. This suggests that three-body coupling can be achieved by replacing the two-port elements with three-port elements.

In fact, it is possible to combine three unit cells using a three-port device. For example, a transistor has three ports: a source, a drain, and a gate. This creates a genuine three-body coupling in the Hamiltonian, and thus a third moment in the corresponding probability distribution being sampled. Varicaps (also known as varactors) provide an alternative to transistors for this purpose.

More generally, sets of 3-port elements can be cascaded together to obtain n -port elements, where n ≥ 3. This allows generating true n- field couplings between the unit cells, and hence n- th moments in the corresponding probability distribution. This is discussed in detail in Section 16.

13.4 Additional Interest Distribution

13.4.1 Univariate distribution

Here we provide some mathematical details about the kinds of probability distributions (other than Gaussian distributions) from which samples can be drawn.

First, for simplicity, we consider univariate distributions. Some notable examples include:

exponential distribution;

gamma distribution;

Laplace distribution;

Generalized normal distribution; and

The exponential family.

Only for non-negative numbers, i.e. the domain Consider two distributions that are defined only for , starting with the exponential distribution:

Next, we consider the gamma distribution, which is a generalization of the exponential distribution:

Here, α > 0 is the shape parameter, β > 0 is the rate parameter, and is called the gamma function. By setting α = 1, the gamma distribution is reduced to an exponential distribution.

The Laplace distribution, unlike the gamma distribution, has a support set which is the entire distribution of the "location" parameter and the scale parameter σ is given as follows:

Next, consider a generalized normal distribution with the following probability density:

If β = 2 is chosen, it corresponds to the standard Gaussian distribution. Also, if β = 1 is chosen, it corresponds to the Laplace distribution. Therefore, equation (74) includes Gaussian distribution and Laplace distribution as special cases.

These distributions are special cases of the exponential distribution family. The exponential distribution family can be written as:

For example, for the exponential distribution, T ( x ) = x, η (Θ) = -β, h ( x ) = 1, and A( Θ) = -logβ. Similarly, other distributions can be reconstructed from the general formula in equation (75).

13.4.2 Multivariate Distribution

In general, we can consider multivariate random variables in which higher-order cumulants affect the shape of the distribution in a statistically significant way. For example, if a multivariate random variable has nontrivial third-order cumulants, the probability distribution will be skewed. Higher-order cumulants define different features of the shape of the distribution. Cumulants are discussed in more detail below.

The exponential family discussed above can be extended to the multivariate case as follows:

Here, h ( x ) is a scaling constant, s is a vector called the natural parameter, t ( x ) is called the sufficient statistic, and A ( s ) is the log partition function.

13.4.3 About distribution Cumulent Approach

Section 2.1 discusses circuits that sample from a normal distribution. To proceed with sampling from non-Gaussian distributions, we can consider the problem of sampling from a distribution that is a small perturbation away from the Gaussian distribution. To what small parameter should we apply this perturbative expansion?

The cumulants of the distribution are defined as follows:

It turns out that this is closely related to the moments of the distribution. In fact, the cumulants are polynomial functions of the moments with integer coefficients, and the first two cumulants are simply the mean and variance of the distribution. However, the cumulants turn out to be more useful than the moments of the distribution, since the normal distribution has the special property that k _n = 0 for n ≥ 3. Hence the answer to the question of the perturbation expansion parameter, i.e. k _{≥ 3} .

The formalism for this expansion is known as the Gram-Charlier A series (closely related to the Edgeworth series), which is expressed as:

Here μ and σ ² are the mean and variance of p ( x ), respectively, and D is the differentiation operator. For simplicity, we assume μ = 0 and σ = 1. Then we expand the exponent and collect similar terms. Furthermore, we use the fact that the derivative of the normal distribution is equal to the normal distribution times the Hermite polynomial He _n , and we obtain:

Here, B _n is a Bell polynomial and He _n is a Hermitian polynomial. Now, if we omit the addition constant, the log probability is as follows.

To understand what this means, let's include only the first two amendments.

Here, maintaining only the first two modification terms cannot be valid for all x , since it is possible for the expression inside the log to be negative. However, we assume that k ₃ and k ₄ are small anyway, so these x are far enough away that we can ignore this pathology. For a completely rigorous expansion, we can consider the Edgeworth series. However, as follows, log(1 + x ) Proceed by performing x expansion:

So we see that the ability to add higher-order polynomials to the Hamiltonian allows us to sample from non-Gaussian regions. This comes with the caveat that we must avoid regions where the log(1 + x ) expansion (and the truncation of the series expansion) is invalid. However, in practice, if k ₃ and k ₄ are small, the quadratic force x ² should form a large enough potential barrier to keep the problematic regions out.

14 Non-Gaussian Maxwell's Demon Approach to Distribution Sampling

14.1 Introducing the Maxwell Daemon

Sampling from a non-Gaussian distribution involves changing the potential energy function of the system. For this purpose, we introduce the concept of Maxwell's Demon. Figure 34 illustrates the basic concept of Maxwell's Demon. Historically, James Clerk Maxwell considered an experiment in which an intelligent observer observes gas particles in a two-chambered box and selectively opens a door only for the fast-moving particles, thereby separating the fast and slow particles over time. This decreases the entropy of the gas system over time, but does not violate the second law of thermodynamics because entropy is created elsewhere. The intelligent observer is called Maxwell's Demon.

As will be discussed in detail below, the Maxwell Daemon allows HMC hardware to go beyond Gaussian distributions. As discussed below, there are a variety of ways to physically construct a Maxwell Daemon (MD), including digital and analog approaches.

14.2 Maxwell Demon in the HMC Context

14.2.1 Potential Changing the energy landscape

Performing HMC sampling for various distributions is done by It involves changing the form of the potential energy U ( x ) function in , which is illustrated in Fig. 35.

In Gaussian sampling, this potential energy function U ( x ) is of quadratic form, represented by the covariance matrix ∑ , as shown in the left panel of Fig. 35. The coupled LC oscillator system (described above) is equipped with this kind of quadratic potential, allowing it to draw samples from a Gaussian distribution P ( x ).

Introducing a Maxwell demon device into the conventional LC oscillator framework modifies the Hamiltonian circuit by changing the potential energy landscape with respect to the circuit variables x . An "ideal" Maxwell demon can force the conventional HMC cell according to a suitably smooth potential energy function U ( x ). The dynamical state x of the LC unit cells is fed to the Maxwell demon, which efficiently computes and transmits the result back to the LC cell as a force. This allows adjusting the probability distribution P ( x ) sampled by the HMC. The probability distribution is an exponential function of the potential energy up to normalization as follows:

14.2.2 Voltage Control As a voltage source MD

Figure 36 illustrates that the Maxwell Demon device can be viewed as a voltage-controlled voltage source (VCVS). In this case, the drawing is shown for a single unit cell, but the concept can be generalized to multiple unit cells.

As shown, the voltage x _i across the capacitor of the i -th unit cell is fed to a Maxwell Demon (MD) device, which can be either a digital or analog device. The MD device processes this input by applying some function f to the input. The output f ( x _i ) is fed back to the voltage source, which is applied as the voltage inside the original unit cell. Therefore, the MD device acts as a VCVS, since the voltage output is controlled by the voltage across the capacitor.

14.2.3 Modified for one unit cell Potential

The mathematical analysis of how the MD device modifies the potential energy function begins by considering the case of a single unit cell, as illustrated in Fig. 36.

The MD device is assumed to draw no current from the unit cell. This assumption simplifies the analysis, but can be relaxed through a slightly more complex derivation.

Under this assumption, Kirchhoff's voltage law can be applied to the circuit of Fig. 36 to obtain the following equation:

Here, is the current passing through the circuit, is the voltage across the capacitor, is the output voltage of the Maxwell Demon device.

This equation can be rewritten in Newtonian notation. That is, let x _i = V _i be the position, let p _i = I _i ( m _i / C _i ) be the momentum of mass m _i , and let f _i ( x _i ) be the force. This gives the following equation:

Then, by taking the indefinite integral of the force, we obtain the following potential energy function:

From this equation we can obtain the desired potential energy function U ( x _i ), because Because we can freely choose the function for . Specifically, We can choose , where h _i (x _i ) is any function. Then we arrive at:

So the output of the Maxwell Demon device By appropriately selecting the potential energy function, the potential energy function can be designed and implemented as a desired function.

14.2.4 Modified for multiple uncoupled unit cells Potential

Extending the above analysis to multiple unit cells can be done in two steps. First, we consider the case where capacitive coupling between unit cells is neglected, and then we include this coupling in a subsequent subsection.

With multiple unit cells, one can imagine that the MD device takes as input the full state vector x = { x ₁ , x ₂ , ..., x _d }. It then outputs a voltage vector g(x) . g can act on the full state vector x , and thus effectively couples the unit cells (even in the absence of direct capacitive coupling).

In this case, we obtain the following equations, similar to those described above for the evolution of each unit cell:

The total force can be written as a vector as:

Here is an element is a diagonal matrix with . Then, we can compute the potential energy as the (component-wise) indefinite integral of f(x) with respect to x . Then, the potential energy is given by:

Therefore, by appropriately choosing g ( x ), we can essentially set up the desired potential energy function U ( x ).

For example, we can choose g ( x ) = x - F ^-1 h ( x ). This leads to the potential energy given by:

So, a proper choice of h ( x ) allows us to design and implement the function U ( x ). After all, from the equation P ( x ) = exp(- U ( x )), this means that we can design and implement the probability distribution from which the samples are drawn.

14.2.5 several Combined Modified for unit cell Potential

For multiple coupled unit cells, the analysis is more complicated. Nevertheless, at the conceptual level, there is not much difference from the uncoupled case. The derivation given above for the uncoupled case provides essential intuition about how the coupled case works. Therefore, for simplicity, the derivation of how the potential is modified in the coupled case is omitted.

14.3 Digital Maxwell Demon

Figure 37 illustrates how a digital Maxwell Demon device interacts with an analog LC oscillator system. This fits into the paradigm illustrated in Figure 36, which views the MD device as a voltage-controlled voltage source.

The digital MD device can be stored and processed on the CPU or FPGA. ADC and DAC are employed to convert analog signals and digital signals to each other. Specifically, the voltages across the capacitors (denoted as x ) in the unit cells can be converted into digital signals and fed to the digital MD device. The digital MD device then outputs a digital signal, which is then converted into an analog voltage vector g ( x ), and the components of this vector are applied as voltage sources in the appropriate unit cells.

The main advantage of the digital approach to MD devices is their flexibility and programmability. Digital devices provide the flexibility to choose essentially any function g ( x ) for the output. In this sense, the digital approach to MD devices is programmable and flexible, allowing for sampling from a wide range of probability distributions P ( x ).

A major drawback of the digital approach to MD devices is that it does not fully exploit the potential speedup that can be achieved with analog hardware. Digital devices can compute the slope of the potential energy function, which can be difficult for complex potential energy functions. In contrast, in the analog approach to MD devices, the slope of the potential energy function corresponds to a physical force, which is applied naturally to the LC oscillator system without any computational effort (this is discussed in more detail below). Therefore, there is a potential for computational speedup using analog MD devices.

14.4 Analog Maxwell Demon's Motivation

The Maxwell Demon (MD) device, if implemented in an analog fashion, offers the opportunity to significantly accelerate the computation. This is because the analog MD device can apply forces to the LC oscillator system without performing any computations to compute the forces. This is in contrast to digital MD devices, which compute forces from the potential energy function. This computation is usually expensive and can lead to non-trivial scaling, such as linear or quadratic scaling, depending on the dimensionality of the problem. Therefore, the advantage of the analog MD approach is that there is no need to compute the gradient of the potential energy function.

14.5 Analog Maxwell To the demon About Korea Emergency Approach

14.5.1 Overview

Next, we consider an approach to constructing an MD device that does not correlate the elementary unit cells, before introducing an alternative approach that correlates the elementary unit cells in Section 14.6 below. Each LC oscillator in the HMC hardware discussed above is called a elementary unit cell. In the following, we introduce additional circuitry called auxiliary unit cells.

The uncorrelated strategy for constructing an analog Maxwell daemon involves several building blocks, which are illustrated in Figure 38 and described as follows:

1. Voltage measurement of the basic unit cell capacitor, which is applied as a voltage source to supply current through the auxiliary unit cell.

2. Nonlinear circuit elements that allow complex I-V relationships in auxiliary unit cells.

3. Measure the voltage across the resistor of the auxiliary unit cell (to read the current in this cell), which is then inverted and applied as a voltage source to the primary unit cell.

Figure 38 illustrates these components in a single unit cell. In reality, the same five structures can be replicated d times for d unit cells, possibly with different hyperparameters (e.g., different nonlinear elements) for each copy.

The voltage output by the auxiliary unit cell is reversed before being applied to the primary unit cell. This is because the system must have negative feedback (instead of positive feedback) to maintain stability. Negative feedback provides a restoring force, similar to how a spring in physics pulls a mass back toward its equilibrium point. In short, the voltage reversal helps maintain the stability of the system.

However, Fig. 39 shows that voltage reversal can be prevented by disconnecting the resistor in the auxiliary unit cell from ground. This allows physically reversing (or swapping) the wires that carry the voltage difference across this resistor as shown in Fig. 39.

There are a number of additional components that can be inserted into the circuit to increase flexibility. These optional additional components include:

An amplifier with adjustable gain inserted into the coupling between the two unit cells. For example, the amplifier can amplify the voltage read from the capacitor of the primary unit cell to boost the voltage source applied to the auxiliary unit cell. Alternatively, the amplifier can amplify the voltage read from the resistor of the auxiliary unit cell to boost the voltage source applied to the primary unit cell.

An additional voltage source inserted into the basic unit cell, which will be constant (independent of the state variables) but adjustable. This provides flexibility by providing a constant offset to the voltage output by the Maxwell Demon, and a constant offset to the force applied to the basic unit cell.

14.5.2 Auxiliary Unit Cells

As illustrated in FIGS. 38 and 39, the basic components of the auxiliary unit cell are a voltage source (provided from the basic unit cell), a nonlinear element, and a resistor.

One strategy is to use the IV characteristics of a nonlinear element (NLE) as a function generator. The input V _i and the output We want to generate a complex functional relationship between, for example, as mentioned in the discussion of equation (81). The IV property of NLE provides one way to design and implement this complex functional relationship. Thus, the desired function The problem of generating is transformed into the problem of designing and implementing the desired IV characteristics for an NLE.

In the following, we discuss explicit strategies for constructing NLEs based on individual circuit elements or on combinations of elements.

14.5.3 Single peak Nonlinear elements for distribution

Nonlinear elements (NLEs) can be used to extract samples from a unimodal distribution. An NLE with a unimodal I-V characteristic mathematically generates a unimodal probability distribution. Therefore, we will only consider NLEs with a unimodal I-V characteristic here.

Let I _i ( V _i ) be the current through the ith NLE as a function of the voltage from the ith elementary unit cell. The output of Maxwell's demon is , which can be set equal to the force f _i ( V _i ) acting on the i-th basic unit cell for simplicity. This assumption is This can be justified by assuming that is large compared to V _i , as shown in equation (82) (where x _i is used instead of V _i ). The resulting potential energy U _i ( V _i ) is then the negative integral of the force f _i ( V _i ), and the associated local probability distribution P _i ( V _i ) is then the negative exponent of the potential energy. To simplify notation, the current, force, potential energy, and probability distribution can be written as I ( x ), f ( x ), U ( x ), and P ( x ), respectively.

Figure 40 shows a schematic diagram of these four functions for the special case where the NLE is a forward-biased diode. A diode is a nonlinear device that conducts current roughly in one direction, with the current increasing exponentially with the forward-bias voltage. However, there is a breakdown voltage that allows current flow in the opposite direction. This gives rise to a somewhat asymmetric IV characteristic, as shown in the first panel of Figure 40. The second panel shows the force f ( x ), which is the negative of I ( x ) (due to the voltage inverter shown in Figure 38). The third panel shows U ( x ), and the fourth panel shows P ( x ). Since P ( x ) is asymmetric, it is non-Gaussian and has a steeper decay than Gaussian.

Instead of considering a single diode, we consider two diodes oriented in parallel with opposite directions. Choosing these two diodes in parallel with the NLE will result in a symmetric distribution P ( x ), since the circuit elements are symmetric with respect to positive and negative voltages. Choosing the NLE makes it possible to design and implement a symmetric version of the probability distribution shown in Fig. 40.

Transistors also exhibit nonlinear IV characteristics. For example, in the common-base configuration, the input characteristic is a convex nonlinear function, while the output characteristic is a concave nonlinear function that saturates for large voltages. The former case (the input characteristic) produces a distribution P ( x ) similar to that of a diode, while the latter case (the output characteristic) produces a distribution P ( x ) with a long, slowly decaying tail. Thus, the latter allows for non-Gaussian distributions with long tails. Again, it is possible to make these distributions symmetrical about zero voltage by orienting the two transistors in parallel in opposite directions. Other configurations for transistors are also possible, such as common-emitter.

14.5.4 Multi mode Nonlinear elements for distribution

We now discuss NLEs for sampling from multimodal distributions. To obtain multiple modes, the potential energy must have multiple local minima (e.g., local and global minima). This means that the derivative of the potential energy, i.e., the force, must change sign more than once. If the force changes sign more than once, it must be a nonmonotonic function. Therefore, the I-V characteristic of an NLE must be a nonmonotonic function (since the property is the negative of the force). Therefore, we discuss NLEs with nonmonotonic I-V characteristics.

First, consider the tunnel diode as an NLE. In this case, Fig. 41 shows the associated curves for the current, force, potential energy, and probability distribution. The current looks somewhat cubic, is clearly not monotonic, and has local maxima and minima in the first quadrant. Then, assuming that we shift the force by a constant negative amount, we obtain a force f ( x ) that changes sign many times. This is achieved by subtracting a constant voltage from the output of Maxwell's demon, which physically corresponds to applying a constant negative voltage source to the basic unit cell. The resulting f ( x ) curve is shown in the second panel of Fig. 41 . This gives rise to a potential energy U ( x ) with two local minima, and consequently to a multimodal probability distribution P ( x ) with two modes.

As an alternative to the tunnel diode, the Gunn diode is also considered as a NLE candidate. The Gunn diode also has a nonmonotonic IV characteristic under forward bias, with local maxima and local minima in the first quadrant. It also has a nonmonotonic IV characteristic under reverse bias. Therefore, the curve is nonmonotonic for both positive and negative voltages. This can lead to a more complex probability distribution P ( x ) than that of the tunnel diode. In particular, for the Gunn diode, P ( x ) can have multiple modes in both the positive and negative voltage regions.

Other NLEs such as memristors can also be considered. In addition, complex I-V characteristics can be generated by combining various nonlinear elements in series or in parallel. For example, when NLEs are combined in series, the overall I-V curve becomes a composite of the individual I-V curves. In this way, complex I-V curves can be designed and implemented by combining various circuit elements to form an overall NLE block.

14.6 Analog Maxwell To the demon Korean Correlation Approach

14.6.1 Overview

We now consider an approach to constructing an analog Maxwell demon involving correlation of multiple unit cells. The outline of this strategy is depicted in Fig. 42 and includes the following building blocks:

1. Voltage measurements of each basic unit cell capacitor whose collection forms a voltage vector v .

2. An analog coupling operation, where the elements of the voltage vector v are coupled together, so that the output of the coupling operation is another vector h(v) . This coupling may correspond to an affine transformation h ( v ) = A v + b , or it may correspond to a higher-order tensor, such as a second- or third-order form.

3. A nonlinear operation acting on each component of the vector h ( v ). This may involve driving a current through a nonlinear circuit element using the vector components of h ( v ), and then reading the current through this nonlinear element (NLE) by selecting the voltage across a resistor connected in series with the NLE. Examples of NLEs, including various types of diodes and transistors, are given in Sections 14.5.3 and 14.5.4. The resulting voltage vector is written as n ( h ( v )), where n corresponds to the nonlinear operation.

4. The previous two steps can be repeated multiple times so that there are multiple layers of combinations and nonlinear operations. If there are L layers, the final output vector will have the following form:

Here, the symbol "o" indicates that the operations are composed together sequentially. Here, h _j and n _j are the coupling and nonlinear transformations for the j -th layer, respectively.

5. The components of the v _output are inverted to achieve negative feedback, the reason being explained in Section 14.5.1. This can be done using a voltage inverter or by disconnecting the resistor of the auxiliary unit cell from ground and swapping the wires associated with the voltage difference across this resistor.

6. After inverting the signs of the elements of the v _output , these elements are applied as voltage sources inside the corresponding basic unit cells.

In the following, we consider various methods of constructing the analog combination operation shown in Fig. 42.

14.6.2 Analog Affin Conversion

The analog combination operation mentioned above can take various forms, including, for example, affine transformations. Fig. 43 provides an analog circuit implementing an affine transformation for a vector v . Suppose we want to implement the A v transformation for some matrix A . Each row of A is multiplied by v , and this multiplication can be accomplished using layers of resistors as illustrated in Fig. 43 . Assuming that A is a dense matrix, this involves d ² resistors, since there are d resistors associated with each row of A and d rows of A . If A is not dense, fewer resistors can be employed. As an optional feature, adding a voltage vector b to the output yields the general affine transformation A v + b .

14.6.3 Higher-order combination operations

The analog affine transformation mentioned above can be viewed as a first-order associative operation, since it produces a weighted sum of terms that are linear in the components of v .

Alternatively, one could consider analog combination operations that include higher-order terms, such as quadratic or cubic terms, in the components of v . One advantage of this approach is that it allows representing more complex probability distributions P ( x ), including distributions that are difficult to sample using standard digital methods.

One way to obtain higher-order terms is to replace each resistor in Fig. 43 with another circuit element, say an n- port element where n > 2, to achieve higher-order couplings.

For example, consider a three-port device, a transistor. Figure 44 illustrates one way to create higher-order couplings using transistors. Here, a voltage component v _j can be used as a gate voltage for a transistor in the k -th arm of the circuit, where k ≠ j . In this sense, v _j controls the resistance encountered by the voltage v _k to create higher-order couplings.

More generally, other higher-order couplings can be created using other multiport devices, including chains of multiple transistors.

14.7 Maxwell's Demon Approach to Distributions Represented by Neural Networks

In many applications, the log probability being sampled is proportional to the output of the task. In this case, the probabilities can be expressed as:

Specifically, the posterior probability distribution sampled in Bayesian inference is

, which is the prior probability distribution Wow, too is defined as the product of , is a set Conditioned on the observations of is the probability of observing the gradient. In this case, the prior probability distribution and/or the likelihood can be parameterized by the neural network. Here, the runtime of the digital HMC can be limited by the computation of the gradient, which is typically computed using automatic differentiation. The runtime of the gradient computation in the neural network depends on the number of parameters and the number of layers in the network. For example, for a dense feedforward neural network with the same number of neurons in each layer m , the runtime of the gradient computation using automatic differentiation is O ( dm2 ). In many cases ^, m can be chosen to scale with the input size, and ^thus the runtime for neural networks employed in modern applications is O ( d3 ).

In this case, the model If this analog (e.g. using an analog neural network as described above) is compiled and its gradient can be estimated, this bottleneck can be eliminated by using analog hardware. Such a device is hereinafter referred to as a "gradient calculator" (4114). The gradient calculator (4114) takes as input a voltage vector corresponding to the output of the analog neural network device and a voltage vector corresponding to the negative gradient. is an analog device that outputs a gradient calculator (4114). One possible implementation of the gradient calculator (4114) includes an analog adder and subtractor to estimate the gradient using finite differences. Solving the Hamiltonian equations uses the output voltage vector of the gradient calculator (negative slope) can be performed in hardware by feeding the d integrator units. In this way, the voltages obtained at the integrator output represent the momentum variables (momentum integrator unit), and this integration step corresponds to the Hamiltonian equation for momentum:

Therefore, after some time T , we obtain:

Here we explicitly show the time dependence of the position, and are the resistors and capacitors that characterize the momentum integrator device. Then the position is given by the integral of the momentum, which is the voltage vector This can be done by feeding the second integrator device (position integrator device):

Inverse mass matrix after time increment T The multiplied voltage vector is fed to the gradient calculator. In most cases the mass matrix is diagonal, so this corresponds to multiplying the voltage vector by a constant value, which can be done in an analog fashion using an affine layer consisting of d resistors, as illustrated in Fig. 43. The output can also be rescaled using the 1/ RC term to obtain the exact position. Here all values of R and C can be set to the same. If the components are analog and there is no latency between the gradient calculator, the momentum integrator and the position integrator, then the integration of each of the Hamiltonian equations will be exact (up to noise sources such as the Johnson-Nyquist noise discussed above).

Figure 45 illustrates a schematic diagram of a Maxwell daemon device (4110) for extracting samples from a distribution represented by a neural network. The user (41) first compiles into a hardware analog neural network (ANN) (4112) (including gradients if closed form exists). Then, the gradient output is fed to a momentum integrator unit (4120), and the output of this unit is then fed to a position integrator unit (4130). The momentum integrator unit (4120) is a unit composed of d analog integrators that output voltages representing momentum as illustrated in FIG. 45. The position integrator unit (4130) is also composed of d analog integrators and represents position variables. A voltage vector corresponding to the position can be fed to an affine layer, which multiplies the voltage vector by M ^-1 . If M is arbitrarily dense, the inversion can take O ( d ³ ) steps. The matrix inverter block (4140) is an affine layer as shown in the diagram, where the elements of matrix A are elements of M ^-1 .

14.7.1 Slope calculation

Computing the gradient in an analog way is difficult. A simpler way to estimate the gradient of a function is to perform finite differences. This can introduce noise into the gradient computation. In the case of HMC, since there is a Metropolis-Hastings step, this noise does not degrade performance, as it would be another source of uncorrelated noise, as discussed below for stochastic gradient HMC. Specifically, the finite difference gradient can be written as the exact gradient with added Gaussian noise of zero mean and variance V :

In this way, slopes can introduce errors in the dynamics. These errors can be considered as an additional source of noise added to the noise coming from the hardware.

15 Through the precision matrix Gaussian Hardware for sampling

The device depicted in Fig. 45 can be slightly modified to perform Gaussian sampling. One difference is that the slope calculator uses voltage , which feeds , that is, the negative slope of the multivariate Gaussian distribution. In some exemplary applications, such as Bayesian regression, the matrix Σ ^-1 = Q , known as the precision matrix, is provided directly, because it participates in the process. In fact, the multivariate Gaussian can be generally defined by the mean and the precision matrix.

Fig. 46 illustrates an alternative arrangement for sampling a multivariate Gaussian using a given mean and precision matrix. The slope calculator can be replaced by a resistive layer (4210) whose resistors correspond to the elements of the matrix Q. To achieve this, it is possible to add an operational amplifier to each output port of the resistive layer so that the value does not disappear. Then, the negative slope The voltage corresponding to can be fed to the parallel integrators to form a momentum integrator device (4120), the output of which corresponds to the momentum. As above, the following is established:

and

This has the quantities defined in Section 14.7. As before, the outputs for both the momentum integrator unit (4120) and the position integrator unit (4130) must be rescaled by the 1/ RC term. These units are completely parallel; they are not coupled. If, i.e., the operations of the gradient calculator (4210), the momentum integrator unit (4120) and the position integrator unit (4130) are analog and there is no latency between them, then the integration of each of the Hamiltonian equations should be accurate (up to noise sources such as the Johnson-Nyquist noise described above). The advantage of this approach is that it is scalable since no inductors are required. It is not limited by the coupling constraints as in the case of the coupled LC resonator unit disclosed above. In fact, this hardware can handle dense covariance or precision matrices without all-to-all coupling.

16 Non-Gaussian For distribution sampling Multi-body Combined approach

16.1 Motivation and Overview

16.1.1 2 bodies Beyond the Union

The previous section described one way to go beyond the Gaussian distribution based on the adoption of Maxwell's Demon (MD). This approach assumed that the direct coupling between unit cells is a two-body coupling, like the capacitive coupling discussed above.

In this section, we describe an alternative approach beyond Gaussian distributions that involves changing the direct coupling between unit cells. Specifically, this approach designs and implements the direct coupling between unit cells to be a multi-body coupling, such as a 3-body or, more generally, an n- body coupling.

This approach gives rise to n- body terms in the potential energy function U ( x ), which in turn leads to higher-order cumulants in the probability distribution p ( x ), for example k ₃ or k ₄ , as discussed in Section 13.4.3. There is a mathematical connection between introducing many-body coupling and introducing higher-order cumulants into the distribution p ( x ), as discussed in Section 13.4.3.

To construct 3-body couplings, we use electrical components with three ports (instead of the two ports associated with capacitive coupling). More generally, to construct n -body couplings for some integer n > 2, we use electrical components with n ports. The voltages associated with n different unit cells (i.e., n different LC oscillators from the HMC hardware) are fed to n ports of the electrical component responsible for the coupling. This n- port component is called the coupling element.

Below we detail specific physical circuit elements that may be used to construct an n- port coupling device. For example, networks of transistors or networks of varicaps (also known as varactors) may be employed for this purpose, as detailed below.

This n- body coupling approach can be used independently or in conjunction with the Maxwell demon approach. Also, this approach can be used independently or in conjunction with the two-body capacitive coupling described in the previous sections.

Changing the direct coupling between unit cells has the advantage that it is difficult to simulate digitally, as detailed below.

16.1.2 Effective Multi-body Digital simulation of bonding

In principle, many-body couplings can be difficult to simulate on a standard digital computer. One way to understand this is that a typical k- th order coupling is mathematically represented by a k ^-th order tensor with dk elements. Therefore, computing the coupling at a given point in time may involve dk multiplication operators. For example, a cubic coupling would involve ^d3 operations ^that a digital computer can compute. Therefore, many - body couplings like this can be difficult to simulate on a digital computer.

However, there are situations where digital devices can efficiently simulate such couplings. That is, consider the case where a many-body coupling is generated by a series of two-body couplings. In this case, the many-body coupling is not real but effective, meaning that many many-body couplings can be efficiently represented by a small number of two-body couplings. (In contrast, a true many-body coupling is a coupling that cannot be efficiently represented by a small number of two-body couplings.)

The Maxwell Demon approach described in Section 14 is largely applicable to the case of valid many-body couplings, but not true many-body couplings. One exception is the approach described in Section 14.6.3, which involves transistors to create higher-order couplings that are actually true.

If we ignore this exception and instead consider the analog affine transformations of Section 14.6.2, we expect that a digital computer will be able to simulate the analog Maxwell's demon with a complexity that scales as d ² , since the matrix-vector multiplication of the affine transformation involves d ² operations. This d ² complexity is maintained even if the analog Maxwell's demon generates many-body couplings via nonlinear circuit elements, since this version of the analog Maxwell's demon involves valid many-body couplings, not true many-body couplings.

When aiming to design analog systems that are difficult to simulate digitally, it is helpful to use true many-body coupling, rather than valid many-body coupling. Here, when we say "difficult to simulate digitally", we mean that the complexity of employing a digital computer scales by more than d ² (e.g., d ³ scaling).

16.1.3 True Multi-body Combine

Let us briefly explain what we mean by true many-body coupling. True many-body coupling can be viewed as a many-body coupling of the relevant circuit elements. This is in contrast to many-body coupling, which consists of circuit elements with few-body coupling (e.g., two-body coupling).

From a computational point of view, a true many-body coupling can be thought of as a case where constructing a many-body coupling does not entail additional computational resources. For example, one can combine several two-body couplings with auxiliary systems to form a three-body coupling, but this entails using additional resources (e.g., auxiliary systems) and is therefore not considered a true three-body coupling. On the other hand, it is possible to construct a three-body coupling with only a single transistor and without additional auxiliary systems, which is considered a true three-body coupling.

16.1.4 Capacitance High resistance Multi-body Combine

As mentioned above, creating a multi-body coupling element requires employing circuit elements with many input and/or output ports.

One goal is to design and implement a probability distribution P ( x ) from which to draw samples. This corresponds to designing and implementing a potential energy function U ( x ), Because there is a relationship called .

From the perspective of physical circuit elements, potential energy is associated with capacitive elements in the circuit, since capacitors can store energy. On the other hand, resistive elements in the circuit are associated with damping or friction, which technically do not directly contribute to the potential energy function.

Therefore, for designing and implementing P ( x ), we are more interested in capacitive elements than resistive elements. Purely resistive circuit elements are not very useful for designing and implementing probability distributions P ( x ).

Therefore, in addition to wanting n- port devices, we also want them to have some capacitance. In particular, we want that capacitance value to be voltage controlled, i.e., controlled by the other input ports of the device. In this sense, we are interested in n- port devices that essentially act as voltage-controlled capacitors.

16.1.5 Elements that act as voltage control capacitors

There are several approaches that can be taken to construct a voltage controlled capacitor.

Transistors have capacitance that can be affected to some extent by the gate voltage. Therefore, transistors are one option for voltage-controlled capacitors. A strategy must be considered to ensure that none of the three terminals of the transistor are grounded and consequently do not take up too much space inside the chip. Fortunately, transistor technology based on Fully Depleted Silicon on Insulator (FDSOI) provides a way to keep the space inside the chip small in this case.

Varicaps (also known as varactors) provide another option. Although each varicap is a two-port device, two varicaps can be placed in a back-to-back configuration facing each other, as illustrated in Fig. 19. By inserting a voltage line between the two varicaps, the overall capacitance of the system can be adjusted. The device illustrated in Fig. 19 acts as a voltage-controlled capacitor with a total of three ports. It thus has the desired characteristics for the purposes of this application. It is referred to as a three-port varicap device.

Another strategy is to use a switched capacitor bank, as shown in Fig. 18. An external voltage signal controls the capacitance of this capacitor bank, and it can therefore be viewed as a voltage-controlled capacitor. One drawback of this approach is that the circuit devices of Fig. 18 need to be turned on and off in real time, which may result in energy loss (i.e., attenuation) during the HMC process.

Transistors, 3-port varicap elements, and switching capacitor banks can all be viewed as voltage-controlled capacitors. Each of these voltage-controlled capacitors can act as a building block of higher-order coupling. That is, by chaining sets of these 3-port elements, an n -port element can be formed. These n- port elements then enable the design and implementation of an n- body coupling.

16.1.6 Adding Tunability or Switches to a Bond

Adding tunability to 3-body or n -body combinations would allow hardware devices to represent a wide variety of probability distributions with n -body terms.

Adding tunability to an n- body combiner can be challenging, but a simpler approach is to add switches to the combiner. The switches turn the n- body combiner on or off. The user can specify the probability distribution by choosing the direction of these switches. This provides a coarse-grained degree of tunability that is relatively easy to implement.

16. 2 years old Dog's Combined unit cell

Figure 47 illustrates a schematic diagram of a three-body coupling mediated by a voltage-controlled capacitor.

Here, we extract the voltages associated with each capacitor (with respect to ground) from each of the three unit cells. These three voltages are fed as inputs to the three ports of the three-body coupler. For example, if the three-body coupler is a transistor, one voltage is fed to the gate of the transistor, one voltage is fed to the source, and one voltage is fed to the drain. We denote the three-body coupler as a voltage-controlled capacitor (VCC), and assume that there are various physical implementations of this VCC (as discussed above).

In the following subsections, we analyze how this combination affects the probability distribution being sampled.

16.2.1 Mathematical description of three combined cells

For two capacitively coupled unit cells, the contribution to the potential energy has the form U = V ₁ V ₂ C ₁₂ . Now assume that C ₁₂ is a function of the voltage V ₃ of the third unit cell, which is the case depicted in Fig. 47. Then the contribution to the potential energy is U = V ₁ V ₂ C ₁₂ ( V ₃ ), where C ₁₂ ( V ₃ ) is a function determined by the physical properties of the voltage-controlled capacitor.

For example, suppose that this is a linear function that is often true over a small voltage range. In this case, Then, the contribution to potential energy is:

So in this case we get two terms in the potential energy: (1) a three-body term proportional to C ₁₂₃ , and (2) A two-body term proportional to .

In Section 13.4.3, we established a mathematical connection between the three-body terms in the potential energy and the third-order cumulants (denoted by k ₃ ) in the corresponding probability distribution P ( x ). Thus, equation (98) generates a probability distribution P ( x ) that is non-Gaussian with a third-order cumulant. Therefore, we can design and implement the probability distribution P ( x ) to be non-Gaussian using the three-body coupling approach.

16.3 Combining four or more unit cells

Figure 48 illustrates a schematic of how four unit cells can be combined using a four-body combiner. The voltages of the third and fourth unit cells, V ₃ and V ₄ , are fed to a voltage multiplier, the output of which is the product of V = V ₃ V ₄ . The voltage multiplier may be analog or digital, as described in detail below. This voltage V is then applied to the gate of a voltage-controlled capacitor, which induces a capacitance C ₁₂ ( V ₃ V ₄ ) between unit cells 1 and 2.

C ₁₂ ( V ₃ V ₄ ) can be a complex function, but this Assuming that it is a linear function over some domain of the form, in this case we obtain a contribution to the potential energy of the form:

Thus, through this implementation we obtain both four-body bonding and two-body bonding.

In Section 13.4.3, we established a mathematical connection between the four-body terms in the potential energy and the fourth-order cumulants (denoted by k ₄ ) in the corresponding probability distribution P ( x ). Hence, using this approach, we obtain a non-Gaussian probability distribution P ( x ) with a fourth-order cumulant.

The approach illustrated in Fig. 48 can be straightforwardly extended to combine more than five cells ( n > 4). This involves feeding an additional input to the voltage multiplier, such that the n -2 voltages are multiplied together to obtain the voltage ultimately applied to the voltage-controlled capacitor gate of the n- body coupler. Due to the mathematical connection between n- body combinations of potential energy and n -th order cumulants of the probability distribution P ( x ), the resulting distribution is non-Gaussian and has an n- th order cumulant (including other possible lower order cumulants).

16.4 Multi-body Analog vs. Digital-Analog Approach to Coupling

16.4.1 Fully Analog Approach

In principle, the components shown in Figs. 47 and 48 can be analog components, since both the voltage control capacitors and the voltage multipliers can be constructed from analog components. This has the advantage of low latency, since no communication with a digital device is involved. Furthermore, this can lead to greater speedups, since multiplying voltages incurs some computational overhead in a digital device (but not in an analog device).

However, the fully analog approach has more circuit elements than the two-body coupling case. For example, in the three-body coupling case, these circuit elements include the infrastructure associated with wiring the third unit cell to the inputs of the capacitors of the two other unit cells. In the four-body coupling case, there is also infrastructure associated with the voltage multiplier. Therefore, the fully analog approach can increase the circuit complexity.

16.4.2 Digital-to-Analog Approach

A convenient approach to multibody coupling is to use existing infrastructure that is already in place to implement HMC hardware that exploits two-body coupling. That is, Section 8 describes the use of voltage-controlled capacitors in the two-body case, where digital hardware is used to set the capacitance value.

In the case of a 3-body combination, the voltage measurement V ₃ can be taken from the third unit cell and sent as an input to the digital hardware, which outputs a voltage proportional to V ₃ and applies this voltage to the capacitor connecting unit cells 1 and 2. Thus, the digital hardware acts as an intermediary between the capacitor between unit cells 1 and 2 and unit cell 3.

For the four-body combination, the same approach can be taken, except that the digital hardware computes the product of the voltages V ₃ and V ₄ . The digital hardware outputs a voltage proportional to V ₃ V ₄ and applies this voltage to the capacitor connecting unit cells 1 and 2.

Again, one advantage of this approach is that it leverages infrastructure that already exists to implement HMC using two-body coupling.

16.5 Local Connection Strategy

In principle, hardware can be designed to allow full connectivity (global connectivity). The number of connections increases rapidly in this case, since n- body coupling leads to d ⁿ scaling with respect to the total number of connections. There is a tradeoff to be considered between hardware complexity and speedup, with more connections potentially leading to greater speedup, but also more complex hardware.

This tradeoff may make it desirable to consider local connectivity. One approach to local connectivity is to leverage the existing infrastructure used for two-body coupling. That is, we can adopt a digital-to-analog approach to many-body coupling (discussed immediately above), and use whatever interconnect structure is already in place for two-body coupling as the interconnect structure for n- body coupling. In this case, the interconnect structure of the n -body coupling will be identical to that of the two-body coupling. So if the two-body coupling is local, then the n- body coupling is also local.

An alternative approach is to design an n- body connection structure from scratch without utilizing a two-body connection structure. Here, we present several strategies for this approach for the special case of a three-body connection, as illustrated in Figures 49 and 50. Figure 49 shows how to fit a three-body connection into a square lattice, and Figure 50 shows how to fit a three-body connection into a hexagonal lattice.

17 Speed Improvement Analysis

17.1 Digital HMC’s Runtime analyze

An analysis of the computational complexity of the HMC process can be found in Neal et al. , "MCMC using Hamiltonian dynamics," in Handbook of Markov Chain Monte Carlo , 2.11 (2011), which is hereby incorporated by reference in its entirety for all purposes. As noted, the system performs HMC by integrating the Hamiltonian equations as follows:

These are vector differential equations, and therefore and There are d equations to solve using d- dimensional vectors.

A numerically common integration method is the leapfrog integrator, which provides updates to position and momentum as follows:

In the optimal case, the numerical complexity of integrating these equations is given by O ( dn _t ), where n _t is the number of time steps. However, computing the right-hand side of each equation to perform numerical integration can be more expensive. In practice, computing the log-probability depends on the shape of the probability distribution being sampled. The mass matrix M is in the general case d × d matrix, and therefore the cost of computing the inverse matrix scales by d ³ steps. However, in most implementations, M is considered to be a diagonal matrix, and the cost is reduced by d .

In order to keep the numerical integration error sufficiently low without being tied to a specific theory, the time step size used for numerical integration should be scaled as d ^1/4 . Thus, the gradient of the log probability is obtained efficiently, and we reach a numerical scaling of d ^5/4 by assuming that the mass matrix M is diagonal and therefore its inverse is computed in d steps. (In fact, if the probability distribution were represented by a neural network, the complexity of finding the gradient would grow polynomially with the number of parameters in the model and the size of the input.)

In this sense, the d ^5/4 scaling is the optimal scaling for digital implementation of HMC. In cases where the sparsity structure of M is not known, dense matrix inverse computations such as Cholesky factorization, which costs O ( d ³ ) steps, can be used. In the following sections, we examine the specific speedup that can be achieved using the hardware of the present invention.

17.2 Analog HMC Runtime

In general cases, the speedup obtained from implementing differential equations in an analog system comes from two main areas: (i) there is no need to discretize the equations and therefore no step size, which means that the numerical error introduced by the step size disappears. (In electrical hardware, there may be additional sources of noise that can affect performance, such as thermal and shot noise.) (ii) the differential equations are solved numerically, for example using frog-hop integration to compute the update terms. In analog, the system naturally follows the differential equations, so there is no need to perform matrix multiplications and inversions, which greatly reduces the number of operations.

An additional subtlety is the gradient calculation. Depending on the shape of the probability distribution, this can affect the speedup achieved by analog HMC devices. This is discussed in detail in Section 14.7 above.

In the ideal case, the number of digital operations to implement an analog SDE solver for a d- dimensional object is O (1). However, one must communicate digitally with the analog solver, for example, to set inputs and measure outputs. For sufficiently small data (tens of megabytes), the data can be loaded into the analog system in O (1) steps. However, as the data size increases, there is a linear scaling with the dimensionality d . Therefore, we believe that the best-case scaling for an analog solver is O (1).

17.3 Gaussian Increased sampling rate

If the probability distribution being sampled is multivariate Gaussian as suggested above, the form of the slope can be computed analytically. Recall that the probability distribution is defined as follows:

The slope of the log odds can be expressed in affine form as follows:

Therefore, computing the gradient is equivalent to performing a vector-matrix multiplication, where the relevant matrix is ^the inverse of the covariance matrix. This shows that the running time of digital HMC for Gaussian distribution is O ( d3 ).

However, if the distribution is multivariate Gaussian, there is a more efficient sampling method known as Cholesky sampling. The trick with Cholesky sampling is to It is parameterized as (random variable) (recall that the mean is 0), but The covariance matrix can be written as follows:

Here the Cholesky decomposition is identified as the last term. Therefore, can be sampled directly at random from (since u _i are independent) and then multiplied by the Cholesky factor L after sampling. can be sampled from. Therefore, Cholesky sampling is generally at least as efficient as computing the full inverse of the covariance matrix. This depends greatly on the structure of the covariance matrix.

17.3.1 Speedup of dense covariance matrices

For dense covariance matrices, Cholesky decomposition is performed in O ( d3 ) steps. In many applications ^, such as Gaussian processes, the covariance matrix is updated and resampled many times, thus changing the structure of the covariance matrix. This means that in many applications ^, the covariance matrix is assumed to be dense. In this case, a speedup of O ( d3 ) can be achieved if the analog system is close to fully connected, i.e., the connectivity k scales as d . Otherwise, when updating the capacitance values so that the capacitance matrix matches the covariance matrix, most of the values in the covariance matrix are suppressed due to the limited connectivity, and thus the acceptance rate of the HMC proposal may be scaled unfavorably. If the covariance matrix is dense but the difference between the largest and smallest elements is many orders of magnitude, the dense covariance matrix can be sparse by a preprocessing step. This step requires O ( d2 ) operations, and can obtain sparse matrices and latent band matrices. This step can be performed not ^only in digital mode but also in analog mode, and can be a bottleneck. Therefore, for the case of connectivity k∼O (1), using analog HMC for Gaussian sampling does not yield any speedup.

17.3.2 Speedup for diagonal and sparse covariance matrices

If the covariance matrix is diagonal, or equivalently the random variables If are uncorrelated, the cost of computing the inverse is simply O ( d ), since the covariance matrix is equivalent to a d -dimensional vector. Drawing samples from this distribution is equivalent to drawing samples from d independent Gaussian distributions, and therefore costs O ( d ). Therefore, we expect a speedup of O ( d ) by using analog hardware where the capacitance matrix is diagonal (so the cells are decoupled).

If the covariance matrix is sparse, the most efficient strategy is to convert the sparse matrix into a banded matrix, which has a main diagonal and k nonzero central diagonals (if k = 1, the matrix is tridiagonal). By sorting the n _nz nonzero elements of the sparse matrix, the banded matrix can be obtained in linear runtime using the Cuthill-McKee process.

This leads to the special case that sparse matrices are already banded matrices, whose main diagonal and the k central diagonals are all non-zero (for k = 1 the matrix is tridiagonal). In this case, Cholesky factorization can be done in O ( k ² d ) steps. The Cholesky factorization preserves the band structure of the matrix, and therefore the final matrix-vector product used in Cholesky sampling is O ( k d ). Due to the interconnectivity of HMC analog devices, banded covariance matrices are the most natural ones to implement. Hence, in this case we expect a speedup given by analog HMC of O ( k ² d ).

17.4 Higher order correlations The tensors Speedup for distributions with Tensor

We now consider potential speedups in the context of perturbation distributions with nontrivial higher-order cumulants (order 2 and higher); for this purpose, we focus on d-dimensional distributions expressed in terms of the Gram-Charlier A series introduced in (78). Taking the logarithm of such a distribution, the potential energy function is a polynomial of d variables. Every k-th order cumulant involved in the perturbation expansion introduces a homogeneous polynomial of order k into the potential energy function. The coefficients of this polynomial are expressed in terms of a k-th order tensor; thus, if the highest nontrivial cumulant is of order K, then the potential energy function is a polynomial of leading order K.

The potential speedup due to higher-order cumulants is expected to come from the fact that these cumulants are introduced by higher-order couplings or, equivalently, by higher-order tensor interactions. Computations involving these higher-order tensors can be used to compute the gradient of the log probability for HMC.

In analog hardware, a typical k- th order cumulant is introduced by a k-th order coupling, represented as a k - th order tensor. Therefore, to compute the gradient of the log probability digitally, the number of operators scales as O ⁽ dk ). In contrast, in analog hardware, the dynamics introduced by these interactions will naturally occur in O (1). This gives a potential speedup ^of O ( dk ). This analysis applies to the case where there are genuine many-body couplings, as explained above.

When multibody coupling is introduced by chaining two-body couplings, the operations to digitally simulate these dynamics proceed in O ( d2 ). Again, in the analog case, the physical dynamics can be realized in O (1), leading to ^the expected O ( d2 ⁾ speedup.

17.4.1 Speedup for distributions represented by neural networks

This approach, discussed in Section 10, yields a number of speedups:

The d ^1/4 speedup is possible because there is no need to discretize time as is done in digital HMC, and the number of time steps is adjusted according to the dimensionality.

Gradient computation for dense neural networks can cost up to O ( d3 ⁾ steps (if the number of neurons per layer m is scaled with the input dimension d ), and is performed in an analogous way by finite differences. Here, thanks to the Metropolis-Hasting step, the error introduced by finite differences is not as detrimental as in other applications.

Finally, both the accept/reject step computations and the data upload/download can be performed in an analog fashion, avoiding the large number of operations O ( d ) as discussed in Section 16.2.

Therefore, overall, we expect a speedup of d ^13/4 for analog HMC with probability distributions represented by neural networks.

18 Running an Alternative Monte Carlo Process on Hardware

18.1 Realistic analysis of hardware using damping

The protocols for HMC sampling presented so far exploit trajectories of analog circuit dynamics that are assumed to be conservative; these trajectories are then borrowed from a computational protocol that generates samples from a target probability distribution. We now provide a more realistic analysis of physical circuits, and present alternative processes that can be used when the actual circuit dynamics deviate significantly from the intended Hamiltonian dynamics.

Consider the LC oscillator cell of Fig. 3 used for one-dimensional Gaussian sampling. As a first approximation to the actual dynamics of this circuit, resistors are connected in series with both the capacitor and the inductor. If these resistors have resistances R _c and R _L , the circuit will behave as a series RLC circuit with effective resistance R = R _C + R _L . Also, assume that a voltage signal V ( t ) is applied to the circuit and that only V (0) = 0.

The charge q on the capacitor C will serve as the position coordinate of the system, while the magnetic flux LI through the inductor will be the conjugate momentum coordinate.

The time derivative of the momentum LI ( t ) of the circuit is given by the following differential equation:

Until time t The integral of is equal to the charge q ( t ) stored in the capacitor at this point. Thus, the ODE can be rewritten as:

From this relationship between the integral of (108) and the charge on the capacitor, we can see that the time derivative of q ( t ) is I ( t ), as follows:

By rearranging the coupled dynamic equations (109) and (110), we arrive at:

Labeling the momentum p = LI and using dot notation, the above combined equations are reduced to:

Terminology- Rp represents a linear dissipative force in momentum, which does not exist in conservative systems.

We now elaborate on the characteristics of the voltage source V ( t ). Although this voltage source signal can be a deterministic signal provided by the user, we focus on the stochastic voltage signals applied to the system through baseline noise treatment. As a first approximation to a realistic noise model for the RLC cell dynamics, V ( t ) comes from Johnson thermal noise; this noise is caused by thermal agitation of the macroscopic electron collection through an effective resistor with resistance R maintained at a certain temperature T.

Johnson noise can be appropriately modeled using the one-dimensional Brownian motion W ( t ) and the diffusion constant D , i.e., with a little abuse of notation. When we introduce this into equation (113), we see that the RLC dynamics is modeled by a coupled stochastic differential equation (SDE) as follows:

All that remains is to derive the diffusion constant D; for this, we apply the variational dissipation theorem. According to the variational dissipation theorem, which the circuit system follows, the diffusion constant D of the Johnson noise is related to the coefficient of dissipation R (which is the effective resistance), the inductance L , and the temperature of the resistor. In fact, D is given by:

Here k is the Boltzmann constant. With this relation we finally arrive at a realistic dynamic description of the circuit given as an underdamped Langevin SDE;

Therefore, the circuit dynamics is a stochastic process, not a deterministic process, as assumed above. Moreover, stochastic dynamics converge to equilibrium. Suppose that the circuit is initialized at time t = 0 with an arbitrary position q ₀ and an arbitrary momentum p ₀ ; by converging to equilibrium, this initialization This means that it will be distributed according to the stationary distribution π( q , p ) of SDE(118). This distribution can be solved as follows:

When marginalizing with respect to the momentum coordinate p , the equilibrium distribution of the position coordinate (charge on the capacitor) is given by the Boltzmann distribution:

The Boltzmann distribution of an RLC system is a univariate Gaussian distribution with zero mean and kTC variance; if the dynamical trajectory is allowed to converge to equilibrium, repeated measurements of the position coordinates of the circuit (charge on the capacitor) are equivalent to sampling from this univariate distribution.

We now present a protocol to extract samples from an arbitrary multivariate Gaussian distribution using the stochastic circuit dynamics used above.

18.2 Univariate Gaussian Sampling Protocol

Recall that an RLC system converges to equilibrium regardless of the initial values of the positions and momentums. The equilibrium of the system is characterized by a stationary distribution π ( q , p ) in phase space. We try to prepare the RLC circuit in this equilibrium state, which is equivalent to placing the circuit in contact with a hot water bath for a sufficiently long time. When the circuit is in equilibrium, repeated measurements of the position q and momentum p (the charge on the capacitor and the flux through the inductor, respectively) of the circuit are distributed according to π ( q , p ). Measuring only the position effectively marginalizes the momentum and samples from the Boltzmann distribution of the system, which is a zero-mean Gaussian distribution with variance kTC. Adjusting the value of the capacitance C changes the variance of the Boltzmann distribution from which the samples are drawn; for any σ² > 0. If we set , the sampled distribution will have variance σ ² . Since a system in equilibrium will remain in equilibrium unless it is affected by an external source, we can extract an infinite number of samples from the Boltzmann distribution by repeated measurements of the positions of the circuit (charges on the capacitors).

The variance of a sampled Boltzmann distribution can be adjusted, but the mean of this distribution remains zero. If set to , the univariate normal distribution is achieved by adding a fixed mean μ to each location sample. can extract samples from; this can be implemented in analog fashion by applying a constant voltage bias of μ to the RLC circuit, but can also be done digitally. This allows us to have a protocol for extracting samples from arbitrary univariate distributions.

18.3 Multivariate Gaussian Sampling Extension

Although the realistic analysis presented in 18.1 was performed for a single RLC cell, a similar analysis is followed for a network of capacitively coupled RLC cells. Assuming that N RLC cells are coupled in pairs via capacitive bridges, the system has an N-dimensional position vector and N-dimensional momentum vector can be specified by . Then, we define matrices C and L as the capacitance and inductance matrices of the lossless coupled oscillator of Section 3, respectively (i.e., similar to what was done in the coupled LC circuit without resistors).

As the coordinates ( q , p ) of a single unit cell follow the underdamped Langevin dynamics in two-dimensional phase space, the coordinate system follows the underdamped Langevin dynamics in a 2N- dimensional phase space. We assume that the dynamics are given by the following coupled stochastic differential equations:

Here, R is a drag matrix that depends on the resistances of the network, and D is a diffusion matrix multiplied by the N-dimensional Brownian motion W ( t ). Just as the friction coefficient is related to the diffusion constant in the case of a single cell, the drag matrix in the coupled case is related to the diffusion matrix, since the variation-dissipation is still applicable.

Fortunately, one does not need to know the explicit form of the drag and diffusion matrices to know that the dynamics of the coupled RLC network converges to equilibrium. Just as the equilibrium distribution (120) of a single cell does not depend on the friction and diffusion constants, the equilibrium distribution of the dynamics (122) does not depend on the drag or diffusion matrices. If each resistor in the RLC network is maintained at temperature T, the equilibrium distribution over the phase space is defined as:

If we marginalize in the momentum coordinate, the resulting Boltzmann distribution over the N-dimensional position space is given by:

This distribution has a mean of 0 and a covariance matrix of is a multivariate Gaussian distribution. The position vectors of equilibrium RLC networks are Repeated sampling of the multivariate normal distribution It is like extracting samples from .

Then, the protocol for sampling from an arbitrary multivariate Gaussian distribution can be performed similarly to the protocol of Section 18.2, but instead of adjusting a single capacitor, we adjust the value of the capacitance matrix C . If the capacitance matrix is set to C = (1/ kT )∑ , then the sampled Boltzmann distribution will have zero mean and covariance ∑ . As with the univariate protocol, by applying a constant voltage bias μ _i to each of the i cells throughout the dynamics, we can obtain the mean of the sampled distribution as a vector Change to . Extract a sample without voltage bias and vectorize it to that sample. Adding them digitally would have the same effect. The probabilistic dynamics of such a combined RLC network can be used to sample from any multivariate Gaussian distribution.

18.4 Probabilistic Gradient Hamiltonian Monte Carlo

We present a realistic analysis of the single-cell oscillator cell used in the HMC protocol and show that the circuit dynamics undergo underdamped Langevin dynamics, in contrast to the Hamiltonian dynamics. We present a process that exploits the underdamped Langevin dynamics, digitally simulated. A similar sampling protocol can be performed using the underdamped circuit dynamics presented in Section 18.1.

HMC involves directly computing or accessing the gradient of the potential energy function to digitally simulate the Hamiltonian dynamics of the system. For analog hardware, access to the exact gradient may also be limited, depending on whether the user has a digital or analog neural network representing the probability distribution p ( x ).

Stochastic gradient HMC (SGHMC) addresses these challenges and improves the scalability of HMC. SGHMC introduces the use of stochastic gradient computation, which can be performed with less than the full data set. However, the resulting dynamics are no longer those induced by the desired distribution. Therefore, an additional friction term is included to counteract the effects of noise injected by the stochastic gradient.

Explicitly in SGHMC the slope is the entire data set It is calculated using mini-batches, are data points is uniformly sampled and used to compute:

Assuming that the observations are independent, we apply the central limit theorem

We can reach the following expression:

Here, V is the covariance of the stochastic gradient noise from the stochastic gradient approximation. As mentioned above, we can include a friction term to prevent the appearance of this term. This leads to the following dynamic equations defining the SGHMC including friction:

Here, M is the mass matrix, is the diffusion matrix, is the power and is a Wiener process. The influence of stochastic noise is reduced due to the inclusion of friction terms. This is the target probability distribution Obtained over a long period of time corresponding to leads to a normal distribution.

18.4.1 At SGHMC The Korean Approach

Consider performing Gaussian process regression, which learns the posterior probability distribution from mini-batches of training data rather than the entire training set. If we map the log of this mini-batch posterior probability distribution to the potential energy function of the actual oscillator circuit, we obtain a noisy version of the actual potential energy function. Since the potential energy function of the oscillator circuit is a quadratic function defined by the inverse capacitance matrix C ^-1 , the noise will appear as a misspecification of the capacitance values that determine this matrix. When the training set size is large, the central limit theorem suggests that the misspecification of any single capacitance value will be approximately Gaussianly distributed around the actual value.

To attenuate the effects of this modeling error noise and converge to the actual target distribution, the friction (resistance) of the physical oscillator circuit cannot be adjusted to match the noise profile, as is done in digital SGHMC. The friction can be adjusted, but the fluctuation-dissipation relation states that adjusting the friction will proportionally adjust the thermal noise of the system; this is not necessarily true for digital dynamics, but may be inevitable for circuit dynamics. Therefore, the friction term of the physical oscillator should attenuate the effects of both thermal noise and modeling error noise, but can only attenuate the contributions from thermal noise.

To avoid this problem, one can work in a high friction/resistance regime. When this happens, the thermal fluctuations will be much larger than the fluctuations induced by the model-mismatch noise. In this way, the noise profile describing both the thermal noise and the model-mismatch noise approximately satisfies the fluctuation dissipation with respect to the friction term, since the thermal noise is predominantly thermal. When this relation is approximately satisfied, equilibrium is reached. Since the noise part of the potential has been absorbed by the noise term, the Boltzmann distribution is the negative exponent of the true nonminibatch potential, i.e., the true posterior probability distribution. Then, in the high-resistance limit, the true multivariate posterior probability is sampled using the protocol defined in Section 18.3.

19 Langevin Hardware for Monte Carlo

In a variation of HMC, the number of frog-jump steps is set to L = 1. This is known as Langevin Monte Carlo and is defined by

Here is a fixed time step It is a discretization of the following stochastic differential equation:

Here W is the standard Brownian motion. In the limit as t→∞, x approaches a normal distribution, which is exactly the same as the distribution defined by p ( x ). This is the distribution A very simple way to extract samples from: we can use the gradient. In fact, this idea motivated a component of diffusion models known as Noise Conditional Score Networks, because these networks can perform generative modeling simply by learning the gradient of the underlying density function of the data. In fact, Langevin dynamics are increasingly popular in applications involving high-dimensional vector spaces (e.g., parameters of deep neural networks).

19.1 LMC Unit cell of hardware

A unit cell (i.e., a building block) of the LMC hardware is illustrated in Figure 51.

This is different from the unit cell used in HMC, because it can operate without inductors. This is a useful feature of LMC hardware, as inductors are difficult to implement on chip (e.g., they take up a significant amount of space on the chip).

Meanwhile, the unit cell of LMC has a stochastic noise source. This corresponds to a voltage source whose voltage output is Gaussian white noise. For this purpose, there are various ways to implement the stochastic noise source, such as the following possible methods:

1. Analog thermal noise from resistors (also known as Johnson-Nyquist noise).

2. Analog shot noise from diodes such as Zener diodes.

3. Digital pseudo-random noise generated from field-programmable gate arrays (FPGAs) using analog filtering.

Each of these possible noise sources is shown in Fig. 51 as δv A voltage source can be generated. Both analog thermal noise and analog shot noise can be amplified using amplifiers to increase the magnitude of the associated voltage fluctuations.

Each unit cell may also include a voltage source (although this voltage source is not explicitly shown in Fig. 51). This may be useful to shift the mean vector μ associated with the probability distribution p ( x ) away from the value μ = 0 .

19.2 LMC Combining unit cells for hardware

We consider two alternatives for joining the unit cells. One approach is to use a resistive bridge, as shown in Fig. 52a. Another approach is to use a capacitive bridge, as shown in Fig. 52b.

For resistive coupling, the equations of motion for the voltage vector v of the two coupled cells are:

Here

Here, the self-resistance matrix R , the capacitance matrix C , and the conductance matrix J are introduced.

For capacitive coupling, the equation of motion is:

Here the capacitance matrix is as follows:

19.3 Acceptance or rejection of samples

Similar to the HMC case, a Metropolis-Hastings (MH) step can be added to the LMC process. This involves the following steps:

Measure the voltages across each capacitor in the unit cells of the LMC hardware.

This voltage vector v is fed to the digital processor.

Quantifies the energy associated with this voltage vector on a digital processor.

Based on this energy value, the acceptance/rejection step of the MH process is applied to this sample.

This allows the digital processor to be used as a judge of sample quality (via the MH condition) to accept or reject a proposed sample.

19.4 LMC Using hardware Gaussian Distribution sampling

For Gaussian distributions, the slope of the logarithm of the probability distribution is given by an affine function. Specifically, it is given by:

Here, μ is the mean vector and ∑ is the covariance matrix.

Using this formula, equation (131) for the special case of Gaussian can be rewritten as:

This equation can be compared with equations (132) and (133), which are differential equations for LMC hardware with resistive coupling and capacitive coupling, respectively. In fact, equation (135) has a similar form to these equations. Specifically, these equations can be used to obtain the relationship between the distribution parameters and the circuit parameters. That is, the following relationship can be obtained for the case of resistive coupling:

Wow, in the case of capacitive coupling

can be established. These relations allow the user to convert the distribution parameters into circuit parameters of the LMC hardware. For resistive coupling, the following precise relations are provided:

Equations (132) and (133) assume that the mean vector is zero and μ = 0. However, non-zero μ can be allowed by adding voltage sources to each unit cell of the LMC hardware.

19.4.1 Encoding the covariance matrix with resistive coupling

For the special case of resistive coupling, we now provide an analysis of how to encode the covariance matrix as a coupling.

For additional tuning flexibility, a resistor branch is added in parallel with each unit cell, as shown in Figure 53.

Now let's analyze the circuit of Fig. 53. Voltage vector Taking into account, here is the voltage across the i-th capacitor. The following differential equation for v is obtained,

Here

Consider the special case where R = RI and C = CI are proportional to the identity matrix I, and consider the vector w. is defined to be . Then, the dynamics of equation (138) can be rewritten as follows:

Here dw = wdt. Comparing this equation with equation (135), we see that:

So, using the above equation, we can find the covariance matrix can be encoded as a J matrix.

This provides a recipe for how to store the covariance matrix in the parameters of an electrical circuit.

The above analysis generalizes from two unit cells to d unit cells. This essentially entails generalizing the J matrix to the following formula for the matrix elements:

19.5 LMC Using hardware Non-Gaussian Distribution sampling

Beyond Gaussian sampling discussed in Sections 14 and 16, techniques developed for HMC are also applicable to LMC.

19.5.1 Maxwell Demon Approach

The Maxwell demon approach discussed in Section 14 can be applied to LMC hardware to allow sampling from non-Gaussian distributions.

Specifically, the Maxwell Demon can be viewed as a voltage-controlled voltage source acting on the unit cell of the LMC hardware. The same circuit structure as illustrated for the examples in Figs. 36, 38, and 42 can be used, for example, assuming that the HMC unit cell is replaced by an LMC unit cell. (This involves removing the inductor and adding a stochastic noise source to the unit cell, as discussed above.)

19.5.2 Multi-body Combined approach

The multibody coupling approach discussed in Section 16 can also be applied to LMC hardware to allow sampling non-Gaussian distributions.

Specifically, either a voltage-controlled capacitor (VCC) or a voltage-controlled resistor (VCR) can be used as the coupling element between the unit cells. If the HMC unit cell is replaced with an LMC unit cell, the discussion in Section 16 for HMC can be directly applied to the LMC hardware. This includes a discussion on how to design and implement 3-body coupling, 4-body coupling, or n-body coupling using voltage-controlled capacitors. Also, the discussion in Section 16 on using either an analog approach or a digital-to-analog approach for multi-body coupling is relevant to the LMC hardware.

19.6 In the integrator Based on Alternative LMC Hardware

An alternative device to the LMC hardware may be used here, as follows. Specifically, this alternative device employs integrators instead of unit cells, as illustrated in Fig. 54.

The analog device for the LMC of Fig. 54 can be used to sample any probability distribution, including those parameterized by neural networks as described in Section 14.7. To do so, the user must compile a representation of the probabilities into hardware. In the case of neural networks, this would be equivalent to compiling the network into analog hardware and then computing the gradient. Alternatively, if one has a formula for estimating the gradient, one can compile this into analog and the Maxwell daemon can then compute the gradient using the input If fed, it will output the slope of the log probability.

20 Effective temperature increase with noise injection

One question regarding the feasibility of a Gaussian sampling device is whether the attenuation (caused by realistic electrical components) is so large compared to the thermal noise that the measured voltages are very small. The measured values can always be amplified, but adding an analog amplifier would increase the complexity of the device, and there is a concern that the large amplification achieved by digital postprocessing would cause a loss of precision. A physically motivated way to overcome this problem is to simply increase the operating temperature of the device until the thermal fluctuations become so large that amplification is unnecessary; this approach has its own problems, such as the time it takes to reach the desired temperature and the additional energy costs.

Increasing the effective temperature of the device through noise injection can be done by generating white noise using a digital pseudorandom number generator (e.g. generated by an FPGA) and adding the noise to each cell independently. The key question is whether this approach has the same effect as increasing the physical temperature, especially if uncorrelated noise is injected into each cell.

Figure 55 shows a circuit diagram with noise injection into two unit cells. Assuming uncorrelated noise injection, the circuit equations take the following form:

Here And is a Maxwell matrix. The injected noise amplitude is controlled by the parameter κ ₀ , which must have the dimension of charge ² time ⁻¹ for the unit to work, but we leave the dimensionality ambiguous for now.

In fact, the effective temperature can be increased by using uncorrelated noise in each cell. This means that the circuits discussed in the previous sections of this paper can be used to obtain samples without additional amplifiers or temperature control. Assuming uncorrelated white noise injection with the same amplitude in each cell, the resulting covariance matrices for voltage and current are

Here is the Maxwell capacitance matrix, L is the inductance matrix, R is the resistance matrix (which is proportional to the identity matrix), and κ0 is a number that is scaled by the injected noise amplitude and determines the effective temperature.

In this notation, the Maxwell capacitance matrix encodes the inverse of the covariance matrix, not the covariance matrix itself, and is therefore best suited for applications where the precision matrix, rather than the covariance matrix, is specified. However, it is still possible to use this approach in cases where a numerically invertible covariance matrix is given, for example. There are some subtleties in trying to scale the capacitance matrix directly to the covariance matrix, but here we consider the case where the capacitance matrix is scaled by the precision matrix.

21 Thermodynamic systems for solving linear algebraic primitives

Here we show that the same thermodynamic hardware described above can also be used to accelerate key primitives of linear algebra. We exploit the fact that the mathematics of harmonic oscillator systems is affine (i.e., linear), and thus linear algebraic primitives can be mapped onto such systems. We show that a variety of linear algebraic problems can be solved by sampling from the thermal equilibrium distribution of coupled harmonic oscillators. Specifically, we develop thermodynamic algorithms for the following linear algebraic primitives: (i) solving the linear system Ax = b, (ii) estimating the inverse A ^-1 , (iii) solving Lyapunov equations of the form A∑+∑A ^T = I, and (iv) estimating the determinant of a symmetric positive definite matrix A. When implemented on thermodynamic hardware, these methods are shown to scale favorably with problem size compared to digital algorithms.

21.1 Solving linear equation systems

The famous linear system problem is as follows: I'm looking for

Just an invertible matrix and not 0 This is given. Without loss of generality, we can assume that the matrix A in equation (145) is symmetric and positive definite (SPD). To check this, assume that A is not SPD. Then, we define A' = A ^T A and b' = A ^T b and solve the equation as follows.

The transpose of an invertible matrix is invertible, which means that A' is invertible and b' is nonzero, so equation (146) is a nonsingular linear system. The matrix A' is SPD by construction, so if a method is available to solve a symmetric linear system, then x' satisfying equation (146) can be found. Then, multiplying both sides of equation (146) by (A ^T ) ^-1 gives Ax' = b , which means that x' is a solution to the original linear system. This may impact the overall runtime, but still allows for an asymptotic scaling improvement over digital methods (constructing an SPD system from a general system in this way results in the condition number being squared, which impacts performance). Therefore, in the following, we assume A is SPD.

Now we connect this problem to thermodynamics. Consider a macroscopic device with d degrees of freedom, as described by classical physics. The device is assumed to have a potential energy function:

Here . This is a second-order potential that can be physically realized through a harmonic oscillator system, where the coupling between the oscillators is determined by the matrix A, and the b vector represents a constant force acting on each individual oscillator.

We assume that the device is allowed to reach thermal equilibrium with the environment at reciprocal temperature β = 1/k _B T . In thermal equilibrium, the Boltzmann distribution describes the probability that the oscillators will have given spatial coordinates as follows: . Since V(x) is quadratic, p(x) corresponds to a multivariate Gaussian distribution. Therefore, in thermal equilibrium, the spatial coordinate x is a Gaussian random variable as follows.

The unique minimum of V(x) occurs at Ax - b = 0, which also corresponds to the unique maximum of p(x). For the Gaussian distribution, the maximum of p(x) is at the first moment 〈x〉. Therefore, in thermal equilibrium, the first moment is the solution of the linear system of equations:

From this analysis, a thermodynamic protocol can be constructed to solve the linear system, which is depicted in Figure 56. Namely, the protocol involves realizing the potential of equation (147), waiting until the system reaches equilibrium, and then sampling x to find the mean of the distribution. is estimated. This average can be approximated using the time average, which is defined as follows:

Here t ₀ must be large enough to allow equilibrium, should be large enough for the mean to converge to the desired precision. This eventual convergence of the time average to the mean is a key feature of the ergodic hypothesis, which is often assumed in very general thermodynamic systems. The mean can be approximated as the average of a sequence of samples. However, the integration approach can be conveniently implemented in a completely analog way (e.g., using an integrator circuit), which eliminates the need to transmit data from the device until the end of the protocol.

The overall protocol can be summarized as follows:

1. Given a linear system Ax = b, the potential of the device is as follows:

At this time, time t = 0.

2. Balance tolerance parameters Select and choose the equilibrium time as follows:

Here is computed from the physical properties of the system and is calculated by Eq. (155) or (152). The system is allowed to evolve under dynamics until t = t ₀ . δμ and It guarantees.

3. Select the error tolerance parameter δx and the success probability Pδ, and select the integration time as follows:

Here is calculated from the physical characteristics of the system and is calculated by Equation (155) or (152). The time average is measured using an analog integrator as follows:

This is at least With a probability of satisfies.

To implement the above protocol and The desired values of have to be identified, which entails a more quantitative description of the equilibrium and ergodicity. To obtain this description, a model of the microscopic dynamics of the system can be introduced. Considering that the system under consideration consists of a harmonic oscillator in contact with a bath, it is natural to allow for damping (i.e., energy exchange with the bath) and stochastic thermal noise, which is always accompanied by damping due to the fluctuation-dissipation theorem. The Langevin equations account for these effects, and in particular, two general formulations are considered: the overdamped Langevin (ODL) equations and the underdamped Langevin (UDL) equations.

In the overdamped case, we arrive at the following formulas which can be used in the above protocol:

For the case of low damping, we derive somewhat different formulas for the parameters:

21.1.1 Hardware Implementation

Here we describe an electronic device consisting of d coupled RC cells that can be mapped to an overdamped Langevin process and used to implement the thermodynamic linear system algorithm described above.

One potential approach to this hardware implementation is to use the building blocks illustrated in Fig. 17, which employ RC unit cells and resistive coupling between the cells. The equations of motion for the voltage across the capacitors of each cell, v = (v1, v2, ..., vd), are:

Here, C = diag(C1, C2, ... , Cd), R = diag(R1, R2, ..., Rd), w is uncorrelated Brownian motion, and the elements of J are given by

Therefore, there are two resistors within the cell, allowing degrees of freedom for the diagonal elements of the J matrix independently of the R matrix, and there is one resistor that couples each cell. Since R and C are diagonal matrices, And, we can set γ = 1/RC. Let x = v. By selecting , you will be taken to:

This reduces to the overdamped dynamics of the linear system solver when the noise has a mean of γb and a variance of 2γ/β. Note that to implement this in hardware, the values of the intra- and extra-cell resistors must be computed based on the A matrix, which therefore incurs a cost of O(d2) operations at initialization. For electrical hardware, the effective temperature is given by the Johnson-Nyquist noise. , and Δf is the bandwidth of the system. For a bandwidth of 1 MHz, for T = 300 K, Get .

21.2 Estimating the inverse of a matrix

The results in the previous section rely on estimating the mean of x and do not exploit the variation of x around equilibrium. By using the second moments of the equilibrium distribution, it is possible to do more than solve linear systems. For example, it is possible to find the inverse of a symmetric positive definite matrix A. As mentioned, the normal distribution of x is , which means that the inverse matrix of A can be obtained by finding the covariance matrix of x. This can be achieved in a completely analog manner using a combination of analog multipliers and integrators. By setting b = 0 for this protocol, , and so the static covariance matrix is by definition as follows:

To estimate this, time averaging is performed again after the system reaches equilibrium.

The analog component computes the product x _i (t) x _j (t) for each pair (i, j), requiring d2 analog multiplier components. Each of these products is fed into an analog integrator component, which computes one element of the time-averaged covariance matrix as follows.

Although the equilibrium time is the same as the linear system protocol, the integration time is different, because the covariance matrix generally converges more slowly than the mean. Now, assuming ODL dynamics, we give a detailed description of the backestimation protocol as follows.

1. Given a positive definite matrix A, the potential of the device at time t = 0

Set to .

2. Balance tolerance parameters Select and find the equilibrium time.

Allow the system to evolve under dynamics until t = t ₀ , and It guarantees.

3. Select the error tolerance parameter δ _∑ and the success probability P _δ and obtain the integration time as follows:

Using analog multipliers and integrators, we measure time averages as follows:

This is at least with probability P _δ . satisfies.

21.3 Lyapunov Solving equations

This section concerns a device having a controllable noise source such that the covariance matrix of the noise term can be chosen to be an arbitrarily symmetric positive definite matrix. We do not include a linear b ^T x term in the potential, and therefore obtain the following overdamped Langevin equations:

Here R is symmetric and positive definite. In this case, the normal distribution has mean zero and covariance matrix ∑ _S , which is a solution to the following Lyapunov equation.

1. Given two symmetric positive definite matrices A and R, set the potential of the device as

The noise term of the overdamped Langevin equation at time t = 0 is set to . That is, the system evolves under the dynamics of equation (167).

2. Balance tolerance parameters Select and find the equilibrium time as follows.

Allow the system to evolve under the dynamics of equation (167) until t = t ₀ , and is guaranteed.

3. Select the error tolerance parameter δ _∑ and the success probability P _δ and obtain the integration time as follows:

Using analog multipliers and integrators, we measure time averages as follows:

This is with a probability of at least P _δ satisfies.

21.4 Estimating the determinant of a matrix

The determinant of the covariance matrix appears in the normalization factor of the multivariate normal distribution, and the density function of this distribution is as follows.

Hardware that can generate Gaussian distributions can estimate the determinant of a matrix, because this problem is equivalent to estimating free energy differences, an application of probabilistic thermodynamics. The free energy difference between the equilibrium states of potentials V ₁ and V ₂ is:

Assume that the potentials are second-order potentials with V ₁ (x) = x ^T A ₁ x and V ₂ (x) = x ^T A ₂ x . Then each integral simplifies to the inverse of the Gaussian normalization factor,

so

This suggests that the determinant of matrix A ₁ can be found by comparing the free energies and potentials V ₁ and V ₂ of the equilibrium states (where A ₂ has a known determinant), and then calculating

Fortunately, assuming that we can measure the work done on the system as the potential V(x) changes from V ₁ to V ₂ , the free energy difference ΔF can be found. According to Jarzynski's equation, the free energy difference between the initial and final potential (equilibrium) states is

Here is at time t = 0 and time In between, represents the weighted average of the probabilities of all possible trajectories of the system. This can be approximated by

Table 1: Comparison of asymptotic complexity of linear algebra algorithms, where d is the matrix dimension, k is the condition number, and is the error. For thermodynamic algorithms (TA), the complexity depends on the dynamic regime, i.e., whether the dynamics are overdamped or underdamped. For the digital SOTA case, the complexities of symmetric positive definite linear systems, inverse matrices, Lyapunov equations, and matrix-determinant problems are those for conjugate gradient methods, fast matrix multiplication/inversion, Bartels-Stewart algorithm, and Cholesky decomposition-based algorithms, respectively. represents the matrix multiplication constant.

The average over N repeated attempts is as follows:

In summary, the determinant of A ₁ is approximated by

In practice, we might be interested in the log determinant to avoid computational overflow. This would be:

To estimate the log determinant within the error δ _determinant with a probability of at least P _δ , the waiting time is:

21.5 Runtime analyze

The table summarizes various results obtained by comparing with state-of-the-art digital methods for dense symmetric positive definite matrices. These results are based on bounds obtained on correlation time, burn-in time, and physical latency when approaches using integrators and multipliers are employed to ensure convergence of the present method.

22 Conclusion

While various embodiments of the present invention have been described and illustrated in the present disclosure, those skilled in the art will readily envision various other means and/or structures for performing the functions and/or obtaining the results and/or one or more of the advantages described herein, and each of such variations and/or modifications are considered to be within the scope of the embodiments of the present invention described herein. More generally, those skilled in the art will readily appreciate that all parameters, dimensions, materials, and configurations described in the present disclosure are exemplary and that the actual parameters, dimensions, materials, and/or configurations may vary depending upon the particular application for which the teachings of the present invention are used. Those skilled in the art will recognize, or be able to ascertain using no more than routine experimentation, many equivalents to the specific embodiments of the invention described herein. It is therefore to be understood that the foregoing embodiments are presented by way of example only and that, within the scope of the appended claims and their equivalents, embodiments of the present invention may be practiced otherwise than as specifically described and claimed. Inventive embodiments of the present disclosure relate to each individual feature, system, article, material, kit, and/or method described herein. Additionally, any combination of two or more of these features, systems, articles, materials, kits, and/or methods is included within the inventive scope of the present disclosure, so long as these features, systems, articles, materials, kits, and/or methods are not mutually inconsistent.

Additionally, the various inventive concepts may be implemented as one or more of the methods provided herein. The acts performed as part of the methods may be sequenced in any suitable manner. Thus, embodiments may be constructed in which the acts are performed in a different order than illustrated, including performing some of the acts concurrently, even if they are depicted as sequential acts in the exemplary embodiments.

All definitions defined and used in this disclosure should be understood to control over dictionary definitions, definitions in documents incorporated by reference, and/or the ordinary meaning of the defined term.

The expressions corresponding to the use of the indefinite articles “a” and “an” as used in this specification and claims should be understood to mean “at least one,” unless otherwise expressly indicated.

The phrase "and/or" as used in this specification and claims should be understood to mean "either or both" of the elements so linked together, i.e., elements that are conjunctively present in some cases and disjunctively present in other cases. Multiple elements listed with "and/or" should be construed in the same manner, i.e., "one or more" of the elements so linked together. In addition to the elements specifically identified by the "and/or" clause, other elements may optionally be present, related or unrelated to the elements specifically identified. Thus, as a non-limiting example, reference to "A and/or B," when used in conjunction with open-ended language such as "comprising," may in one embodiment refer to A alone (optionally including elements other than B); in another embodiment refer to B alone (optionally including elements other than A); in another embodiment refer to both A and B (optionally including other elements); and so forth.

As used in this specification and claims, “or” should be understood to have the same meaning as “and/or” as defined above. For example, when separating items in a list, “or” or “and/or” should be interpreted to include more than one of a plurality of elements or a list of elements, or optionally, additional items not listed. Only terms that clearly indicate the opposite, such as “only one of,” “exactly one of,” or, when used in the claims, “consisting of,” mean the inclusion of exactly one element of a plurality of elements or a list of elements. In general, the term “or” as used herein is interpreted to indicate an exclusive alternative (i.e., “either one, but not both”) only when preceded by an exclusive term such as “either,” “one of,” “only one of,” or “exactly one of.” “Consisting essentially of,” when used in the claims, has its ordinary meaning in the patent law art.

The phrase "at least one," as used in this specification and claims, when referring to a list of one or more elements, should be understood to mean at least one element selected from one or more of the elements in the list of elements, but does not necessarily include at least one of every element specifically listed in the list of elements, and does not exclude any combinations of elements in the list of elements. Furthermore, this definition allows that, in addition to the elements specifically identified in the list of elements to which the phrase "at least one" refers, there may optionally be elements, related or unrelated, to the elements specifically identified. Thus, as a non-limiting example, "at least one of A and B" (or, equivalently, "at least one of A or B," or, equivalently, "at least one of A and/or B") means, in one embodiment, at least one of A alone, optionally including one or more of A without B being present (and, optionally, including elements other than B); in another embodiment, at least one of B alone, optionally including one or more of A without A being present (and, optionally, including elements other than A); In another embodiment, at least one comprising optionally one or more As, and optionally one or more Bs (and optionally other components); etc.

In the claims, all transitional phrases such as "comprising," "including," "carrying," "having," "containing," "involving," "holding," "composed of," and the like are to be understood as open-ended, meaning including but not limited to. Only the transitional phrases "consisting of" and "consisting essentially of" are closed or semi-closed transitional phrases, respectively, as specified in United States Patent Office Manual of Patent Examining Procedures §2111.03.

Claims (32)

In a system that samples a multivariate distribution,
An analog electrical circuit network for sampling the multivariate distribution, wherein each analog electrical circuit of the analog electrical circuit network includes at least one adjustable passive electrical component; and
A system comprising a digital controller operably coupled to said analog electrical circuit network, said digital controller adjusting said adjustable passive electrical component in accordance with parameters of said multivariate distribution and sampling voltages at points within said analog electrical circuit network. A system in accordance with claim 1, wherein the at least one adjustable passive electrical component comprises at least one of an adjustable capacitor or an adjustable resistor. A system in accordance with claim 1, wherein the at least one adjustable passive electrical component comprises an adjustable inductor. In the first paragraph,
A system further comprising a voltage-controlled voltage source operably coupled to said analog electrical circuit network, said voltage-controlled voltage source modifying said multivariate distribution. In the fourth paragraph, the voltage-controlled voltage source is a system that uses an artificial neural network to relate an input voltage to the voltage-controlled voltage source and an output voltage of the voltage-controlled voltage source. In the fifth paragraph, the input voltage is based on voltages sampled by the digital controller, and the voltage control voltage source is configured to apply the output voltage to the analog electric circuit network to change the potential energy distribution of the entire analog electric circuit network. In the first aspect, the multivariate distribution has a non-zero cumulant of at least the third order, and the system,
A system further comprising a multiport element operably coupled to at least three analog electrical circuits of said analog electrical circuit network, said multiport element generating a k -field term in the potential energy distribution of said analog electrical circuit network as a whole, wherein k is an integer greater than or equal to 3. A system in accordance with claim 7, wherein the multiport element comprises at least one of a transistor or a voltage-controlled capacitor. A system in accordance with claim 1, wherein each analog electrical circuit further includes a probabilistic noise source. A system in claim 9, wherein the probabilistic noise source is one of a thermal noise source, a shot noise source, or a digital pseudo-random noise source. A system in accordance with claim 1, wherein the analog electrical circuits within the analog electrical circuit network are coupled to each other via capacitive and/or resistive coupling. In the first aspect, the system is configured to generate samples from the multivariate distribution via a Langevin Monte Carlo algorithm and/or a Metropolis-adjusted Langevin algorithm. A system in accordance with claim 1, wherein the multivariate distribution is a Gaussian distribution in which a covariance matrix is encoded in at least some of the parameters of the analog electrical circuit network. In a system that samples a multivariate distribution,
A coupled analog unit cell network characterized by a total energy function mapped to the multivariate distribution, wherein each coupled analog unit cell of the coupled analog unit cell network comprises a digitally adjustable capacitor and a digitally adjustable voltage source and/or a digitally adjustable current source; and
A digital controller operably coupled to the above combined analog unit cell network, said digital controller comprising:
Initializing the digitally adjustable capacitors, the digitally adjustable voltage sources and/or the digitally adjustable current sources;
Iteratively reading voltages and currents at nodes of the coupled analog unit cell network and updating the digitally adjustable voltage sources and/or the digitally adjustable current sources, wherein the voltages and the currents correspond to positions and momentums of the total energy function; and
A system that converts the above voltages and the above currents into samples of the above multivariate distribution. A system in accordance with claim 14, wherein the total energy function is Hamiltonian, and the coupled analog unit cell network has dynamics that are used to execute a Hamiltonian Monte Carlo protocol. In claim 14, the coupled analog unit cell network has dynamics used to execute the Langevin Monte Carlo protocol, the system. A system in accordance with claim 14, wherein the coupled analog unit cell network comprises one analog unit cell for each dimension of the multivariate distribution. A system in accordance with claim 14, wherein each analog unit cell of the coupled analog unit cell network comprises a digitally adjustable inductor. A system in accordance with claim 14, wherein the analog unit cells of the coupled analog unit cell network are capacitively coupled to each other. A system in accordance with claim 14, wherein the analog unit cells of the coupled analog unit cell network are resistively coupled to each other. A system in which the analog unit cells of the analog unit cell network in the 14th paragraph are inductively coupled to each other. A system in claim 14, wherein the multivariate distribution is a multivariate Gaussian distribution. A system in claim 14, wherein the multivariate distribution has a non-zero cumulant higher than the second order. In a system that samples a multivariate distribution,
slope calculator;
a momentum integrator device operably coupled to said slope calculator; and
A system comprising a position integrator device operably coupled to said momentum integrator device. In claim 24, the system wherein the slope calculator includes an analog neural network. In claim 24, the slope calculator is a system comprising a network of resistive circuit elements. In a thermodynamic system that samples a target probability distribution,
An analog dynamical system configured to evolve through the Langevin equations and/or Hamiltonian equations of motion; and
A thermo-dynamic system comprising a digital controller operably coupled to said analog dynamical system and configured to receive, at each time step of said Langevin equations and/or Hamiltonian equations of motion, a proposed sample of a target probability distribution from said analog dynamical system and to accept or reject said proposed sample. In the method of finding the inverse of a matrix,
A step of uploading elements of a matrix to each value of the adjustable circuit elements of an analog unit cell network;
a step of allowing the above analog unit cell network to reach thermal equilibrium; and
A method comprising the step of calculating a covariance matrix based on dynamic variables of the analog unit cell network, wherein the covariance matrix is an inverse matrix of the matrix. In the 28th paragraph, the step of calculating the covariance matrix comprises:
A step of extracting samples from the thermal equilibrium distribution of the dynamic variables of the above analog unit cell network; and
A method comprising the step of calculating the covariance matrix of the above samples. A method according to claim 28, wherein the step of calculating the covariance matrix comprises the step of integrating the dynamic variables over time using a plurality of analog integrators. In a method for solving a system of linear equations expressed by matrices and vectors,
A step of uploading elements of a matrix to each value of the adjustable circuit elements of an analog unit cell network;
A step of uploading the elements of the vector to the respective current sources or voltage sources of the analog unit cells;
a step of allowing the above analog unit cell network to reach thermal equilibrium; and
A method comprising the step of calculating average values of dynamic variables of the analog unit cell network, wherein the average values represent solutions to simultaneous linear equations. A method according to claim 31, wherein the step of calculating the average values comprises the step of integrating the dynamic variables over time using a plurality of analog integrators.