A thermodynamic computing system for sampling high-dimensional probability distributions

Abstract

The thermodynamic computing platform samples multivariate distributions (e.g., multivariate Gaussian distributions) using a network of coupled analog unit cells controlled by a digital processor. The analog unit cells are capacitively, resistively, and / or inductively coupled to each other. Each analog unit cell contains an adjustable 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 network of coupled analog unit cells is then evolved through a series of time steps until thermal equilibrium is reached. As 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

[Technical Field]

[0001] (CROSS-REFERENCE TO RELATED APPLICATIONS) This application claims the benefit of priority under 35 U.S.C. §119(e) to 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, each of which is incorporated herein by reference in its entirety for all purposes. [Background technology]

[0002] Traditional artificial intelligence (AI) and machine learning (ML) work best in the absence of uncertainty. However, in reality, uncertainty often exists. Noisy data adds uncertainty to AI and ML predictions, causing them to become unreliable. This impairs the AI's ability to determine whether to defer to human input. This is particularly evident in high-stakes applications of AI, such as criminal justice, healthcare, energy, finance, defense, manufacturing, and robotics. For example, consider a melanoma diagnosis assistance model. Ideally, such a model should defer to the physician when in doubt. Similarly, the AI ​​in a self-driving vehicle should hand over control to the human driver when it is uncertain whether a pedestrian is on the road. However, if uncertainty impairs the AI's predictive capabilities, the AI ​​may not effectively defer to human judgment.

[0003] Uncertainty can be divided into two main categories. First, there is ambiguity that typically arises from noisy data or missing features ("allelic ambiguity"). Alelic ambiguity remains regardless of exposing an AI or ML system to more data. For example, the perception systems of autonomous vehicles include noisy sensors and image processing. Second, there is uncertainty that results from finite training data ("epistemic uncertainty"). Model inputs may look different from the training data or may have few supporting examples in the training data. Epistemic uncertainty decreases as the model is trained on more data.

[0004] Traditional AI deals with uncertainty by assigning a confidence level to its judgments (e.g., 90% confidence that an image shows a particular person or object). By assigning a confidence level to a judgment, the AI ​​quantifies the uncertainty associated with the judgment. Unfortunately, traditional AI and ML generally lack a meaningful formalism for assigning properly quantified uncertainty.

[0005] Figure 1 illustrates how uncertainty affects mainstream deep neural networks' ability to assign meaningful confidence levels to their predictions under real-world conditions. It shows idealized and realistic plots of uncertainty (confidence) in deep neural network predictions, with darker shading indicating higher confidence. Deep neural networks were trained on training data (dots) distributed along the semicircle shown in both plots of Figure 1. The trained deep neural networks were then presented with "out-of-distribution" input data (dots) clustered from x ≈ 0 to x ≈ 1.5 at y ≈ 1.5 in both plots. Ideally, deep neural networks should recognize that their predictions about the input data are highly uncertain or low-confident (light shading), as shown in the left plot. However, in reality, even when the input data is far from the training data, deep neural networks still assign high confidence (dark shading) to their predictions, as shown in the right plot. Unfortunately, common AI and ML models tend to be particularly overconfident in their predictions when operating on unfamiliar input data.

[0006] Fortunately, probabilistic AI and ML can quantify and account for uncertainty in predictions. Probabilistic AI and ML combine AI and ML with probabilistic inference, a mathematical process that determines uncertainty as a function of prior beliefs, biases, and data. For any level of stakes, probabilistic inference should be optimal as long as uncertainty exists. Probabilistic AI and probabilistic inference naturally involve continuous distributions (e.g., Gaussian distributions) to reflect the state of knowledge. Unfortunately, when using standard digital computers, performing probabilistic inference with these continuous distributions tends to be computationally intractable as a function of model complexity. Summary of the Invention [Problem to be solved by the invention]

[0007] Nature performs probabilistic reasoning through the physical principle of thermodynamic equilibrium. Probabilistic AI models are trained and estimated by probabilistic reasoning, typically operating on continuous random variables. Thermodynamic computing platforms can utilize thermodynamic equilibrium to perform probabilistic reasoning for probabilistic artificial intelligence and machine learning. The thermodynamic equilibrium of a system describes the dynamics that evolve a physical system until the system can be described by a stationary probability distribution over its possible states. The technology of the present invention uses such thermodynamic equilibrium in computations with dynamics designed so that the final probability distribution over states matches the target probability distribution being computed. Furthermore, utilizing inherently continuous physical systems naturally fits into the mathematical framework of continuous distributions used in probabilistic reasoning. [Means for solving the problem]

[0008] The techniques of the present invention can be implemented as a thermodynamic system for sampling a multivariate distribution (e.g., a multivariate Gaussian probability distribution). The system may include a network of coupled harmonic oscillators operatively 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 in the network of coupled harmonic oscillators includes a digitally adjustable inductor, a digitally adjustable capacitor, a digitally adjustable voltage source, and a digitally adjustable current source. The digital controller initializes the digitally adjustable inductor, the digitally adjustable capacitor, the digitally adjustable voltage source, and the digitally adjustable current source. The digital controller repeatedly reads voltages and currents at nodes in the network of coupled harmonic oscillators and updates the digitally adjustable voltage source and the digitally adjustable current source based on the voltages and currents. These voltages and currents correspond to the position and momentum of the Hamiltonian. The digital controller also converts the voltages and currents into samples of the multivariate distribution.

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

[0010] Another implementation of the present technology is a thermodynamic system for sampling a target probability distribution. The thermodynamic system includes both an analog dynamic system and a digital controller. The analog dynamic system is configured to evolve via Hamilton's equations over a series of time steps to provide samples of the target probability distribution. The digital controller is operatively coupled to the analog dynamic system and configured to receive, at each time step, a proposed sample of the target probability distribution and to accept or reject the proposed sample.

[0011] Yet another implementation of the present technology is a thermodynamic system for sampling a multivariate distribution (e.g., a multivariate Gaussian probability distribution). The system may 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 tunable passive electrical component (e.g., a digitally tunable inductor, a digitally tunable capacitor, and / or a digitally tunable resistor). In operation, the digital controller adjusts the tunable passive electrical components according to the multivariate distribution and samples voltages at points in the network of analog electrical circuits.

[0012] Such a thermodynamic system may also include a voltage-controlled voltage source operatively coupled to the network of analog electrical circuits to modify the multivariate distribution. The voltage-controlled voltage source may use an artificial neural network to relate an input voltage to the voltage-controlled voltage source to an output voltage of the voltage-controlled voltage source. In this case, the input voltage may be based on a voltage sampled by the digital controller, and the voltage-controlled voltage source may be configured to apply an output voltage to the network of analog electrical circuits to change the distribution of potential energy across the network of analog electrical circuits.

[0013] The multivariate distribution can have at least a third order non-zero cumulant, in which case the system includes a multiport device operably coupled to at least three of the analog electrical circuits in the network of analog electrical circuits. The multiport device can generate k-body terms in the potential energy distribution over the network of analog electrical circuits, where k is an integer greater than or equal to 3. For example, the multiport device can include a transistor, a voltage-controlled capacitor, or both.

[0014] Each analog electrical circuit in a thermodynamic system may also include a stochastic noise source (e.g., a thermal noise source, a shot noise source, or a digital pseudorandom noise source), and the analog electrical circuits in the network of analog electrical circuits may be coupled to each other via capacitive and / or resistive coupling.

[0015] The thermodynamic system can generate samples from multivariate distributions via the Langevin Monte Carlo algorithm and / or the Metropolis adjusted Langevin algorithm.

[0016] The covariance matrix of the multivariate distribution can be encoded in at least some of the parameters of the network of analog electrical circuits.

[0017] Yet another system of the present invention for sampling a multivariate distribution (e.g., a multivariate Gaussian distribution or a distribution with non-zero cumulants higher than second order) includes a network of coupled analog unit cells operatively coupled to a digital controller. The network of coupled analog unit cells is characterized by a total energy function that maps to the multivariate distribution. Each coupled analog unit cell includes 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. During 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 in the network of coupled analog unit cells and updates the digitally adjustable voltage source and / or the digitally adjustable current source. The voltages and currents correspond to the position and momentum of the total energy function, and the digital controller converts the voltages and currents into samples of the multivariate distribution.

[0018] The total energy function can be Hamiltonian, in which case the network of coupled analog unit cells has dynamics used to implement a Hamiltonian Monte Carlo protocol. Alternatively, the network of coupled analog unit cells can have dynamics used to implement a Langevin Monte Carlo protocol.

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

[0020] Yet another system of the present invention for sampling a multivariate distribution includes 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.

[0021] Yet another inventive thermodynamic system for sampling a target probability distribution includes an analog dynamical system operably coupled to a digital controller, the analog dynamical system configured to evolve through the Langevin equations and / or Hamilton's equations of motion, and the digital controller configured to receive, at each time step of the Langevin equations and / or Hamilton's equations of motion, a proposed sample of the target probability distribution from the analog dynamical system and to accept or reject the proposed sample.

[0022] The method of the present invention includes a method for inverting a matrix. The method includes uploading elements of the matrix to respective values ​​of adjustable circuit elements of a network of analog unit cells, allowing the network of analog unit cells to reach thermal equilibrium, and calculating a covariance matrix based on dynamic variables of the network of analog unit cells. The covariance matrix is ​​the inverse of the matrix. Calculating the covariance matrix may include taking samples from the thermal equilibrium distribution of the dynamic variables of the network of analog unit cells and calculating a covariance matrix of the samples. Calculating the covariance matrix may also include integrating the dynamic variables over time with a plurality of analog integrators.

[0023] Another method of the present invention includes a method for solving a system of linear equations represented by a matrix and a vector. The method includes uploading elements of the matrix to respective values ​​of adjustable circuit elements of a network of analog unit cells, uploading elements of the vector to respective current or voltage sources within the analog unit cells, allowing the network of analog unit cells to reach thermal equilibrium, and calculating average values ​​of dynamic variables of the network of analog unit cells. The average values ​​represent solutions to the system of linear equations. Calculating the average values ​​can include integrating the dynamic variables over time with a plurality of analog integrators.

[0024] All combinations of the foregoing concepts and additional concepts described in more detail below (provided such concepts are not mutually inconsistent) are considered to be part of the inventive subject matter disclosed herein. In particular, all combinations of subject matter recited in the claims appearing at the end of this disclosure are considered to be part of the inventive subject matter disclosed herein. Terminology explicitly used herein, and which may also appear in any disclosure made a part of this specification by reference, should be given the meaning most consistent with the specific concepts disclosed herein.

[0025] As will be appreciated by those skilled in the art, the drawings are primarily for illustrative purposes and are not intended to limit the scope of the inventive subject matter described herein. The drawings are not necessarily to scale. In some cases, various aspects of the inventive subject matter disclosed herein may be shown exaggerated or enlarged in the drawings to facilitate an understanding of different features. In the drawings, like reference characters generally refer to like features (e.g., functionally similar and / or structurally similar components). [Brief explanation of the drawings]

[0026] [Figure 1] FIG. 1 shows uncertainty plots for a deep neural network in ideal and real cases. [Figure 2]An overview of a thermodynamic computing system is shown: a digital controller extracts samples from a time-evolving analog system (also known as an analog dynamical system) and implements a Hamiltonian Monte Carlo (HMC) method to sample a target probability distribution. [Figure 3] FIG. 1 is a circuit diagram for a basic LC oscillator that is mapped to a one-dimensional normal (i.e., Gaussian) probability distribution. [Figure 4] FIG. 10 is a circuit diagram for a pair of capacitively coupled oscillators in a network of capacitively coupled harmonic oscillators that are mapped to a multivariate normal probability distribution. [Figure 5] FIG. 10 shows a U-turn displayed while performing HMC on a 2D Gaussian distribution. [Figure 6] System-level diagram of the UL hardware implementation. The PC downloads setup information to the on-board controller and consequently reads sampled data from the cells. [Figure 7] This diagram shows two coupled LC resonators forming a simple two-cell system. The resonant frequency of each resonator can be adjusted independently, allowing for control of capacitive coupling. The cells are excited by a Gaussian noise current source. [Figure 8] FIG. 1 illustrates the general all-to-all connections in the UL hardware setup. [Figure 9a] FIG. 1 shows a tunable cell whose resonant frequency can be digitally controlled by a current noise source realized using a digital controller followed by an LPF. [Figure 9b] FIG. 10 illustrates a tunable cell with an initial condition achieved by adding an additional switch that separates the inductor from the capacitor. [Figure 10] Instead of initializing the capacitor and inductor in each cell, the capacitor can be initialized during disconnection. At a given time after connecting it to the inductor, the required initial condition is achieved. In this way, one component is initialized. [Figure 11] Figure 1 shows a coupling unit between two cells. The cells are connected at nodes A and B. Using control signals, coupling can be enabled or disabled and the polarity is also determined. [Figure 12] Figure 1 shows the FPGA implementation of the controller. (a) Diagram of the high-level control where the download and upload operations take place, as well as the setup of each cell and coupling unit in the FPGA. (b) The cell interface reads the voltage from the ADC and then stores the data. It also controls the switches of the tunable resonators to control the noise amplitude. (c) The coupling unit interface toggles the appropriate switches to control the coupling and its polarity. (d) Diagram of the format of the data downloaded to the FPGA. [Figure 13A] The current noise source of Figure 7 can be implemented by a thin-bit PN generator with a series resistor to emulate the current source. The RC circuit also functions as an LPF. [Figure 13B] FIG. 1 illustrates an example of a linear feedback shift register (LFSR) used in a single-bit pseudo-noise (PN) generator. [Figure 14] 1 shows a visualization of the entire process: the shaded boxes are steps that occur on the digital device, while the white boxes are steps that occur on the analog device. [Figure 15] FIG. 1 illustrates coupling between HMC unit cells via mutual inductance. [Figure 16] FIG. 1 illustrates direct inductive coupling between HMC unit cells. [Figure 17A] FIG. 11 is a circuit diagram of a possible physical realization of a single unit cell consisting of a noisy resistor and capacitor, showing direct inductive coupling between HMC unit cells. [Figure 17B] FIG. 10 is a circuit diagram of a possible means of coupling two unit cells using a coupling resistor. [Figure 18]FIG. 1 illustrates an example of a switching capacitor bank controlled by digital logic, where each physical capacitor adds an additional bit of precision to the overall capacitance. [Figure 19] Two varicaps in opposite directions 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] FIG. 1 is a circuit diagram of a digital-analog system for realizing bipolar coupling. [Figure 21] 1A and 1B show a gyrator-based active inductor (top) and an equivalent circuit with an inductor (bottom). [Figure 22] FIG. 1 is a schematic diagram of a basic symmetric RC unit cell used for overdamped Langevin sampling of one-dimensional normal (i.e., Gaussian) probability distributions. [Figure 23] FIG. 1 is a circuit diagram of a basic coupling cell used to achieve negative coupling between two RC unit cells without an inductor. [Figure 24] FIG. 1 is a circuit diagram of a basic coupling cell used to achieve positive coupling between two RC unit cells without an inductor. [Figure 25] FIG. 1 illustrates various strategies for achieving local connectivity, including square lattices, hexagonal lattices, and King's graphs, where boxes represent unit cells. [Figure 26] FIG. 1 illustrates a further strategy for achieving local connectivity based on a cluster graph including square and hexagonal clusters, where the boxes represent unit cells. [Figure 27] Figure 1 shows a semi-local approach to connecting unit cells (represented by 5-boxes), where unit cells within a row are fully connected, but connections between rows are looser. For D=9, the number of nearest neighbors is either n=3 or n=4, depending on the row. For D=16, the number of nearest neighbors is either n=4 or n=5, depending on the row. In general, n is a number that is used in this connection scheme.

number

[0027] 1. System Overview Figure 2 shows a thermodynamic computing system 200 (also referred to as a thermodynamic computing platform, thermodynamic system, or thermodynamic platform) capable of sampling multivariate Gaussian distributions using Hamiltonian Monte Carlo (HMC), 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 the Gibbs distribution of f to simulating particle trajectories according to the laws of Hamiltonian dynamics, where f acts as a potential energy function. Samples generated by a thermodynamic computing system can be used in probabilistic AI, for example, to weight connections between neurons in different layers of a probabilistic neural network.

[0028] The thermodynamic computing system 200 of FIG. 2 is a hybrid digital-analog system that includes a digital system 210 that functions in conjunction with an analog system 220 (also referred to as an analog dynamic system or a time-evolving analog system). The digital system 210 can be implemented in appropriately 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 being estimated by the thermodynamic system, and a controller 214 for functioning in conjunction with the analog system. The controller can be implemented in firmware. During operation, the controller 214 initializes control parameters 211 for the analog dynamic system 220 and accepts and rejects samples 221 of a target multivariate distribution from the analog system 220. Additionally, the digital system 210 interfaces with a user, either directly or via another digital computing device (e.g., a personal computer (PC)).

[0029] Analog system 220 includes a dynamic physical system 222. Dynamic physical system 222 progresses through Hamilton's equations over a series of time steps to provide samples of a target multivariate distribution (e.g., a multivariate Gaussian, Bernoulli, Binomial, Poisson, or other exponential distribution). More generally, analog system 220 can sample any distribution that can be specified in terms of an energy function that can be realized in (analog) electronic circuits. The target multivariate distribution can be a high-dimensional distribution (e.g., a distribution with 10, 100, 1000, 10,000, or more dimensions) representing probability, heat, population, or any other suitable quantity.

[0030] At each time step in the sequence of time steps, the analog system 220 proposes a sample 221 of the target probability distribution to the controller 214, which either accepts the proposed sample or rejects it if it is an outlier. The natural time evolution of the analog system 220 is mapped to Hamiltonian dynamics, which means that the analog system naturally performs the integration steps of the HMC sampling method. In other words, the natural physical dynamics of the analog system replaces the numerical integration traditionally performed on digital hardware. As a result, a coupled digital-analog thermodynamic system solves these calculations through time evolution.

[0031] A controller 214 in the digital system 210 sets analog system input parameters 211 that represent a multivariate distribution (e.g., a normal (Gaussian) distribution or other distribution that the analog dynamic system 220 is sampling). The input parameters 211 may include the inductance, capacitance, voltage, and current of adjustable inductors, capacitors, voltage sources, and current sources 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 progress for a predetermined time (also referred to as a time step, integration time, or sampling time). After the time has elapsed, the controller 214 reads the analog system's output parameters (which may include voltages, currents, and / or other analog values ​​that represent the final state of the analog system) via the serial bus. These output parameters represent samples 221 of the multivariate Gaussian probability distribution. The controller 214 rejects erroneous or out-of-range parameter readings and corrects analog imperfections in the analog system 220.

[0032] The analog component 220 of the thermodynamic system 200 simulates physical (Hamiltonian) dynamics to sample a multivariate Gaussian or other target probability distribution. The Hamiltonian H of the physical system 222 specifies the total energy (including both kinetic and potential energy) of the physical system. In the context of sampling, the Hamiltonian represents the target probability distribution in terms of q sample space variables in D dimensions. HMC methods traditionally involve taking derivatives of the Hamiltonian with respect to the sample space variables and with respect to the generalized momentum p, and then integrating the trajectory with respect to time t over a series of time steps of duration ε. In this hybrid digital-analog approach, the dynamics of the physical system replaces both the derivative taking and the trajectory integration.

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

[0034] 2. Sampling from a one-dimensional normal probability distribution A single lumped circuit element (e.g., 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 the 1D normal probability distribution in a thermodynamic computing platform that performs HMC calculations.

[0035] 2.1 Hamiltonian Monte Carlo (HMC) on 1D normal probability distribution Mean μ and variance σ 2 Consider a normal distribution with , so that the probability distribution can be written as

number

[0036] The energy function corresponding to this system is:

number

[0037] Let p be the canonical momentum corresponding to x, and choose some kinetic energy function independent of x. A standard choice is

number

number

[0038] After specifying the Hamiltonian, HMC relies on two main components to sample from the target distribution. Hamiltonian dynamics Given a Hamiltonian H and initial conditions x0 and p0 for x and p, HMC integrates Hamilton's equations over time.

number

number

[0039] The structure of the HMC benefits from the property that x and p are independent random variables, which follows from the fact that U(x) and K(p) have non-overlapping supports.

number

[0040] Connecting the two components is the choice of initial conditions for the integration. The initial values ​​x0 and p0 set up a surface of constant probability density along which the Hamiltonian dynamics propagate. Up to normalization, the probability density along the trajectory is:

number

[0041] Because the marginals are independent, resampling of p can be used to choose a new constant value for the joint probability density. Hamiltonian dynamics then exchanges this resampling of momentum for a resampling of position. To quote Neal et al., “MCMC using Hamiltonian dynamics,” Handbook of Markov Chain Monte Carlo, 2.11 (2011), “[x] is unchanged, and p is drawn from its correct conditional distribution given [x] (which is the same as its marginal distribution, due to independence), so this step clearly leaves the canonical joint distribution unchanged.”

[0042] 2.2 Hamiltonian and Electric Circuits In the next subsection, we outline three steps to derive a Hamiltonian from a lumped-element circuit.

[0043] 2.2.1 Coordinate Selection Suppose we start with an LC oscillator like the one in Figure 3. The LC oscillator has an inductor with inductance L and a capacitor with capacitance C. We want to write down the Lagrangian for this circuit so that we can derive the corresponding Hamiltonian.

[0044] The method begins by selecting coordinates for all branches and nodes in the circuit, where branch coordinates I L and I C We chose as the currents in the inductors and capacitors, and marked one node as the reference (ground) and the other node with the coordinate V (the voltage difference from that node to ground).

[0045] 2.2.2 Derivation of the Lagrangian We then write the total energy of the circuit in terms of these chosen coordinates. For a circuit, the total energy is the sum of the energies of the branches. From Kirchhoff's law, I L =-IC From the capacitor definition, I C = C(dV / dt). Therefore, the energy of a capacitor or inductor can be expressed as the relationship between voltage V and its time derivative.

number

[0046] Since voltage is treated as a position coordinate, capacitor energy is potential energy U and inductor energy is kinetic energy K.

[0047] We can now write down the Lagrangian for the system.

number

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

number

number

[0049] 2.2.3 Legendre transformation to obtain the Hamiltonian Since we are attempting to perform Hamiltonian Monte Carlo, the Hamiltonian of the physical system must match the Hamiltonian corresponding to the probability distribution. To that end, we can now derive the Hamiltonian for the LC oscillator:

[0050] The Hamiltonian function for a 1D system is defined as follows:

number

number

[0051] The momentum conjugate p is p=LCI in terms of the current through the capacitor. C It can be written as:

[0052] 2.3 Connecting a Gaussian distribution to an LC oscillator Recall the two Hamiltonians we obtained: one is an abstract Hamiltonian corresponding to a Gaussian distribution of the virtual mass m.

number

number

[0053] First, note that capacitance C is equal to the inverse variance. The mean value can be implemented in hardware using a constant voltage offset, or in software by simply adding the mean value to an otherwise zero-mean sample. The mass is LC 2 Therefore, after choosing the mean and variance, the inductance L becomes a free parameter that allows us to choose the mass.

[0054] This is an example of how to map a probability distribution onto an equivalent electrical circuit.

[0055] 3 Capacitively Coupled Harmonic Oscillators and Multivariate Normal Probability Distribution 3.1 Distribution Specification An N-dimensional multivariate normal distribution has an N-dimensional mean vector

number

[0056] Here,

number

number

[0057] Taking the logarithm gives us the "energy function" corresponding to this distribution.

number

number

[0058] 3.2 Capacitively Coupled Oscillators 3.2.1 Circuit Lagrangian Consider a set of N harmonic oscillators, each capacitively coupled to all the others. Any pair of harmonic oscillators, i, j, can be represented by the subcircuit shown in Figure 4.

[0059] The total inductive energy in the circuit is:

number

[0060] If we define the recursive matrix as follows:

number

number

[0061] The total capacitance energy is:

number

[0062] To simplify the subsequent calculations, we define Maxwell's capacitance matrix C as follows:

number

[0063] Then, the capacitive energy can be rewritten as:

number

[0064] Suppose we treat the current through the inductor as a position coordinate. To write down the Lagrangian, we rewrite the voltages in terms of the time derivatives of these position coordinates. For each voltage,

number

number

[0065] Now we can write the Lagrangian for the system as follows:

number

[0066] 3.2.2 Circuit Hamiltonian The Hamiltonian is given by N positions qi , N momentum p i , and the Lagrangian L relationship is defined as follows:

number

number

[0067] In this example, the momentum is:

number

number

number

[0068] The Lagrangian can then be rewritten as:

number

[0069] This means that the Hamiltonian of the system is

number

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

number

number

number

[0071] And the Hamiltonian can be calculated as follows:

number

[0072] 3.3 Changing coordinates The existing code base for HMC and its extensions assumes that the target distribution is encoded by potential energy. Therefore, if the distribution could be similarly encoded in potential energy, the engineering workload could be significantly reduced. This involves changing the coordinate system used to describe the analog electrical system. The original coordinate system

number

number

[0073] 3.3.1 General theory First, consider whether the original Hamiltonian system satisfies the following:

number

number

[0074]

number

number

[0075] 3.3.2 Proof of new coordinates Building on the results from the previous section, we have a procedure for checking that a transformation is canonical.

number

number

[0076] Consider the original Hamiltonian for the current system.

number

[0077] Assume new coordinates.

number

number

[0078]

number

number

number

[0079] 3.3.3 Final Physical and Abstract Relationships The provisions are summarized below: We have an abstract Hamiltonian.

number

number

[0080] 3.3.4 Solving inductor currents in terms of capacitor currents Consider again the circuit of Figure 4. The current at node i is:

number

[0081] The current through the capacitor is

number

number

[0082] Define the code matrix.

number

[0083] Therefore, we can combine the sums by writing:

number

[0084] Going into matrix form, we get:

number

[0085] 3.4 Connecting to the multivariate normal distribution Setting the values ​​in the circuit so that C=Σ, the potential energy of the electrical system in the transformed coordinate system is equal to the energy of the target multivariate normal.

[0086] Recall that Maxwell's capacitance matrix is ​​defined as follows:

number

[0087] This means that for discrete capacitors, the following relationships hold equivalently:

number

[0088] Given a covariance matrix Σ, we invert the relationship and set the individual capacitances by first setting the coupling capacitors as follows:

number

[0089] The capacitance to ground can then be set as follows:

number

[0090] Consider an instance of an HMC process running on a 2D Gaussian distribution shown in Figure 5. This figure illustrates the importance of carefully choosing the integration time: integrating for too long means returning to places already visited, resulting in wasted computation.

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

[0092] To do this, we assume we have a system of coupled oscillators with a mass matrix M and a spring constant matrix K. In terms of a physical mass-spring system,

number

[0093] Then the oscillator position vector

number

number

[0094] By solving the equation of motion, the unknown eigenvectors

number

number

[0095] This can be rearranged as follows:

number

[0096] A=K -1 M and λ=ω -2 and we obtain the standard eigenvalue equation:

number

[0097] Then, for each eigenvalue λ, f=ω / (2π)=λ -1 / 2 The standard frequency is obtained by setting / (2π).

[0098] Given these standard frequencies, we know how long to integrate before starting a U-turn in the slowest moving (lowest frequency) direction: f min Let be the minimum normal mode frequency, and after one period we are back where we started.

number

number

[0099] 3.6 Frequency Scaling The oscillator frequency is

number

[0100] 3.6.1 General theory average

number

number

[0101] A bijective and continuously differentiable transformation φ is defined as follows.

number

[0102] Let us define a new function f as follows:

number

[0103] If X is some random variable with density f, then φ(X) has density g. Let us see if it is convenient to sample from density f.

[0104] The following applies.

number

number

[0105] The synthesis is as follows:

number

[0106]

number

number

[0107]

number

number

number

[0108] 3.6.2 Circuit parameter perspective Thus we obtain the scaling procedure. We consider that increasing the frequency of a physical oscillator is done by decreasing either the inductance or the capacitance, or both. For a capacitively coupled oscillator, Σ physical = C, so for a fixed inductance L, Σ physical Setting <Σ increases the frequency.

[0109] In other words, if η>1, the frequency is high.

number

number

number

[0110] 3.6.3 Quantified frequency scale changes How much does the frequency change? Recall that the eigenvalue equation for a system of coupled oscillators is:

number

[0111] In terms of the system covariance matrix, this becomes:

number

[0112] Substituting the target covariance Σ and the scaling factor η, we obtain:

number

[0113] If ω represents the normal mode frequency of the original target distribution, then

number

number

[0114] Therefore, the factor η used in the definition of φ is exactly the amount by which the normal mode frequencies are multiplied.

[0115] 4 Physical Implementation FIG. 6 shows a hardware implementation that primarily realizes the sampling portion of the overall computational process. The hardware implementation includes a digital system, here implemented as an FPGA controller 610, operably coupled to an analog system, here shown as an analog-to-digital converter (ADC), a digital-to-analog converter (DAC), a noise source, and an LC resonator 620 coupled to a switch 630. The sample is the voltage level across the LC resonator 620 (also referred to as a unit cell or cells, and described in more detail in Section 4.1) when excited by the noise source 630. A user can interact with the hardware implementation via a personal computer (PC) 602, which is operably coupled to the FPGA controller 610.

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

[0117] The interface between the FPGA controller 610 and the PC 602 shown in Figure 6 is a virtual COM port over USB. Other interface protocols are possible (e.g., PCIe, or CXL).

[0118] 4.1 Coupled resonator Figure 7 shows a coupled resonator in the form of two capacitively coupled parallel LC circuits. The coupling and cell capacitances are determined by a controller. The capacitors can be adjusted in an analog or digital manner, while the inductance can remain constant.

[0119] Coupling more resonators improves system performance. The nature of the coupling (e.g., all-to-all, linear, etc.) is determined by the desired processing and hardware limitations. Figure 8 shows a conceptual all-to-all coupling between N cells.

[0120] The exact implementation of the cells and coupling sections is technology dependent. Sections 4.1.1 and 4.1.2 provide examples of how to implement cells in technologies compatible with printed circuit boards (PCBs).

[0121] 4.1.1 Tunable resonators Figure 9a(a) shows the structure for a switch-tunable resonator. The capacitor's value can be changed using a shunt switch (e.g., a transistor). For the smallest capacitance value (C1), no switch is required. The top plate of the capacitor represents the voltage that is sampled and coupled to other cells.

[0122] The noise source can be implemented by a digital controller to reduce or minimize the reliance on analog circuitry, and a low-pass filter (LPF) then filters out unwanted and / or undesired high frequency components in the output signal, as will be discussed in more detail in Section 4.4.

[0123] 4.1.2 Controlled Binding The coupling unit between cells can be implemented as shown in Figure 11. A series switch (controlled by a coupling enable signal) allows the coupling to be turned on or off. The polarity of the coupling is controlled by the coupling polarity of the signal. Negative and positive polarity is achieved by using a transformer setup as shown in the figure.

[0124] 4.2 Analog-to-Digital Converter (ADC) Voltage readings from each cell are taken using an ADC. This configuration uses a single ADC channel per cell. Alternatively, if multiple cells are sampled simultaneously, they can share a single ADC. The input impedance of the ADC needs to be considered with respect to the resonant frequency of the cell.

[0125] 4.3 Controller (FPGA) The controller within the PCB handles the following functions: · Download and upload interface with PC. ·Interface with ADC (reading and saving data, initialization, etc.). · Control the cell switch to adjust the resonator. Controlling coupling switches with the appropriate polarity to turn them ON or OFF. Generate pseudorandom noise (Section 4.4).

[0126] Figure 12 shows a block diagram of the FPGA's functions, along with the protocol for the setup data from the PC.

[0127] 4.4 Noise injection Each unit cell has an independent noise source. The noise source in Figure 7 is a current source because the cells are set up as a parallel LC circuit. 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 Figure 13A. The presence of the resistor makes the PRN voltage source look like a current source (Thevenin to Norton conversion).

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

[0129] 4.4.1 PN polynomial Single-bit pseudorandom sequences are implemented by a linear feedback shift register (LFSR) configuration that represents a primitive polynomial. Because each cell starts with a different initial condition, the cells generate uncorrelated pseudo-noise sequences. The length of the shift register must be large enough to ensure that the noise in each unit cell is sufficiently separated from other sequences. A typical shift LFSR is shown in Figure 13B.

[0130] 5. How it works 5.0.1 System initialization At this stage, the device is prepared to sample from a particular multivariate normal distribution according to the process shown in FIG. 14 (starting with system initialization 1410). 1. Define a probabilistic model 1411 (digital device): target covariance matrix Σ and target mean vector

number

number

[0131] 5.0.2 Extracting Samples In the sampling stage 1420, the digital device (digital controller) and the physical device (network of coupled analog unit cells) work together to advance the Markov chain for voltage. The iterative portion of the following procedure is an adaptation of the Hamiltonian Monte Carlo (HMC) process. It uses momentum instead of position as the chain state, and the physical device integrates via its natural dynamics rather than the digital device via numerical integration. The accept / reject step of the loop is a new 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. Select V1421 (digital device). The latest state in the Markov chain S is the next voltage value.

number

number

number

number

number

number

number

number

number

number

number

[0132] 5.0.3 Termination Sequence The end sequence stage 1430 converts back from physical quantities to target distribution samples. 1. Exit 1431 (digital device). Once N samples have been drawn, stop adding values ​​to the Markov chain. 2. Calculate momentum 1432 (digital device). Change the variable from an electrical quantity (voltage) to an algorithmic quantity (momentum). This is done by

number

number

number

number

[0133] 6 Inductive coupling 6.1 Inductively Coupled Harmonic Oscillators for Sampling Multivariate Distributions 6.1.1 Circuit Lagrangian A set of N harmonic oscillators, each with a mutual inductance L ij Any pair of oscillators i, j can be represented by the following subcircuit:

[0134] The total inductive energy in the circuit is:

number

number

[0135] To simplify the subsequent calculations, the inductance matrix L is defined as follows:

number

[0136] The inductive energy can then be rewritten as:

number

[0137] Suppose we want to treat the voltages between nodes in a circuit as position coordinates. To write down the Lagrangian, we rewrite the currents in terms of the time derivatives of these position coordinates. For each current,

number

[0138] Therefore, the inductive energy can be written as:

number

[0139] Now, the Lagrangian for the system can be written as follows:

number

[0140] 6.1.2 Circuit Hamiltonian The Hamiltonian is given by N positions q i , N momentum p i , and in terms of the Lagrangian L, is defined as follows:

number

number

[0141] In this example, the momentum is:

number

number

number

[0142] The Lagrangian can then be rewritten as:

number

[0143] This means that the Hamiltonian of the system is

number

[0144] 6.1.3 Direct inductive coupling Consider a set of N harmonic oscillators, each directly inductively coupled to all the others. Any pair of oscillators, i, j, can be represented by the subcircuit shown in Figure 16.

[0145] The total inductive energy in the circuit is:

number

number

[0146] This means that the inductive energy does not couple into adjacent cells.

[0147] 6.2 Combining inductive and capacitive coupling More generally, a multivariate normal probability distribution, whether its coupling is inductive, capacitive, or inductive and capacitive, is mapped onto a network of coupled harmonic oscillators, again with 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.

[0148] In this general case of inductive 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 according to the derivation presented above. Therefore, for simplicity, we omit the derivation and state that inductive and capacitive coupling can be combined.

[0149] 7 Resistive coupling An alternative means of combining unit cells is to use resistors. In this case, we can consider a unit cell containing a noisy resistor at a 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 resistor's terminal capacitance C is also added to the equivalent resistor model. The voltage at node 1 (here simply labeled v(t)) is a state variable whose dynamics follows:

number

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

[0151] When building a system of many unit cells, one may wish to couple them together to express correlations and geometric constraints. As an example, two unit cells can be coupled through a resistor, as shown in Figure 17(B). The coupled unit cells (represented by the voltages at nodes 1 and 2) are coupled through their drift terms as follows:

number

number

number

[0152] Here we introduce the self-resistance matrix R, the capacitance matrix C, and the conductance matrix J. We can then construct a system of many coupled unit cells using the basic building blocks shown in Figure 17.

[0153] Adjustable configuration for 8 components To be able to represent any Gaussian distribution, hardware components must be changed or tuned for each distribution, and the choice of components must be made with precision vs. range trade-offs in mind.

[0154] 8.1 Adjustable configuration for capacitors Two ways that tunable capacitors can be implemented in HMC hardware are (1) using switching capacitor banks and (2) using variable capacitor (varactor) diodes.

[0155] 8.1.1 Switching Capacitor Bank Figure 18 shows a switching capacitor bank. A switching capacitor bank is made up of multiple capacitors that can be connected in parallel in any combination to achieve many different capacitance values. Digital signals are used to actuate the switches to connect or disconnect the capacitors. Thus, the device effectively behaves like a tunable capacitor, and the capacitance value is adjusted by a digital voltage signal, as shown in Figure 18.

[0156] Because the device can include any number of capacitors of any capacitance value, the device can theoretically have any range and accuracy at the expense of a large physical area, losses, and parasitic capacitors. Furthermore, because the device is constructed from regular capacitors, the device has good overall linearity and operates normally over a wide range of voltage swings.

[0157] 8.1.2 Varactor diode Varicap diodes are semiconductor devices with voltage-dependent capacitance values. In this type of diode, the depletion layer is particularly sensitive to voltage bias. Therefore, the capacitance value can be changed by changing the voltage bias across the diode. These diodes are specifically designed for a desired tunability range. Varicaps generally have a narrower tunability range than switching banks of capacitors, but the tunability is continuous.

[0158] Another consideration for these semiconductor devices is their nonlinearity, which means that the device behaves as a linear capacitor only for very small fluctuating voltages. At larger AC fields, nonlinear distortion becomes significant. Finally, because these devices are diodes, they are polarized; that is, they behave as a capacitor only in one direction of DC voltage bias.

[0159] Figure 19 shows a possible approach to constructing a tunable capacitor based on varicaps, where two varicaps are oriented in opposite directions and a voltage is applied between the two varicaps by an external source. This approach allows the overall capacitance to be controlled using an external voltage (e.g., a voltage signal from a digital device).

[0160] 8.2. Configuration for Bipolar Coupling Capacitors Any target covariance will have some elements that are negative, representing negative correlations between variables. Bipolar (positive and negative) capacitive coupling allows us to represent these elements in the hardware capacitance matrix.

[0161] 8.2.1 Transformer system As previously discussed with respect to Figure 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, and the transformer changes its sign. The transformer is wound to effectively cancel the sign of the voltage from one side to the other, thus achieving an effective negative capacitive coupling. In addition, a switch is placed between the capacitor and the transformer that can be switched between the wires going to the transformer or adjacent cells to turn the negative coupling on or off.

[0162] 8.2.2 Hybrid digital-analog system An alternative way to implement bipolar coupling is the hybrid digital-analog approach shown in Figure 20. In this approach, the hardware has uncoupled unit cells, each equipped with a means to measure voltage, which are capacitively coupled to an arbitrary voltage source. This voltage source measures the voltage V = (v 1, v 2, ...,v N-1, v N ) depending on the voltage f i (v) can be controlled to supply cell i. This applied voltage can be digitally cancelled, effectively implementing negative coupling.

[0163] 8.3 Adjustable configuration for inductors Real inductors have the disadvantage of taking up a large amount of space on-chip, so having an alternative circuit element that acts like an inductor can be useful to reduce the area used on-chip.

[0164] Active inductors can be used. For example, Figure 21 shows a circuit diagram for a gyrator-based active inductor. Here, two Gm cells and a capacitor provide a transfer function that mimics the behavior of an inductor within a given frequency range. The overall inductance can be adjusted by changing the capacitance of the capacitor.

[0165] This active inductor approach has several limitations, including a limited frequency range, a low quality factor, increased power consumption, and nonlinearity. Nevertheless, this circuit can be useful for reducing the on-chip area for inductors.

[0166] 9 Strategies for Eliminating Inductors 9.1 Introduction As mentioned previously, analog computation primitives for Gaussian sampling (and others) involve encoding variables of interest to the application into degrees of freedom in circuit unit cells. Correlations between variables of interest can be encoded via capacitive coupling of the unit cells. However, as discussed previously in Section 8.2, simple capacitive coupling is not suitable for encoding both negative and positive correlations, or bipolar correlations.

[0167] One way around this problem is to use a transformer in series with a capacitor wired to cancel the voltage, if you want a coupling of opposite sign to pure capacitive coupling. However, this solution has the drawback of using two inductors per coupling element. Furthermore, in an all-to-all connection graph of unit cells, the number of coupling elements is d times the number of unit cells, d. 2 This is a significant problem because inductors take up a large physical area when placed on a silicon chip, limiting the number of unit cells per chip.

[0168] To address this problem, we use a unit cell architecture based on RCR circuits. The degree of freedom of interest is the voltage difference across the capacitors. Bipolar coupling can then be achieved 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).

[0169] 9.2 Unit Cell A proposal for a unit cell circuit is shown in FIG.

[0170] Analysis of this circuit using Kirchhoff's current laws yields the following equations describing the current in each branch:

number

[0171] Summing these equations gives:

number

number

[0172] Equation 24 can be expanded in terms of node voltages as follows:

number

number

[0173] To simplify the equation, we add a new set of variables and

number

number

number

[0174] To characterize the asymmetry between the resistances within the unit cell, we define the asymmetry parameter α i is defined as follows:

number

[0175] Equation 30 then becomes:

number

number

[0176] Finally, substituting Equation 34 in Equation 33 gives:

number

[0177] The dynamics of this unit cell has one degree of freedom

number

[0178] 9.3 Merged Cells Bipolar capacitive coupling between unit cells of the type described in Section 9.2 is described below. The circuit contains two switches, which are operated so that there are only two possible configurations: one for positive coupling and one for negative coupling.

[0179] 9.3.1 Negative Coupling In the first case, consider the case where the effective coupling is negative and the capacitive coupling is between like nodes of the unit cell (e.g., node ia to node ja, node ib to node jb). The circuit diagram is shown in Figure 23. Analysis of this circuit using Kirchhoff's current law yields the following equations describing the current in each branch for cell i:

number

number

[0180] Following the analysis from the unit cell, we take the sums and differences of the current equations for cell i,

number

number

[0181] These equations can be expanded in terms of voltage, and for cell i we get

number

number

[0182] The above expansion used the coordinates defined in 28 for each cell, and used the resistance asymmetry parameter from equation 32, and a new capacitance asymmetry parameter defined as follows:

number

[0183] 9.3.2 Symmetric Case In the symmetric case, α i =α j = 0 and β = 0, the following set of equations is obtained for cell i,

number

number

[0184] 9.3.3 Positive Bonds In the second case, consider the case where the effective coupling is positive and the capacitive coupling is between different nodes of the unit cell (e.g., node i a to node j b , node i b to node ja ). The circuit diagram for this case is shown in Figure 24. Following the analysis in Section 9.3.1, the following equations for the sum and difference current equations are obtained in cell i:

number

number

[0185] The equations are very similar to equations 44-47, but the coupling terms between the degrees of freedom of interest are

number

[0186] 9.3.4 Symmetric Case In the symmetric case, α i =α j = 0 and β = 0, we obtain the following set of equations for the cell:

number

[0187] 9.4 Discussion The above derivation shows that this scheme can achieve two degrees of freedom by changing the position of the switches in the coupling cell.

number

number

[0188] 10 Hardware Connectivity Strategy 10.1 Local versus global connectivity As previously mentioned, the inventive device can be used to sample from a Gaussian probability distribution. Conveniently, a Gaussian distribution can be fully characterized by its first two moments (mean and covariance). As explained in Section 3.6, the mean can be incorporated into each oscillator variable as a constant offset in post-processing. However, the covariance matrix must be incorporated into the system dynamics, which is represented by the capacitance matrix in the physical device. Given a target Gaussian probability distribution of d variables, a dxd 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 d self-capacitances on the diagonal and all two-oscillator coupling capacitances on the off-diagonal.

[0189] To build a system of coupled oscillators, we need to determine the degree of connectivity that can be implemented. Connectivity refers to how many connections any oscillator in the system has with other oscillators. In the previous section, we assumed that the system of coupled oscillators is all-to-all connected. This is the highest degree of connectivity that can be achieved for a given coupling scheme, where every oscillator is connected to every other oscillator.

[0190] In this coupling scheme, the hardware provides a coupling of order d between the oscillators. 2 (the dimension of the problem). This can become a space bottleneck in larger systems. This problem can be addressed by reducing the degree of connectivity and having only local coupling between a set number of oscillators. The degree of connectivity can be specified by the maximum number of connections for a given oscillator. For example, a simple connectivity scheme is a degree-2 connectivity architecture, where an oscillator is coupled to only its two nearest neighbors.

[0191] However, for less than full connectivity, we run into the problem of mismatch 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:

number

[0192] When using devices with all-to-all connectivity, the following capacitance matrix is ​​obtained:

number

number

number

[0193] Given limited hardware connectivity, it is necessary to quantify the impact of covariance mismatch on the overall process behavior.

[0194] 10.2 Hardware with Local Connectivity 10.2.1 Overview If a specific application is targeted, the properties of the corresponding covariance matrix can be directly mimicked by hardware. For example, hardware with connectivity k can be designed for a specific class of applications with a banded covariance matrix. Furthermore, the accept / reject step process can overcome some degree of 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 decay exponentially with distance from the diagonal, a process running on banded local hardware should still be able to extract high-quality samples.

[0195] Figures 25 and 26 show some possible local connectivity hardware architectures that can be implemented for the appropriate class of distributions. Implementing these strategies in hardware can provide significant savings in on-chip space and complexity, unlocking the potential for larger systems while efficiently sampling from certain families of covariance matrices.

[0196] 10.2.2 Uniform grid Figure 25 shows graphs of three different examples of uniform lattices: a square lattice, a hexagonal lattice, and a king. 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.

[0197] 10.2.3 Cluster Graph Figure 26 shows an alternative way to achieve local connectivity based on the concept of a cluster graph, where a unit cell is a cluster of some size n c For example, for square clusters, n c = 4, n for hexagonal clusters c= 6. Within a cluster, there is full connectivity, which means that each unit cell is connected to every other unit cell in the cluster. Between different clusters, there is only sparse connectivity, as shown in Figure 26. For square and hexagonal clusters, the number of nearest neighbors is n = 5 and n = 7, respectively. Other types of clusters can also be considered, for example, octagonal or other shapes.

[0198] Cluster graph connectivity has applications in image processing, for example, because the covariance matrices for image processing often have a block-circular structure. Blocks in this block-circular structure may correspond to clusters in the cluster graph connectivity. By associating fully connected blocks in the covariance matrix with fully connected clusters in the hardware connectivity, the hardware connectivity can be matched to the structure of the covariance matrix (e.g., for image processing).

[0199] 10.3 Hardware with Semi-Local or Global Connectivity The local connectivity discussed above has the property that the number of nearest neighbors, n, is independent of the dimension D. Now consider a connectivity structure where n increases with D.

[0200] 10.3.1 Semi-local connectivity At one extreme, global connectivity in hardware is difficult to achieve in practice. At the other extreme, local connectivity may limit the speedup achievable with analog hardware compared to digital approaches. Intermediate connectivity that is neither fully global nor fully local may be an attractive alternative.

[0201] This intermediate connectivity can be said to be semi-local, meaning that the number of nearest neighbors n grows as some monotonic function of D, but the connectivity is not completely global.

[0202] Figure 27 shows the semilocal connectivity between unit cells organized on a square lattice. Each unit cell is fully connected to other unit cells in the same row, and the connections between rows of unit cells are sparse. In this scheme, the total number of connections is

number

number

[0203] This strategy essentially corresponds to cluster graph connectivity (discussed above), where each cluster is a row and the cluster size is

number

number

[0204] 10.3.2 Full Connectivity Finally, consider the extreme case of global connectivity, also known as complete connectivity, where the total number of connections is D 2 and the number of nearest neighbors n increases as D.

[0205] Fully connected hardware can be designed for applications with arbitrarily dense covariance matrices. Figure 28 shows one possible layout of eight fully connected oscillators, following the description in Section 4.

[0206] Fully connected designs are the most expressive and allow for the highest possible speeds on digital hardware, but they are expensive in terms of physical space and on-chip complexity.

[0207] Using the method for realizing bipolar coupling outlined in Section 8.2.2, we can effectively build fully connected hardware with only degree D connections.

[0208] 11 Dealing with hardware device defects 11.1 Thermodynamic error correction One advantage of using thermodynamic hardware to form samples for HMC is the approach's native error resilience. This is due to the presence of a Metropolis-Hastings (MH) step, in which samples are accepted or rejected based on the energy difference between the proposed sample and the previous sample. Incorporating this step, along with the analog hardware that proposes samples, can be thought of as thermodynamic error correction. The accept / reject step serves to filter out any erroneous samples proposed due to hardware limitations (e.g., capacitor inaccuracies or connectivity constraints, as discussed below). Thus, this naturally occurring form of error correction potentially makes the thermodynamic approach to HMC resilient to the presence of imperfections.

[0209] The MH step is typically performed on digital hardware in this methodology, which allows for any errors accumulated in the analog dynamics to be corrected. Therefore, the combined analog-digital interface (during the MH step) corrects for errors.

[0210] In a sense, the MH step serves to keep the analog dynamics on the correct trajectory, meaning that the analog dynamics may deviate away from the desired trajectory due to device imperfections (such as inaccuracies and noise), but the MH step serves to push the position-momentum coordinates back to the desired trajectory.

[0211] One consideration when examining this approach to thermodynamic error correction is the behavior of the acceptance rate (also known as the success rate or acceptance probability), which reflects the quality of the proposed samples. The following subsections provide some deeper investigations into this quantity, highlighting its behavior due to capacitor imprecision, tolerance, and range, and device connectivity constraints.

[0212] 11.2 Dealing with Capacitor Inaccuracies, Tolerances, and Ranges As mentioned earlier, the accept / reject step process acts as a form of error compensation for the hardware. If the distribution described by the physical hardware deviates from the target distribution, samples drawn from this hardware distribution are rejected without compromising the quality of the overall sample. This leads to a trade-off between the overall time to draw the required sample and the precision or accuracy of the hardware. This trade-off can be quantified by calculating the average acceptance rate for the process using imperfect hardware. This error compensation means that the hardware can be constructed using a more conservative level of tunability and still reliably represent a useful Gaussian distribution.

[0213] Devices can then be designed with this concept in mind. The tunability range, along with the number of bits of precision for capacitors and inductors, can be selected by optimizing the acceptance rate and constraints of the hardware architecture. For example, Figure 29 shows the effect of capacitor inaccuracy on acceptance rate, based on numerical simulations. Because acceptance rate is only a weak function of inaccuracy, it is possible to use capacitors represented by a small number of bits in hardware (e.g., 2 or 3 bits) and still maintain good performance.

[0214] An additional consideration during the design phase is component tolerances, which refer to possible deviations of component values ​​from their nominal values. These additional deviations can also be corrected by an accept / reject step (see Figure 30). However, these tolerances can be taken into account by a full characterization of the hardware before it is put into use.

[0215] 11.3 Addressing connectivity limitations Consider the case of limited connectivity (e.g., the local connectivity discussed in Section 10.2). Limited connectivity at the device has several consequences when applying this hardware to HMC. In particular, as we will see later, the expected improvement over digital approaches must be directly related to the device's connectivity for some applications (e.g., generating samples from a multivariate Gaussian distribution).

[0216] As mentioned above, the nature of the accept / reject step again serves to ensure that the accepted samples represent the desired distribution. One metric to consider is the acceptance probability, and we need to be confident that this probability will not decrease significantly.

[0217] Another approach to mitigate the possible negative effects of limited hardware connectivity can potentially be mitigated by classical preprocessing, which is described below for the case of HMC for multivariate Gaussian sampling.

[0218] 11.3.1 Classical Preconditioning for Local Connectivity Constraints To increase the speedup that can be achieved from the hardware even when connectivity is limited, methods can be developed that aim to increase the magnitude of the covariance values ​​within a given bandwidth. A permutation matrix can be used to permute the connected cells, searching for a near-optimal ordering that maximizes the weights for a given bandwidth.

[0219] One way to achieve this permutation is to construct the Laplacian of the weighted graph using the absolute value of the covariance matrix. The second smallest eigenvector of the Laplacian defines the Fiedler vector. The ordering of the elements of this vector contains the approximate optimal order for labeling each of the unit cells.

[0220] Therefore, to obtain the optimal mapping of the covariance matrix values ​​to the hardware, we compute the second smallest eigenvector of the Laplacian constructed from the covariance matrix. The complexity of this operation depends on the properties of the matrix. In the worse case, this can be O(d 2 ) operations. If the covariance matrix is ​​sparse, this should be cheaper to compute.

[0221] Another approach to achieve the same goal is to put the dense matrix into sparse form by direct truncation, and then apply the process cut-hirmakey to put the sparse matrix into band form. Direct truncation can be applied during the construction of the covariance matrix, and if some elements are smaller than some cutoff value, they can be set to 0. This can therefore be completed during the computation of the matrix itself. The cut-hirmakey process then determines how many non-zero elements are in the sparse matrix O(n nz ) depending on how much remains in the original image. This is O(d 2 ) so by applying a simple truncation followed by a cut-hill-make key, we can nz ) step, the dense covariance matrix can be put into band form. The number of non-zero elements depends on the truncation, a parameter that must be set by the user.

[0222] For the specific example of multivariate Gaussian sampling for a random covariance matrix whose elements decay exponentially with distance from the diagonal, Figure 31A shows that the acceptance rate is approximately constant with increasing dimension. When the covariance matrix is ​​a dense random matrix, Figure 31B shows the exponential decay of the acceptance probability with dimension when attempting to generate samples from an approximation of the true covariance matrix with bandwidth k. These two examples highlight the different expected performance for hardware applied to different types of covariance matrices. Furthermore, they also highlight the potential benefit of preprocessing the covariance matrix (discussed above) to move larger elements closer to the diagonal and potentially avoid exponential scaling.

[0223] 12 Thermodynamic error reduction methods Error mitigation can mitigate the effects of errors associated with thermodynamic sampling devices. Specifically, this error mitigation method targets imprecision errors (errors associated with the inaccuracies of analog hardware components). The acronym for this method is THERMIES (Thermodynamic Error Mitigation by Imprecise Ensemble Sampling).

[0224] 12.1 Univariate Protocols We start with an example in the one-dimensional case to illustrate the essential concepts, and then consider the more general case in the appendix (see Appendix: Error Mitigation for Thermodynamic Sampling Hardware). Consider a Gaussian sampling device that can realize a random variable X whose probability density function is

number

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

[0226] In practice, samples may come from distributions with variances that are not multiples of ε. For example, suppose we want to sample a normal distribution N[0,1.5ε] (here and below, N[a,b] denotes a normal distribution with mean a and (co)variance b). This distribution is shown in Figure 32.

[0227] To sample (approximately) this distribution, we follow the steps that form the basis of the Thermy method:

[0228] (Univariate Thermie Protocol) 1. Sample a Bernoulli random variable B∈{0,1} with Pr(B=0)=1 / 2fmo and Pr(B=1)=1 / 2. 2. If the result is B=0, sample the distribution N[0, ε], and if the result is B=1, sample the distribution N[0, 2ε]. 3. Instead of memorizing the results of the Bernoulli trial, record the results as realizations of the random variable X.

[0229] Then the probability density function of the random variable X is:

number

number

[0230] In this example, f a The distribution is plotted as a dashed line in FIG. 32 and matches the target distribution better than either of f1 and f2.

[0231] mixed f a =(1-w)f m +wfm+1 If we form (w∈[0,1]), then f a The variance of is as follows:

number

number

number

number

[0232] In the univariate case, the inaccuracy problem can be avoided by rescaling the random variables to have realizable variances. That is, one can define new random variables Y ∝ X such that Y~N[0,mε] for some m∈{1,2,...}, and therefore error mitigation is not necessary. This can be difficult or impossible in the multivariate case, which is why error mitigation methods are useful. We are currently developing a multivariate generalization of the Tellme method.

[0233] 12.2 Inaccuracy Dependencies It is desirable that the quality of the samples is not affected by the precision of the hardware implementation. Fortunately, the Tellme protocol eliminates the linear dependence of the approximate distribution on ε as ε goes to zero. Figure 33 shows the f with and without Tellme error mitigation. a and f t The L1 distance between ε and ε is shown. When the error is mitigated, the slope disappears at ε=0, but when it is not mitigated, the slope clearly does not disappear at ε=0. This numerical result provides evidence that Thermy (or similar methods) can reduce or eliminate the sensitivity of sample quality to hardware inaccuracies.

[0234] 13 Extension to non-Gaussian distributions 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 interesting. Therefore, we next consider strategies for building analog hardware that can generate samples from non-Gaussian distributions.

[0235] 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). In the case of a Gaussian HMC sampler, a quadratic potential can be specified as a function of the position coordinate x. However, other potentials (other than quadratic potentials) can be realized, and this is the conceptual basis for extending beyond Gaussian sampling.

[0236] The stochastic dynamics of the HMC process guarantees that in the long-time limit, the position and momentum are distributed according to the Boltzmann distribution (given in equation (2)). This is true for any exemplary choice of the potential function U(x). Furthermore, the distributions on x and p are independent, and the marginal distributions on x are

number

[0237] For a quadratic potential, the Boltzmann distribution over the position coordinate is Gaussian, and therefore the sample at position x is sampled from a Gaussian distribution ~exp(-U(x)). As mentioned previously, this quadratic potential is the natural physical potential associated with LC oscillator systems. However, if we can modify the shape of the known quadratic potential energy function U(x) provided by an LC oscillator circuit, we can sample from a distribution associated with a non-quadratic potential.

[0238] Below we present two approaches to modifying the shape of the potential function. An approach based on the thermodynamic concept of Maxwell's demon ·Many-body coupling approach to LC oscillators

[0239] 13.2. Overview of the Maxwell's Demon Approach Section 14 details an approach to sampling from non-Gaussian distributions based on the concept of Maxwell's demon.

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

[0241] Simply put, the voltage source output by Maxwell's demon device artificially creates a new potential energy function U(x) for the existing position coordinate x of the LC oscillator circuit. This allows the LC oscillator system to evolve according to a new Hamiltonian. As a result, by appropriately choosing U(x), the probability distribution being sampled can be programmed.

[0242] 13.3 Overview of the Many-Body Coupling Approach In Section 16, we detail an approach to sampling from non-Gaussian distributions based on designing many-body couplings (e.g., two-body couplings) of unit cells.

[0243] A multivariate Gaussian distribution can be described by its first two moments (i.e., the mean vector and the covariance matrix). Therefore, a distribution with non-trivial higher-order moments (e.g., the third moment) will be non-Gaussian. Therefore, by focusing on generating higher-order moments in a distribution, it is possible to go beyond the Gaussian distribution.

[0244] Within the HMC framework, this corresponds to generating higher-order coupling terms in the Hamiltonian. While two-body couplings are ubiquitous in physical systems, higher-order couplings (e.g., cubic) are rarer. Therefore, these couplings in physical systems need to be generated judiciously.

[0245] Roughly speaking, the reason for realizing two-body coupling in the design of LC oscillator systems is the use of two-port devices (e.g., capacitors) to couple the unit cells, which suggests that three-body coupling can be achieved by replacing the two-port devices with three-port devices.

[0246] In fact, three-port devices can be used to couple three unit cells. For example, a transistor has three ports including a source, a drain, and a gate. This creates a true three-body coupling in the Hamiltonian, and therefore potentially a third moment in the corresponding probability distribution being sampled. Varicaps (also known as varactors) are an alternative to transistors for this purpose.

[0247] More generally, sets of 3-port devices can be chained together to obtain an n-port device (,n ≥ 3). This allows for the generation of true n-body coupling between unit cells and therefore n-th moments in the corresponding probability distributions. This is discussed in detail in Section 16.

[0248] 13.4 Further Distributions of Interest 13.4.1 Univariate Distributions Here we provide some mathematical details about the types of probability distributions (other than Gaussian) that can be sampled.

[0249] First, for simplicity we consider univariate distributions. Some notable examples include: ·Exponential distribution Gamma distribution Laplace distribution ·Generalized normal distribution ·Exponential distribution family

[0250] Consider two distributions that are specified only for non-negative numbers (i.e., the domain x∈[0,∞]). We start with the exponential distribution: p(x)=βexp(-βx). (71)

[0251] Next, we consider the gamma distribution, which is a generalization of the exponential distribution.

number

[0252] where α>0 is the shape parameter, β>0 is the rate parameter, and Γ(.) is said to be the gamma function. The gamma distribution can be turned into an exponential distribution by setting α=1.

[0253] The Laplace distribution, unlike the Gamma distribution, has support all over R. The distribution is given in terms of a "location" parameter μ∈R and a scale parameter σ.

number

[0254] Next, consider the generalized normal distribution with the following probability density:

number

[0255] Choosing β = 2 corresponds to a standard Gaussian distribution, and choosing β = 1 corresponds to a Laplace distribution. Thus, equation (74) includes the Gaussian and Laplace distributions as special cases.

[0256] These distributions are special cases of the exponential family, which can be written as: p(x|θ)=h(x)exp(η(θ)T(x)-A(θ)) (75)

[0257] 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).

[0258] 13.4.2 Multivariate Distributions In general, one 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 a nontrivial third-order cumulant, its probability distribution will be distorted. Higher-order cumulants define other features of the shape of the distribution. Cumulants are described in more detail below.

[0259] The exponential family described above can be extended to the multivariate case as follows: p(x|s)=h(x)exp(s·t(x)-A(s)) (76) where h(x) is a scaling constant, s is a vector called the natural parameter, t(x) is called the sufficient statistic, and A(s) are logarithmic partitioning functions.

[0260] 13.4.3 Cumulant Method for Distributions In Section 2.1 we describe circuits for sampling from normal distributions. To move on to sampling from non-Gaussian distributions, we may consider the problem of sampling from a distribution that is a small perturbation away from the Gaussian distribution. For what small parameters should we take this perturbation expansion?

[0261] The cumulant of the distribution is defined as follows:

number

[0262] It turns out that these are 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, because the normal distribution has kappa functions for n ≥ 3. n = 0. Therefore, the answer for the perturbation expansion parameter is κ ≥ 3.

[0263] The formalism for this expansion is known as the Gram-Charlier A series (closely related to the Edgeworth series), and is as follows:

number

number

number

[0264] To understand what this may mean, we include only the first two amendments.

number

[0265] Now, it's clear that keeping only the first two correction terms won't be valid for all x, since the expression inside the logarithm can be negative. However, we're assuming that κ3 and κ4 are small anyway, and this x is far enough away that this pathology can be ignored. For a fully rigorous expansion, we could consider the Edgeworth series instead. However, we'll just proceed by performing the expansion for log(1+x)≈x.

number

[0266] 14 Maxwell's Demon Approach to Sampling Non-Gaussian Distributions 14.1 Overview of Maxwell's Demon Sampling from a non-Gaussian distribution involves a change in the potential energy function of the system. To this end, we introduce the concept of Maxwell's demon. Figure 34 illustrates the basic idea of ​​Maxwell's demon. Historically, James Clerk Maxwell considered an experiment in which an intelligent observer monitored gas particles in a box with two chambers and selectively opened the door only for the faster particles, thereby separating the faster and slower particles over time. This reduced the entropy of the gas system over time, but did not violate the second law of thermodynamics, since entropy was generated elsewhere. The intelligent observer was called Maxwell.

[0267] As detailed below, Maxwell's Demon allows HMC hardware to go beyond Gaussian distributions. There are various ways to physically construct a Maxwell's Demon (MD), including digital and analog methods, as described below.

[0268] 14.2 Maxwell's Demon in the Context of HMC 14.2.1 Changing the Potential Energy Landscape To perform HMC sampling for different distributions, we need to change the form of the potential energy U(x) function in the circuit Hamiltonian H, as shown in Figure 35.

[0269] In Gaussian sampling, this potential energy function U(x) is a quadratic form given in terms of the covariance matrix Σ, as shown in the left panel of Figure 35. The coupled LC oscillator system (discussed earlier) has this kind of quadratic potential, which can be sampled from a Gaussian distribution P(x).

[0270] Introducing a Maxwell's Demon device within the existing LC oscillator framework modifies the circuit Hamiltonian by changing the potential energy landscape with respect to the circuit variable x. An "ideal" Maxwell's Demon can apply a force to a known HMC cell according to any suitably smooth potential energy function U(x). The dynamic state x of the LC unit cell is passed to the Maxwell's Demon, which effectively evaluates f=-∇U(x) and returns the calculation as a force to the LC cell. This allows us to adjust the probability distribution P(x) being sampled by the HMC. The probability distribution is an exponential function of the potential energy up to normalization. P(x)~exp(-U(x)) (80)

[0271] 14.2.2 MD as a voltage-controlled voltage source Figure 36 shows how a Maxwell's demon device can be viewed as a voltage-controlled voltage source (VCVS). In this case, the diagram shows a single unit cell, but the concept can be generalized to multiple unit cells.

[0272] As shown, the voltage x across the capacitor in the i-th unit cell i is input to a Maxwell's Demon (MD) device, which can be a digital or analog device. The MD device processes this input by applying some function f. The output f(x i ) is sent back as a voltage source and applied as a voltage inside the original unit cell. Thus, the MD device functions as a VCVS, since the voltage output is controlled by the voltage across the capacitor.

[0273] 14.2.3 Modified potential for one unit cell The mathematical analysis of how an MD device modifies the potential energy function begins by considering the case of a single unit cell, as shown in FIG.

[0274] We assume that the MD device does not draw any current from the unit cell. This assumption simplifies the analysis, but can be relaxed with a more complicated derivation.

[0275] Under this assumption, Kirchhoff's voltage law can be applied to the circuit of Figure 36 to obtain the following equation:

number

[0276] where:

number

number

[0277] This equation can be rewritten in terms of Newtonian physics notation: x i =V i Set to the position and p i =I i (m i / C i ) is set to momentum (m i is the mass), and the force is f i (x i ), which gives us the following equation:

number

[0278] The potential energy function can then be obtained by taking the indefinite integral of the force.

number

[0279] From this equation, any desired potential energy function U(x i ) can be obtained because g i (x i ) can be freely chosen as any function. i (x i )=x i -[L i C i / m i h i (x i ) can be selected. i (x i ) is some arbitrary function. Then:

number

[0280] Therefore, the output of Maxwell's demon device g i (x i By appropriately selecting , the potential energy function can be designed to be virtually any desired function.

[0281] 14.2.4 Modified potential for multiple unbonded unit cells Extending the above analysis to multiple unit cells can be done in two steps: first, consider the case where capacitive coupling between unit cells is neglected, and then include this coupling in the next subsection.

[0282] In the case of multiple unit cells, the MD device takes as input the entire state vector x={x 1, x 2, ...,x d}, and then outputs a voltage vector g(x). g can act on the entire state vector x, effectively coupling the unit cells together (even in the absence of direct capacitive coupling).

[0283] In this case, we obtain equations similar to those previously described for the evolution of each unit cell.

number

[0284] The total force can be written as a vector. f(x)=F(-x+g(x)) (86) where:

number

number

number

[0285] Therefore, by appropriately choosing g(x), virtually any desired potential energy function U(x) can be set up.

[0286] For example, g(x)=xF -1 We can choose h(x), which gives us the potential energy as follows:

number

[0287] 14.2.5 Modified potential for multiple coupled unit cells For multiple coupled unit cells, the analysis is more complex. Nevertheless, at a conceptual level, it is not significantly different from the uncoupled case. The derivation obtained above for the uncoupled case provides important intuition for how the coupled case behaves. Therefore, for simplicity, we omit the derivation of how the potential is modified in the coupled case.

[0288] 14.3 Digital Maxwell's Demon Figure 37 shows how a digital Maxwell's Demon device interacts with an analog LC oscillator system. This fits into the paradigm shown in Figure 36 (where the MD device is viewed as a voltage-controlled voltage source).

[0289] The digital MD device can be stored and processed on either a CPU or an FPGA. ADCs and DACs are used to convert between analog and digital signals. Specifically, the voltage across a capacitor in a unit cell (represented by x) can be converted to a digital signal and input to the digital MD device. The digital MD device then outputs a digital signal, which is then converted to an analog voltage vector g(x). The components of this vector are applied as voltage sources in the appropriate unit cells.

[0290] The main advantages of the digital approach to MD devices are their flexibility and programmability. Digital devices offer the flexibility to choose virtually any function g(x) for the output. In this sense, digital approaches to MD devices are programmable and flexible, and can sample from a wide range of probability distributions P(x).

[0291] The main drawback of digital approaches to MD devices is that they do not fully utilize the potential speedup that can be obtained from analog hardware. Digital devices can calculate the gradient of the potential energy function, which can be difficult for complex potential energy functions. In contrast, in analog approaches to MD devices, the gradient of the potential energy function corresponds to a physical force, and this physical force is applied naturally to the LC oscillator system without the need for any computation (this is explained in more detail later). Therefore, there is a greater potential for computational speedup in analog MD devices.

[0292] 14.4 Motivations for Analogue Maxwell's Demon Maxwell's demon (MD) devices, when constructed in an analog fashion, offer the opportunity to significantly accelerate computations. This is because analog MD devices can apply forces to an LC oscillator system without the need for computation to calculate the forces. This is in contrast to digital MD devices, which calculate forces from a potential energy function. This calculation is generally expensive and can have nontrivial scaling (e.g., linear or quadratic scaling) with respect to the dimensionality of the problem. Therefore, avoiding the need to calculate the gradient of the potential energy function is an advantage of analog MD approaches.

[0293] 14.5 Uncorrelated approach to analog Maxwell's demon 14.5.1 Overview Next, we consider approaches to constructing MD devices that do not correlate the main unit cells, before introducing an alternative approach that correlates the main unit cells in Section 14.6 below. Each LC oscillator in the HMC hardware described above is referred to as a main unit cell. In the following, we introduce additional circuits, referred to as auxiliary unit cells.

[0294] The decorrelation strategy for constructing an analog Maxwell's demon involves several components shown in Figure 38 and described as follows: 1. Measuring the voltage on the main unit cell capacitor, which is then applied as a voltage source to drive current through the auxiliary unit cell. 2. Nonlinear circuit elements within the auxiliary unit cell that allow for complex I-V relationships. 3. Measuring the voltage across the resistor in the auxiliary unit cell (to read the current in this cell), which is then inverted in sign and applied as a voltage source in the main unit cell.

[0295] These components within a single unit cell are shown in Figure 38. In practice, the same structure can be copied d times for d unit cells, possibly with different hyperparameters (e.g., different nonlinear elements) for each copy.

[0296] The voltage output by the auxiliary unit cell is reversed in sign before being applied to the main unit cell. This is because the system needs to have negative feedback (rather than positive feedback) to remain stable. Negative feedback provides a restoring force similar to how a spring acts in physics to pull a mass backwards towards an equilibrium point. In summary, voltage reversal helps maintain the stability of the system.

[0297] However, Figure 39 shows that the voltage inverter can be avoided by not grounding the resistor in the auxiliary unit cell, which allows the wires carrying the voltage difference across this resistor to be physically reversed (or swapped), as shown in Figure 39.

[0298] There are several additional elements that can be inserted into the circuit to increase flexibility. These optional additional elements include: An amplifier with adjustable gain inserted in the coupling between the two unit cells. For example, the amplifier can amplify the voltage reading on 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 reading on the resistor of the auxiliary unit cell to boost the voltage source applied to the primary unit cell. An additional voltage source inserted within the main unit cell, which is constant (independent of the state variables) but adjustable. This adds flexibility by providing a constant offset to the voltage output by Maxwell's demon, which in turn causes a constant offset in the force applied to the main unit cell.

[0299] 14.5.2 Auxiliary unit cells As shown in Figures 38 and 39, the basic components of the auxiliary unit cell are a voltage source (coming from the main unit cell), a nonlinear element, and a resistor.

[0300] One strategy is to use the IV characteristic of a nonlinear element (NLE) as a function generator. For example, as mentioned in the discussion of equation (81), the input V i and output g i (V i ) we want to generate a complex functional relationship between the IV characteristics of the NLE. The IV characteristics of the NLE provide one way to design this complex functional relationship. Therefore, the desired function g i (V i ) is transformed into a problem of designing the desired IV characteristics for the NLE.

[0301] Below we describe explicit strategies for building NLEs, based on individual circuit elements or on a combination of multiple elements.

[0302] 14.5.3 Nonlinear elements for unimodal distributions Nonlinear elements (NLEs) can be used to sample from unimodal distributions. We know that NLEs with monotonic IV characteristics mathematically induce unimodal probability distributions. Therefore, we restrict our attention here to NLEs with monotonic IV characteristics.

[0303] For the current through the ith NLE, as a function of the voltage from the ith main unit cell, I i (V i ) and the output of Maxwell's demon is expressed as g i (V i ) for simplicity, we will denote this as the force f acting on the i-th principal unit cell. i (V i ) can be set equal to the equation (82)(x i V i As shown in Figure 1, g i (V i ) is V i This can be justified by assuming that the resulting potential energy U i (V i ) is the force f i (V i) and the associated local probability distribution P i (V i ) minus an exponential function of the potential energy. For simplicity of notation, the current, force, potential energy, and probability distribution can be written simply as I(x), f(x), U(x), and P(x), respectively.

[0304] Figure 40 shows a schematic 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 in a nearly unidirectional manner, with the current increasing exponentially with forward bias voltage. However, there is a breakdown voltage that allows current flow in the reverse direction. This results in a slightly asymmetric IV characteristic, as shown in the first panel of Figure 40. The second panel shows the force f(x), which is simply the negative of I(x) (due to the voltage inverter shown in Figure 38). The third panel shows U(x), and the fourth shows P(x). Because P(x) is asymmetric, it is non-Gaussian and has a steeper drop-off than a Gaussian distribution.

[0305] Instead of considering a single diode, consider two diodes oriented in parallel in opposite directions. Selecting these two diodes in parallel as the NLE results in a symmetric distribution P(x) because the circuit elements are symmetric with respect to positive and negative voltages. This choice of NLE allows us to design a symmetric version of the probability distribution shown in Figure 40.

[0306] Transistors also offer nonlinear I-V characteristics. For example, in the case of a 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 (input characteristic) results in a distribution P(x) similar to that of a diode, while the latter case (output characteristic) results in a distribution P(x) with a long, slowly decaying tail. Thus, the latter allows for non-Gaussian distributions with long tails. Again, these distributions can be made symmetrical around zero voltage by placing two transistors in parallel with opposite orientations. Other configurations (e.g., common emitter) are also possible for transistors.

[0307] 14.5.4 Nonlinear elements for multimodal distributions Next, we describe NLEs for sampling from multimodal distributions. To obtain multiple modes, the potential energy must have multiple minima (e.g., local and global minima). This means that the derivative of the potential energy (the force) must change sign more than once. If the force changes sign multiple times, it must be a non-monotonic function. Therefore, the IV characteristic of the NLE must be a non-monotonic function (because it is the negative of the force). Therefore, we describe NLEs with non-monotonic IV characteristics.

[0308] First, consider a tunnel diode as the NLE. Figure 41 shows the relevant curves for current, force, potential energy, and probability distribution. The current appears somewhat cubic and is clearly nonmonotonic, with local maxima and minima in the first quadrant. Shifting the force by a constant negative amount then results in a force f(x) that switches sign multiple times. This can be 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 main unit cell. The resulting f(x) curve is shown in the second panel of Figure 41. This results in a potential energy U(x) with two minima and a probability distribution P(x) that is multimodal with two modes.

[0309] As an alternative to tunnel diodes, Gunn diodes are also considered as candidates for NLEs. Gunn diodes also have nonmonotonic IV characteristics under forward bias, with local maxima and minima in the first quadrant. In addition, they have nonmonotonic IV characteristics under reverse bias. Therefore, the curve is nonmonotonic for both positive and negative voltages. This can result in a more complex probability distribution P(x) than that of tunnel diodes. In particular, for Gunn diodes, P(x) can have multiple modes in both the positive and negative voltage regions.

[0310] Other NLEs (e.g., memristors) can also be considered. In addition, complex IV characteristics can be generated by combining multiple nonlinear elements in series or parallel. For example, when NLEs are coupled in series, the overall IV curve is a composite of the individual IV curves. In this way, complex IV curves can be designed by combining various circuit elements to form the entire NLE block.

[0311] 14.6 Correlation approach to analogue Maxwell's demon 14.6.1 Overview Here we consider an approach to constructing an analog Maxwell's demon that involves correlating multiple unit cells. This strategy is outlined in Figure 42 and includes the following components: 1. Voltage measurements on each main unit cell capacitor, the collection of which forms a voltage vector v. 2. An analog combining operation in which the elements of a voltage vector v are combined with each other. The output of the combining operation is another vector h(v). This combining can correspond to an affine transformation h(v) = Av + ​​b, or to a higher-order tensor (e.g., a quadratic or cubic form). 3. A nonlinear operation acting on each component of the vector h(v). This involves using the voltage components of h(v) to drive a current in a nonlinear circuit element, and then reading the current through this nonlinear element (NLE) by taking the voltage across a resistor in series with the NLE. See, e.g., Sections 14.5.3 and 14.5.4 for examples of NLEs including various types of diodes and transistors. The resulting voltage vector is written as n(h(v)), where n corresponds to the nonlinear operation. 4. The previous two steps may be repeated multiple times so that there are multiple layers of connections and nonlinear operations. If there are L layers, the final output vector has the form:

number

number

[0312] In the following, we consider different ways to construct the analog combining operation shown in FIG.

[0313] 14.6.2 Analog Affine Transformations The analog combining operation described above could take a variety of forms (including, for example, affine transformations). Figure 43 shows an analog circuit that implements an affine transformation on a vector v. Suppose we want to implement a transformation Av for some matrix A. Each row of A is a multiple of v, and this multiplication can be performed using layers of resistors as shown in Figure 43. Assuming A is a dense matrix, this involves adding d 2 resistors because there are d resistors associated with each row of A, and A has d rows. If A is less dense, fewer resistors can be used. As an optional feature, adding a voltage vector b to the output gives the general affine transformation Av+b.

[0314] 14.6.3 Higher-Order Join Operations The analog affine transformation described above can be viewed as a linear combination operation because it results in a weighted sum of terms that are linear in the components of v.

[0315] Alternatively, one can consider analog combining operations involving higher-order terms (e.g., terms that are quadratic or cubic in the components of v). One advantage of such an approach is that it allows the representation of more complex probability distributions P(x), including distributions that are difficult to sample using standard digital methods.

[0316] One way to achieve higher order terms is to replace each resistor in Figure 43 with a different circuit element (eg, an n-port device, where n>2) to achieve higher order coupling.

[0317] For example, consider a three-port device, the transistor. Figure 44 shows one way to create higher-order coupling using a transistor. Here, the voltage element v j can be used as the gate voltage for the transistor on the kth arm of the circuit (k ≠ j). In this sense, v j is the voltage v k This controls the resistance that occurs and results in higher order coupling.

[0318] More generally, other multi-port devices, including chains of multiple transistors, can be used to generate other higher order couplings.

[0319] 14.7 Maxwell's Demon Approach to Distributions Represented by Neural Networks In many applications, the log probability of being sampled is proportional to the output of the task. In this case, the probability can be expressed as:

number

[0320] Specifically, Bayesian estimation samples the posterior distribution.

number

number

number

number

number

[0321] In such cases, the model

number

number

number

number

[0322] Thus, after some time T, we have:

number

number

number

number

[0323] After a time increment T, the voltage vector is given by the mass inverse matrix

number

[0324] Figure 45 shows a scheme for a Maxwell's demon device 4110 for sampling from a distribution represented by a neural network. A user 41 first selects a model f θis compiled onto a hardware analog neural network (ANN) 4112 (along with its gradient if a closed form exists). The output of the gradient is then input into a momentum integrator device 4120, whose output is then input into a position integrator device 4130. The momentum integrator device 4120 is a device consisting of d analog integrators that output voltages representing momentum, as shown in FIG. 45. The position integrator device 4130 also contains d analog integrators and represents the position variable. The voltage vector corresponding to position

number

[0325] 14.7.1 Calculating the gradient Calculating the gradient in analog is difficult. A simpler way to estimate the gradient of a function is to perform finite differences. This can introduce noise into the gradient calculation. In the case of HMC, this noise does not degrade performance because of the Metropolis-Hastings step, which is another source of uncorrelated noise, as explained below for stochastic gradient HMC. Specifically, the finite difference gradient can be written as the exact gradient with some Gaussian noise with zero mean and variance V added.

number

[0326] Thus, gradients can introduce errors into the dynamics, which can be considered as an additional noise source on top of the noise coming from the hardware.

[0327] Hardware for Gaussian sampling via 15 precision matrices The device shown in Figure 45 may be slightly modified to perform Gaussian sampling. One difference is that the gradient calculator calculates the negative gradient of the multivariate Gaussian distribution.

number

number

[0328] Figure 46 shows an alternative device for sampling a multivariate Gaussian distribution with a given mean and precision matrix. The gradient calculator can be replaced by a resistor layer 4210 whose resistances correspond to the elements of matrix Q. To make this practical, an op-amp can be added to each output port of the resistor layer to prevent values ​​from disappearing. Then, a negative gradient

number

number

[0329] 16 A Multibody Coupled Approach to Sampling Non-Gaussian Distributions 16.1 Motivation and Overview 16.1.1 Beyond two-body bonds In the previous section, we described one method for extending beyond Gaussian distributions based on the use of Maxwell's demon (MD). In this approach, we assumed that the direct coupling between unit cells is a two-body coupling (e.g., the capacitive coupling mentioned above).

[0330] In this section, we describe an alternative approach beyond Gaussian distributions that involves modifying the direct coupling between unit cells. Specifically, in this approach, the direct coupling between unit cells is designed to be many-body (e.g., three-body, or more generally, n-body).

[0331] This approach introduces n-body terms in the potential energy function U(x) and higher-order cumulants (e.g., κ3 or κ4 as discussed in Section 13.4.3) in the probability distribution p(x).As discussed earlier in Section 13.4.3, there is a mathematical relationship between the introduction of many-body couplings and the introduction of higher-order cumulants in the distribution p(x).

[0332] To construct a three-body coupling, an electrical device with three ports (instead of the two ports associated with the capacitive coupling) is used. More generally, to construct an n-body coupling for some integer n>2, an electrical device with n ports is used. The voltages associated with n different unit cells (i.e., n different LC oscillators from the HMC hardware) are input to the n ports of the electrical device responsible for the coupling. This n-port device is called a coupling device.

[0333] The following details specific physical circuit elements that can be used to construct n-port coupled devices. For example, as detailed below, a network of transistors or a network of varicaps (also known as varactors) can be used for this purpose.

[0334] This n-body coupling approach can be used independently or together with the Maxwell's demon approach. In addition, this approach can be used independently or together with the two-body capacitive coupling described in the previous section.

[0335] Varying the direct coupling between unit cells has the advantage that it is difficult to simulate digitally, as will be discussed in more detail below.

[0336] 16.1.2 Effective digital simulation of multibody coupling In principle, many-body coupling can be difficult to simulate by standard digital computers. One way to understand this is that a general k-th order coupling can be expressed as d k The equation is mathematically represented by a k-th order tensor with elements d k For example, a cubic combination may contain d multiplication operators for computation by a digital computer. 3operations are involved. Therefore, such multi-body coupling can be difficult for a digital computer to simulate.

[0337] However, there are situations in which digital devices can efficiently simulate such couplings: namely, when a multi-body coupling is generated by a series of two-body couplings. In this case, the multi-body coupling is effective rather than true, meaning that the multi-body coupling can be efficiently represented by a small number of two-body couplings (in contrast, a true multi-body coupling refers to a coupling that cannot be efficiently represented as a small number of two-body couplings).

[0338] The Maxwell's demon approach described in Section 14 primarily addresses the case of effective, rather than true, many-body coupling. One exception is the approach described in Section 14.6.3, which involves transistors to generate higher-order couplings that are in fact true.

[0339] If we ignore this exception and instead consider the analog affine transformations of Section 14.6.2, we can see that a digital computer can 2 We expect to be able to simulate an analog Maxwell's demon with complexity that scales with d 2 This is because it contains operations. 2 The complexity of (1) persists even when the analog Maxwell's demon generates many-body couplings through nonlinear circuit elements, because this version of the analog Maxwell's demon involves effective, rather than true, many-body couplings.

[0340] When designing an analog system that is difficult to simulate digitally, it is useful to use true multi-body coupling rather than effective multi-body coupling. Here, when we say "difficult to simulate digitally," we mean that the complexity of using a digital computer makes it difficult to simulate digitally. 2 A scale change larger than (for example, d 3This means that the image may need to be scaled.

[0341] 16.1.3 True Many-Body Bonds A brief explanation of what is meant by true multi-body coupling is that true multi-body coupling can be viewed as a multi-body coupling of related circuit elements. This is in contrast to a multi-body coupling that is constructed from circuit elements with few-body couplings (e.g., two-body couplings).

[0342] From a computational perspective, a true many-body bond can be considered when no additional computational resources are required to construct the bond. For example, multiple two-body bonds can be used in combination with an auxiliary system to generate a three-body bond, but this is not considered a true three-body bond because it requires the use of additional resources (e.g., the auxiliary system). In contrast, a three-body bond can be generated using a single transistor and no additional auxiliary systems, which is considered true.

[0343] 16.1.4 Capacitive versus Resistive Multibody Coupling As previously mentioned, creating a multibody coupling device requires the use of circuit elements with many input and / or output ports.

[0344] One goal is to design a sampling probability distribution P(x), which corresponds to designing a potential energy function U(x), since P(x) ∝ exp(-U(x)).

[0345] In terms of physical circuit elements, potential energy is associated with the capacitive elements of a circuit, since capacitors can store energy, while resistive elements in a circuit are associated with damping or friction and technically do not contribute directly to the potential energy function.

[0346] Therefore, for the purposes of designing P(x), we are more interested in capacitive elements than resistive elements: purely resistive circuit elements are not very useful for designing the probability distribution P(x).

[0347] So, in addition to needing n-port devices, we also need some capacitance in those devices. In particular, we want the value of that capacitance to be voltage-controlled, i.e., controlled by the other input port to the device. In this sense, we are essentially interested in n-port devices that function as voltage-controlled capacitors.

[0348] 16.1.5 Elements acting as voltage-controlled capacitors A variety of different approaches can be used to construct a voltage-controlled capacitor.

[0349] Transistors have some capacitance that can be affected to some extent by the gate voltage. Therefore, transistors are one option for voltage-controlled capacitors. None of the three terminals of the transistor are grounded, and as a result, strategies must be devised to avoid occupying excessive space on-chip. Fortunately, transistor technology based on fully depleted silicon-on-insulator (FDSOI) provides a way to keep the on-chip area small in this case.

[0350] Varicaps (also known as varactors) offer another option. Each varicap is a two-port device, but two varicaps can be placed in opposite directions in a back-to-back configuration, as shown in Figure 19. By inserting a voltage wire between the two varicaps, the overall capacitance of the system can be adjusted. As a result, the device shown in Figure 19 functions as a voltage-controlled capacitor with a total of three ports, thus possessing the desired characteristics for our purposes. We will call this a three-port varicap device.

[0351] Another strategy is to use a switching capacitor bank, as shown in Figure 18. An external voltage signal controls the capacitance of this capacitor bank, which can be considered a voltage-controlled capacitor. One drawback of this approach is that it requires real-time switching on and off of the circuit elements in Figure 18, which can result in energy dissipation (i.e., decay) during the HMC process.

[0352] Transistors, three-port varicap devices, and switching capacitor banks can all be considered voltage-controlled capacitors. Each of these voltage-controlled capacitors can function as a building block for higher-order coupling. That is, by chaining together a set of these three-port devices, an n-port device can be constructed. These n-port devices therefore enable the design of n-body coupling.

[0353] 16.1.6 Adding adjustability or switches to couplings. Adding tunability to three-body or n-body couplings allows hardware devices to represent a wide range of probability distributions with n-body terms.

[0354] Adding tunability to an n-body coupler can be difficult, but a simpler approach is to add switches to the coupler. The switches turn the n-body coupling off or on. The user can specify a probability distribution by choosing the direction for these switches. This provides a rough degree of tunability that is relatively easy to implement.

[0355] 16.2 Three coupled unit cells Figure 47 shows a schematic diagram of three-body coupling via a voltage-controlled capacitor.

[0356] Now, in each of the three unit cells, we take the voltage (referenced to ground) associated with each capacitor. 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 input to the transistor's gate, one to the source, and one to the drain. Given that there are many possible physical implementations for this VCC (as discussed above), we will write the three-body coupler as a voltage-controlled capacitor (VCC).

[0357] In the next subsection, we analyze how this coupling affects the sampled probability distribution.

[0358] 16.2.1 Mathematical description of three merged cells For two capacitively coupled unit cells, the contribution to the potential energy has the form U=V1V2C 12 Here, C 12 Assume that V is a function of the voltage V3 of the third unit cell. This is the case shown in Figure 47. Therefore, the contribution to the potential energy is U = V1V2C 12 (V3) where C 12 (V3) is some function determined by the physical properties of the voltage controlled capacitor.

[0359] For example, assume that this is a linear function, which is often true over a narrow range of voltages. In this case,

number

number

[0360] So in this case we arrive at two terms in the potential energy: (1)C 123 and (2)

number

[0361] In Section 13.4.3, we established a mathematical relationship between the three-body terms in the potential energy and the third-order cumulant (denoted by κ3) in the corresponding probability distribution P(x). Thus, equation (98) results in a probability distribution P(x) that is non-Gaussian with a third-order cumulant. Therefore, using the three-body coupling approach, we can design the probability distribution P(x) to be non-Gaussian.

[0362] 16.3 Combining four or more unit cells Figure 48 shows a circuit diagram of how four unit cells can be coupled by a four-body coupler. The voltages V3 and V4 at the third and fourth unit cells are input to a voltage multiplier, the output of which is the product V = V3V4. The voltage multiplier can be either analog or digital, as described in more detail below. This voltage V is then applied to the gate of a voltage-controlled capacitor, creating a capacitance C between unit cells 1 and 2. 12 (V3V4) is produced.

[0363] C 12 (V3V4) can be a complex function, but it is only in a certain domain.

number

number

[0364] Thus, this implementation provides both four-body and two-body coupling.

[0365] In Section 13.4.3, we established a mathematical relationship between the four-body terms in the potential energy and the fourth-order cumulant (denoted by κ) in the corresponding probability distribution P(x). Thus, using this approach, we obtain a probability distribution P(x) that is non-Gaussian with a fourth-order cumulant.

[0366] The approach shown in Figure 48 can be straightforwardly extended to couple more than four cells (n>4). This involves feeding additional inputs to the voltage multiplier, which multiplies n-2 voltages to ultimately obtain the voltage applied to the gate of the voltage-controlled capacitor in the n-body coupler. Due to the mathematical relationship between n-body coupling in the potential energy and n-th order cumulants in the probability distribution P(x), the resulting distribution is non-Gaussian with n-th order cumulants (among other possible lower order cumulants).

[0367] 16.4 Analog versus Digital-Analog Approach to Many-Body Coupling 16.4.1 Fully analog approach In principle, the components shown in Figures 47 and 48 could be analog devices, since both the voltage-controlled capacitor and the voltage multiplier could be constructed from analog components. This has the advantage of low latency, since no communication with a digital device is involved. Furthermore, multiplying voltages incurs some computational overhead in digital devices (but not in analog devices), which could lead to even greater speedups.

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

[0369] 16.4.2 Digital-Analog Approach A convenient approach to multi-body coupling is to use the existing infrastructure already in place to implement HMC hardware with two-body coupling: Section 8 discloses the use of voltage-controlled capacitors in the two-body case, using digital hardware to set the capacitance value.

[0370] In the case of a three-body coupling, a voltage measurement V3 can be taken from the third unit cell, sent as an input to the digital hardware, and the digital hardware can output a voltage proportional to V3, which can be applied to the capacitor connecting unit cells 1 and 2. Thus, the digital hardware acts as an intermediary between unit cell 3 and the capacitor between unit cells 1 and 2.

[0371] For the four-body bond, the same approach can be taken, except that the digital hardware calculates the product of voltages V3 and V4. The digital hardware outputs a voltage proportional to V3V4 and applies this voltage to the capacitor connecting unit cells 1 and 2.

[0372] Once again, one advantage of this approach is that it utilizes the infrastructure already in place to implement HMC with two-body coupling.

[0373] 16.5 Strategies for Local Connectivity In principle, hardware can be designed to allow full connectivity (global connectivity), where n-body connections allow d n The number of connections grows rapidly due to scaling. There is a trade-off between hardware complexity and speedup that must be considered: more connections increase the potential for speedup, but the hardware becomes more complex.

[0374] Because of this trade-off, it may be desirable to consider local connectivity. One approach to local connectivity is to leverage the existing infrastructure used for two-body coupling. That is, one can use the digital-analog approach to many-body coupling (discussed immediately above), and any connectivity structure already in place for two-body coupling can be used as the same structure for n-body coupling. In this case, the connectivity structure for the n-body coupling coincides with the connectivity structure for the two-body coupling. Thus, if the two-body coupling is local, the n-body coupling is also local.

[0375] An alternative approach is to design n-body connectivity structures from scratch (i.e., without using two-body connectivity structures). Here, we present some strategies for this approach for the special case of three-body bonding, as shown in Figures 49 and 50. Figure 49 shows how to fit three-body bonding to a square lattice, and Figure 50 shows how to fit three-body bonding to a hexagonal lattice.

[0376] 17 Speedup Analysis 17.1 Runtime Analysis of Digital HMC 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 incorporated by reference in its entirety for all purposes. As previously mentioned, the system performs HMC by integrating Hamilton's equations.

number

[0377] Since these are vector differential equations,

number

[0378] Numerically, a common integration method is the leapfrog integrator, which gives updates to position and momentum as follows:

number

[0379] In the optimal case, the numerical complexity of integrating such an equation is n t is the number of time steps, and t ) However, computing the right-hand side of each equation to perform a numerical integration can be more expensive. In fact, the computation of the log-probability depends on the form of the probability distribution being sampled. The mass matrix M is a dxd matrix in the general case, so the cost of its inversion is d 3 However, in most implementations, M is considered to be diagonal, reducing the cost to d.

[0380] Without being bound by any particular theory, in order to keep the numerical integration error low enough, the time step size ε used in the numerical integration is d 1 / 4 Therefore, the gradient of the log-probability is efficiently obtained, and its inversion is computed in d steps, given that the mass matrix M is diagonal. 5 / 4 (In practice, when the probability distribution is represented by a neural network, the complexity of evaluating the gradient grows polynomially with the number of parameters of the model as well as the input size.)

[0381] In this sense, d 5 / 4 is the optimal scaling for a digital implementation of HMC. If the sparsity structure of M is unknown, then scaling is O(d 3 ) steps. In the next section, we consider specific speedups that can be achieved by using hardware.

[0382] 17.2 Execution Times for Analog HMC In the general case, the speedup gained by implementing differential equations in analog systems comes from two main areas: (i) there is no need to discretize the equations, and therefore no step sizes, which means that numerical errors resulting from step sizes are eliminated (in electrical hardware there are additional noise sources, e.g., thermal and shot noise, which can also affect performance); and (ii) to solve the differential equations numerically, the update terms are calculated, e.g., by leapfrog integration methods. In analog, the system naturally follows the differential equations; therefore, there is no need to perform matrix multiplications and inversions, significantly reducing the number of operations.

[0383] A further subtlety is the computation of the gradient: depending on the form of the probability distribution, this may affect the speedup obtained by analog HMC devices; this is discussed further in Section 14.7 above.

[0384] Therefore, in the ideal case, the number of digital operations to implement an analog SDE solver for a d-dimensional object is O(1). However, it is necessary to communicate digitally with the analog solver, for example, to set inputs and measure outputs. For sufficiently small data (on the order of tens of megabytes), the data can be loaded into the analog system in O(1) steps. However, as the data size increases, linear scaling in the dimension d emerges. Therefore, we consider the best-case scaling of an analog solver to be O(1).

[0385] 17.3 Speedup for Gaussian Sampling As shown above, if the probability distribution being sampled is a multivariate Gaussian distribution, the form of the gradient can be calculated analytically. Recall that the probability distribution is defined as follows:

number

[0386] The gradient of the log-probability can be expressed in affine form.

number

[0387] The gradient calculation therefore amounts to performing a vector-matrix multiplication, where the matrix involved is the inverse of the covariance matrix. This results in a running time of O(d 3 ) indicates that

[0388] 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

number

number

number

number

number

number

[0389] 17.3.1 Speedup for dense covariance matrices For dense covariance matrices, the Cholesky decomposition is O(d 3 ) steps. In many applications (e.g., Gaussian distributed processes), the covariance matrix is ​​updated and resampled many times, which changes the structure of the covariance matrix. This means that in many applications, the covariance matrix is ​​assumed to be dense. In this case, if the analog system is nearly perfectly connected, i.e., the connectivity k scales with d, then the covariance matrix can be calculated in O(d 3 ) speedup. Otherwise, updating the capacitance values ​​so that the capacitance matrix matches the covariance matrix may suppress most values ​​in the covariance matrix due to limited connectivity, causing the acceptance rate of the HMC proposal to scale 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 sparsed by a preprocessing step. This step takes O(d 2 ) operations are performed, resulting in a sparse and potentially banded matrix. This step may be performed for analog as well as digital approaches and can be a bottleneck. Therefore, if connectivity k~O(1), using analog HMC for Gaussian sampling may not provide any speedup.

[0390] 17.3.2 Speedup for diagonal and sparse covariance matrices If the covariance matrix is ​​diagonal, or equivalently, if the random variables

number

[0391] If the covariance matrix is ​​sparse, the most efficient strategy is to convert the sparse matrix to a band matrix, which has k central diagonals nonzero along with the main diagonal (if k=1, the matrix is ​​tridiagonal). nz By sorting the non-zero elements, a banded matrix may be obtained in linear running time using the Cuthill-Mackie process.

[0392] This leads to the special case where the sparse matrix is ​​a band matrix, which has k central diagonals nonzero along with the main diagonal (for k=1, the matrix is ​​tridiagonal). In this case, Cholesky factorization takes O(k 2 d) steps. The final matrix-vector product used in Cholesky sampling is O(kd) because the band structure of the matrix is ​​preserved in the Cholesky decomposition. Due to the connectivity of HMC analog devices, it seems most natural to implement a covariance matrix that is banded. Therefore, in this case, O(k 2 d) The increased speed achieved by analog HMC is expected.

[0393] 17.4 Speedup for distributions with higher-order correlation tensors Here, we consider potential speedups in the context of perturbation distributions with nontrivial higher-order cumulants (beyond quadratic). For this reason, we focus on d-dimensional distributions represented through the Gram-Charlier A-series introduced in equation (78). Taking the logarithm of such a distribution results in a polynomial in d variables. For each k-th order cumulant included in the perturbation expansion, a simultaneous k-th order polynomial is introduced in the potential energy function. The coefficients of this polynomial are represented through a k-th order tensor. Thus, if the highest nontrivial cumulant is of order K, then the potential energy function will be a leading k-th order polynomial.

[0394] The potential speedup resulting from higher-order cumulants is expected to arise from the fact that these are introduced by higher-order couplings or, equivalently, higher-order tensor interactions. Calculations involving these higher-order tensors can be used to compute the gradient of the log-probability relative to the HMC.

[0395] In analog hardware, the general kth order cumulant is introduced by a kth order combination (represented by a kth order tensor). Therefore, to digitally compute the gradient of the log-probability, the number of operators is O(d k ) In contrast, in analog hardware, the dynamics introduced by these interactions occur naturally in O(1). This results in O(d k ) yields a potential speedup. This analysis applies when true many-body coupling exists, as discussed above.

[0396] When many-body coupling is introduced by chaining two-body couplings, the operation of digitally simulating these dynamics becomes O(d 2 ) again, in the analog case, the physical dynamics may be realized in O(1), rather than the expected O(d 2 ) will be faster.

[0397] 17.4.1 Speedup for distributions represented by neural networks As explained in Section 10, a certain number of speedups result from this approach. Ad 1 / 4 The speedup of comes from not needing to discretize time as is done in digital HMC, leading to adjusting the number of time steps as a function of the dimension. The gradient computation is at most O(d 3 ) steps, but here it is performed analogously by finite differences. Here, the errors introduced by finite differences are not as detrimental as in other applications thanks to the Metropolis-Hastings step. Finally, the computation of the accept / reject steps and the data upload / download may all be performed analogously, keeping the number of O(d) operations down as mentioned in Section 16.2.

[0398] Overall, therefore, for an analog HMC with a probability distribution represented by a neural network, d 13 / 4 It is expected that the speed will increase.

[0399] 18 Running Alternative Monte Carlo Processes on Hardware 18.1 Realistic Analysis of Damped Hardware The protocols for HMC sampling presented so far utilize trajectories of analog circuit dynamics that are assumed to be conservative. These trajectories are then borrowed in a computational protocol that yields samples from a target probability distribution. Here, we provide a more realistic analysis of physical circuits and suggest alternative processes that can be utilized when the actual circuit dynamics deviate significantly from the intended Hamiltonian dynamics.

[0400] Consider the LC oscillator cell of Figure 3 used for one-dimensional Gaussian sampling. As a first approximation to the realistic dynamics of this circuit, we connect resistors in series with both the capacitor and the inductor. These resistors are represented by the resistance R C and R L , the circuit has an effective resistance R=R C +R L In addition, assume that a voltage signal V(t) is applied to the circuit and that V(0)=0.

[0401] Let the charge q on the capacitor C serve as the position coordinate of the system, and the magnetic flux LI through the inductor be the conjugate momentum coordinate.

[0402] The time derivative of the momentum LI(t) of the circuit is governed by the following differential equation:

number

[0403] The integral of I(τ) up to time t is equal to the charge q(t) stored on the capacitor at this time. Therefore, the ODE can be rewritten as:

number

[0404] From this relationship between the integral in (108) and the charge on the capacitor, we see that the time derivative of q(t) is I(t).

number

[0405] Putting together the coupled dynamic equations (109) and (110), we arrive at:

number

[0406] Labeling the momentum p=LI and using dot notation, the above system of equations becomes:

number

[0407] The term -Rp represents a dissipative force with linear momentum, which does not exist in a conservative system.

[0408] We now elaborate on the nature of the V(t) voltage source. While this source signal can be a deterministic signal provided by a user, we focus here on the stochastic voltage signal imposed on the system via an underlying noise process. As a first approximation to a realistic noise model for RLC cell dynamics, V(t) arises from Johnson thermal noise. This noise arises from the thermal motion of a collection of macroscopic electrons through an effective resistor with resistance R held at some temperature T.

[0409] Johnson noise can be adequately modeled using a one-dimensional Brownian motion W(t) and a diffusion constant D, i.e., with a slight misuse of notation, V(t)~N(0,Ddt)=DdW. Introducing this into equation (113), we see that the RLC dynamics is modeled by a coupled stochastic differential equation (SDE).

number

[0410] What remains is to derive the diffusion constant D. For this, we apply the fluctuation-dissipation theorem. According to the fluctuation-dissipation theorem that circuit systems obey, the diffusion constant D of Johnson noise is related to the coefficient of dissipation power R (which is simply the effective resistance), the inductance L, and the temperature of the resistor. In fact, D is given as follows:

number

number

[0411] Therefore, as assumed above, the circuit dynamics is a stochastic process, not a deterministic process. Furthermore, stochastic dynamics converge to equilibrium. Suppose the circuit is initialized at time t = 0 with an arbitrary position q0 and an arbitrary momentum p0. Convergence to equilibrium means that the position and momentum coordinates at a long time after this initialization (t >> 0) are distributed according to the stationary distribution π(q, p) of the SDE (118). This distribution can be solved as follows:

number

[0412] Marginalizing with respect to the momentum coordinate p, the equilibrium distribution of the position coordinate (charge on the capacitor) is given as a Boltzmann distribution.

number

[0413] Here we provide a protocol for sampling from arbitrary multivariate Gaussian distributions using the stochastic circuit dynamics employed above.

[0414] 18.2. Protocol for Univariate Gaussian Sampling Recall that RLC systems converge to equilibrium regardless of the initial values ​​of position and momentum. The equilibrium of a system is characterized by a stationary distribution π(q,p) in phase space. We want to prepare the RLC circuit at this equilibrium state, which corresponds to leaving the circuit in contact with a heat bath for a long enough time. Once the circuit is prepared at equilibrium, repeated measurements of the circuit's position q and momentum p (the charge on the capacitor and the magnetic flux through the inductor, respectively) will be distributed according to π(q,p). By measuring only the position, we effectively marginalize over the momentum and sample from the system's Boltzmann distribution, which is a zero-mean Gaussian with variance kTC. Adjusting the value of the capacitance C changes the variance of the Boltzmann distribution we sample. For any σ 2 >0

number

[0415] The variance of the sampled Boltzmann distribution can be adjusted, but the mean of this distribution remains zero.

number

[0416] 18.3 Extension to multivariate Gaussian sampling Although the realistic analysis provided in 18.1 was performed for a single RLC cell, a similar analysis follows for a network of capacitively coupled RLC cells. Assuming N RLC cells are coupled in pairs via coupling capacitance bridges, the system can be expressed as an N-dimensional position vector

number

number

[0417] Just as the coordinates (q,p) of a single unit cell obey the underdamped Langevin dynamics in two-dimensional phase space, the coordinate system

number

number

[0418] Fortunately, we do not need to know the explicit forms of the drag and diffusion matrices to know that the dynamics of the coupled RLC network will converge to equilibrium anyway. Just as the equilibrium distribution of a single cell (120) does not depend on the friction and diffusion constants, the equilibrium distribution of the dynamics (122) does not depend on either the drag or diffusion matrices. If each resistor in the RLC network is held at temperature T, then the equilibrium distribution in phase space is defined as follows:

number

[0419] When we marginalize in the momentum coordinate, we find that the resulting Boltzmann distribution in N-dimensional position space is given by

number

number

number

number

[0420] Then, any multivariate Gaussian distribution

number

number

number

[0421] 18.4 Stochastic Gradient Hamiltonian Monte Carlo We present a realistic analysis of the single-cell oscillator cell used in the HMC protocol and show that this circuit dynamics undergoes underdamped Langevin dynamics as opposed to Hamiltonian dynamics. We demonstrate a process for utilizing digitally simulated underdamped Langevin dynamics. A similar sampling protocol can be implemented using the underdamped circuit dynamics presented in Section 18.1.

[0422] HMC involves directly computing or accessing the gradient of the potential energy function to digitally simulate the Hamiltonian dynamics of a system. With 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).

[0423] 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 using less data than the entire dataset. However, as a result, the dynamics introduced are no longer caused by the desired distribution. Therefore, an additional friction term is included to counteract the effect of the noise injected by the stochastic gradient.

[0424] Explicitly, in SGHMC, the gradient is

number

number

number

[0425]

number

number

number

number

number

number

number

number

[0426] 18.4.1 Approach to SGHMC Consider performing Gaussian process regression, where the posterior distribution is learned from mini-batches of training data, rather than the entire training set. Mapping the logarithm of this mini-batched posterior distribution to the potential energy function of a realistic oscillatory circuit yields a noisy version of the true potential energy function. The potential energy function of an oscillatory circuit is determined by the inverse capacitance matrix C -1 Since the noise manifests itself 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 one capacitance value will be approximately Gaussian distributed around the true value.

[0427] To reduce the impact of this misspecification noise and converge to the true target distribution, it is not possible to adjust the friction (resistance) from a realistic oscillator circuit to match the noise profile, as is done in digital SGHMC. While the friction can be adjusted, due to fluctuation-dissipation relationships, this will proportionally adjust the thermal noise of the system. This is not necessarily true for digital dynamics, but may be unavoidable in circuit dynamics. Thus, the friction term in a physical oscillator should reduce the impact of thermal noise as well as misspecification noise, but it can only reduce the contribution from thermal noise.

[0428] To avoid this problem, we can operate in the high friction / resistance region. When this occurs, a large amount of thermal fluctuation occurs, which becomes dominant compared to the fluctuations caused by the misspecification noise. Therefore, the noise profile describing both the thermal noise and the misspecification noise approximately satisfies the fluctuation dissipation by the friction term, since it is primarily thermal noise. When this relationship is approximately satisfied, equilibrium is reached. Because the noisy portion of the potential has been absorbed into the noise term, the Boltzmann distribution becomes a negative exponential function of the true, unminibatched potential, i.e., the true posterior distribution. The true multivariate posterior is then sampled using the protocol specified in Section 18.3 in the high-resistance limit.

[0429] 19 Hardware for Langevin Monte Carlo There is a variant of HMC where the number of leapfrog steps is set to L = 1. This is known as Langevin Monte Carlo and is defined as follows:

number

number

number

number

[0430] 19.1 Unit cells for LMC hardware The unit cell (i.e., building block) for the LMC hardware is shown in FIG.

[0431] This differs from the unit cells used for HMC because they can operate without inductors. Avoiding inductors is a useful feature of LMC hardware because inductors are difficult to implement on-chip (e.g., they occupy a significant amount of area on the chip).

[0432] On the other hand, the unit cell for LMC has a stochastic noise source, which corresponds to a voltage source whose voltage output is Gaussian white noise. There are various ways to implement a stochastic noise source for this purpose, such as the following possible ways: 1. Analog thermal noise from resistors (also known as Johnson-Nyquist noise). 2. Analog shot noise from diodes (e.g., Zener diodes). 3. Digital pseudorandom noise generated from a field programmable gate array (FPGA) with analog filtering.

[0433] Each of these possible noise sources can act to generate the δv voltage source shown in Figure 51. Both the analog thermal noise and the analog shot noise can be amplified by an amplifier to increase the magnitude of the associated voltage fluctuations.

[0434] Each unit cell may also contain a constant voltage source (although this voltage source is not explicitly shown in Figure 51), which may be useful to move the mean vector μ associated with the probability distribution p(x) away from the value μ=0.

[0435] 19.2 Coupled unit cells for LMC hardware Consider two alternatives for coupling unit cells. One approach is to use a resistive bridge, as shown in Figure 52A. Another approach is to use a capacitive bridge, as shown in Figure 52B.

[0436] For resistive coupling, the equation of motion for the voltage vector v for two coupled cells is:

number

number

[0437] The self-resistance matrix R, the capacitance matrix C, and the conductance matrix J are introduced.

[0438] In the case of capacitive coupling, the equation of motion is:

number

number

[0439] 19.3 Acceptance or Rejection of Samples As with HMC, a Metropolis-Hastings (MH) step can be added to the LMC process, which involves the following steps: ·Measure the voltage across each of the capacitors in the unit cell of the LMC hardware. ·This voltage vector v is input into a digital processor. The energy associated with this voltage vector is quantified in a digital processor. Based on this energy value, the accept / reject step of the MH process is applied to this sample.

[0440] This allows the digital processor to accept or reject the proposed sample as a judgement of sample quality (via the MH condition).

[0441] 19.4 Sampling Gaussians with LMC Hardware For a Gaussian distribution, the slope of the logarithm of the probability distribution is given by an affine function, specifically: ∇logp(x)=-Σ -1 (x-μ) (134) where μ is the mean vector and Σ is the covariance matrix.

[0442] Using this formula, equation (131) can be rewritten for the special case of a Gaussian distribution as follows:

number

[0443] This equation can be compared to equations (132) and (133), which are differential equations for LMC hardware with resistive and capacitive coupling, respectively. In fact, equation (135) has a similar form to these equations. Specifically, these equations can be used to solve for the relationship between the distribution parameters and the circuit parameters.

[0444] That is, the following relationship: (Σ,μ)⇔(R,J,C) (136) In the case of resistive coupling, the following relationship can be established: (Σ,μ)⇔(R,C) (137) For capacitive coupling, these relationships allow the user to convert the distribution parameters into circuit parameters for the LMC hardware. For resistive coupling, the exact relationships are given below:

[0445] Equations (132) and (133) assume that the mean vector is zero (μ=0), but a non-zero μ can be tolerated by adding a constant voltage source to each unit cell of the LMC hardware.

[0446] 19.4.1 Encoding the Covariance Matrix into Resistive Junctions For the special case of resistive coupling, we now present an analysis of how to encode the covariance matrix into the coupling.

[0447] For further tuning flexibility, add an additional resistive branch in parallel with each unit cell, as shown in FIG.

[0448] Now let us analyze the circuit in Figure 53. The voltage vector v = {v i , v i’} T Consider the following: where v i is the voltage across the i-th capacitor. We arrive at the following differential equation for v:

number

number

[0449] Consider the special case where R=RI and C=CI are proportional to the identity element I.

number

number

number

[0450] Therefore, the above equation can be used to encode the covariance matrix Σ into a J matrix.

[0451] This provides a recipe on how to store the covariance matrix in the parameters of an electrical circuit.

[0452] The above analysis generalizes from two unit cells to d unit cells, which essentially involves generalizing the J matrix to the following expressions for its matrix elements:

number

[0453] 19.5 Sampling Non-Gaussian Distributions with LMC Hardware The techniques developed for HMC to go beyond Gaussian sampling, described in Sections 14 and 16, also apply to LMC.

[0454] 19.5.1 Maxwell's Demon Approach The Maxwell's demon approach described in Section 14 can be applied to LMC hardware to allow sampling of non-Gaussian distributions.

[0455] Specifically, Maxwell's demon can be viewed as a voltage-controlled voltage source acting on a unit cell of the LMC hardware. The same circuit structures as shown in Figures 36, 38, and 42, for example, can be used, provided that the HMC unit cells are replaced with LMC unit cells (this essentially involves removing the inductors and adding stochastic noise sources to the unit cells, as previously described).

[0456] 19.5.2 Many-body coupling approach The many-body coupling approach described in Section 16 can also be applied to LMC hardware to allow sampling of non-Gaussian distributions.

[0457] Specifically, either a voltage-controlled capacitor (VCC) or a voltage-controlled resistor (VCR) can be used as the coupling device between unit cells. The discussion in Section 16 for HMC can be directly applied to LMC hardware, provided that the HMC unit cells are replaced with LMC unit cells. This includes the discussion of how to design three-body, four-body, or n-body couplings using voltage-controlled capacitors. In addition, the discussion in Section 16 of using either analog or digital-analog methods for many-body coupling is also relevant to LMC hardware.

[0458] 19.6 Alternative LMC hardware based on integrators Here, an alternative device for LMC hardware can be used as follows: Specifically, this alternative device uses an integrator instead of a unit cell, as shown in Figure 54.

[0459] The analog device for LMC in Figure 54 can be used to sample any probability distribution (e.g., one parameterized by a neural network as described in Section 14.7). To do so, the user must compile their representation of the probabilities into hardware. In the case of a neural network, this amounts to compiling the network into analog hardware and then computing its gradient. Alternatively, if there is a formula for estimating the gradient, this may be compiled into analog, and Maxwell's demon will then scale the input

number

[0460] 20 Increase in effective temperature due to noise injection One question regarding the feasibility of Gaussian sampling devices is whether the attenuation (caused by realistic electrical components) will be too great compared to the thermal noise, resulting in a very small measured voltage. While the measurement can always be amplified, there are concerns that adding an analog amplifier would increase the complexity of the device, and that large amplification achieved through digital post-processing could result in a loss of accuracy. A physically motivated way to overcome this problem is to simply increase the operating temperature of the device until the thermal fluctuations become large enough that amplification is no longer necessary. This approach comes with its own problems, such as the time it takes to reach the desired temperature as well as additional energy costs.

[0461] Increasing the effective temperature of a device through noise injection can be done by generating white noise with a digital pseudorandom number generator (e.g., coming from 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.

[0462] The circuit diagram with noise injection for the two unit cells case is shown in Figure 55. Assuming uncorrelated noise injection, the circuit equations take the form:

number

[0463] In fact, the effective temperature can be increased using uncorrelated noise in each cell. This means that samples can be taken without additional amplifiers or temperature control using the circuit described in the previous section of this paper. Assuming the injection of uncorrelated white noise with the same amplitude into each cell, the resulting covariance matrix for voltages and currents is: Σv=k0RC -1 ,ΣI=k0RL -1 (144) where C is Maxwell's capacitance matrix, L is the inductance matrix, R is the resistance matrix (proportional to the identity matrix), and κ is a number that scales with the injected noise amplitude and also determines the effective temperature.

[0464] In this form, the Maxwell capacitance matrix encodes the inverse of the covariance matrix rather than the covariance matrix itself, making it most suitable for applications where a precision matrix is ​​specified rather than a covariance matrix. However, it is still possible to use this method if, for example, you are given a covariance matrix that can be inverted numerically. There are some subtleties when trying to directly scale the capacitance matrix with the covariance matrix, but here we will cover the case where the capacitance matrix scales with the precision matrix.

[0465] 21 Thermodynamic Systems for Solving Linear Algebra Primitives Here, we also show that the same thermodynamic hardware described above can be used to accelerate key primitives in linear algebra. We exploit the fact that the mathematics of harmonic oscillator systems is affine (i.e., linear), and therefore we can map linear algebra primitives to such systems. We show that by sampling from the thermal equilibrium distribution of coupled harmonic oscillators, we can solve a variety of linear algebra problems. Specifically, we develop thermodynamic algorithms for the following linear algebra primitives: (i) solving the linear system Ax=b; (ii) inverting the matrix A -1 (iii) AΣ+ΣA T = I, and (iv) estimate 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.

[0466] 21.1 Solving systems of linear equations A well-known linear system problem is x∈R such that d of, Ax=b (145) Some invertible matrix A∈R d×d and non-zero b∈R d The goal is to find, given A' and A'. Without loss of generality, we can assume that the matrix A in equation (145) is symmetric and positive definite (SPD). To verify this, suppose A is not SPD. Then, A' = A' T A and b' = A T You can also specify that b and solve the equation. A'x'=b' (146)

[0467] 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. Since matrix A' is SPD by construction, if a method for solving symmetric linear systems is available, it is possible to find x' that satisfies equation (146). Then, by adding (A T ) -1Left-multiplying by x gives us Ax'=b, which means that x' is a solution to the original linear system. This may affect the total running time, but still allows for asymptotic scaling improvements with respect to digital methods (constructing an SPD system from a general one in this way squares the condition number, which affects performance). Therefore, in what follows, we will assume that A is SPD.

[0468] Now we connect this problem to thermodynamics. Consider a macroscopic device with d degrees of freedom described by classical physics. Assume the device has the following potential energy function:

number

[0469] Assume that the device reaches thermal equilibrium with the environment. The inverse temperature is β=1 / (k B T). In thermal equilibrium, the Boltzmann distribution describes the probability that an oscillator has a given spatial coordinate: p(x) ∝ exp(-βV(x)). Since V(x) is a quadratic form, p(x) corresponds to a multivariate Gaussian distribution. Therefore, in thermal equilibrium, the spatial coordinate x is a Gaussian random variable. x~N[A ―1 b, βA -1 ] (148)

[0470] The unique minimum of V(x) occurs at Ax-b=0, which also corresponds to the unique maximum of p(x). For a Gaussian distribution, the maximum of p(x) occurs at the first moment <x>Therefore, in thermal equilibrium, the first moment is the solution of a system of linear equations. <x>=A -1 b (149)

[0471] From this analysis, we can construct a thermodynamic protocol for solving the linear system shown in Figure 56. That is, the protocol involves realizing the potential in equation (147), waiting for the system to come to equilibrium, and then sampling x to find the mean of the distribution. <x>This average value can be approximated using a time average defined as:

number

[0472] The overall protocol can be summarized as follows: 1. Given a linear system Ax=b, set the potential of the device at time t=0 as follows:

number

number

number

number

number

number

number

number

[0473] To implement the above protocol,

number

[0474] In the overdamped case, we arrive at the following equation, which can be used in the above protocol:

number

[0475] In the underdamped case, the equations for the parameters are slightly different.

number

[0476] 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 algorithms mentioned above.

[0477] One potential approach to this hardware implementation is to use the components shown in Figure 17, which uses RC unit cells and resistive coupling between the cells. The voltage across the capacitor of each cell, v = (v 1, v 2, ...,v d ) the equation of motion is: dv=c -1 (-Jvdt+R -1 dw) (157) where C=diag(C 1, C 2, ...,C d ), R=diag(R 1, R 2, ...,R d ), where w is the uncorrelated Brownian motion, and the elements of J are given by

number

[0478] Therefore, there are two intra-cell resistors, and we can free the diagonal elements of the J matrix independently of the R matrix, with one resistor connecting each cell. Since R and C are diagonal matrices, we can set J = JR and γ = 1 / RC. By showing x = v and choosing J = A, we obtain:

number

number

number

[0479] 21.2 Estimating the Inverse of a Matrix The results in the previous section rely on an estimate of the mean of x and do not use the variation of x at equilibrium. By using the second moment of the equilibrium distribution, we can go beyond solving a linear system. For example, we can find the inverse of a symmetric positive definite matrix A. As mentioned earlier, the stationary distribution of x is N[A -1 b, β -1 A -1 ], which means that the inverse of A can be obtained by evaluating the covariance matrix of x. This can be done in a completely analog way, using a combination of analog multipliers and integrators. By setting b=0 for this protocol, <x>= 0, the stationary covariance matrix is, by definition,

number

[0480] To estimate this, time averaging is performed again after the system has come to equilibrium.

number

[0481] Analog components are expressed as a product x i (t)x j (t) can be evaluated for each pair (i,j), and d 2 An analog multiplier component is obtained. Each of these products is then input into an analog integrator component to calculate one element of the time-averaged covariance matrix.

number

[0482] The equilibration time is the same as in the linear system protocol, but the integration time is different because the covariance matrix generally converges slower than the mean. Here we present a detailed description of the inverse estimation protocol assuming ODL dynamics. 1. Given a positive definite matrix A, set the potential of the device at time t=0 as follows:

number

number

number

number

number

number

[0483] 21.3 Solving Lyapunov Equations This section concerns devices with controllable noise sources such that the covariance matrix of the noise terms can be chosen to be any symmetric positive definite matrix. T Since the x term is not included, we obtain the overdamped Langevin equation:

number

number

number

number

number

number

[0484] Analog multipliers and integrators are used to measure the time average.

number

number

[0485] 21.4 Estimating the determinant of a matrix The determinant of the covariance matrix appears in the normalization coefficients of the multivariate normal distribution, whose density function is:

number

[0486] Hardware capable of generating Gaussian distributions can estimate the determinant of a matrix, since this problem is equivalent to estimating the free energy difference, an application of stochastic thermodynamics. Recall that the free energy difference between the equilibrium states of potentials V1 and V2 is:

number

[0487] The potential is quadratic, V1(x)=x T A1x and V2(x)=x T Assume that A2x. Each integral is then simplified to the inverse of the Gaussian normalization factor.

number

number

[0488] This suggests that the determinant of matrix A1 can be found by comparing the equilibrium free energy with the potentials V1 and V2 (A2 has a known determinant) and then calculating: |A1|=e -2βΔF |A2| (177)

[0489] Fortunately, assuming we can measure the work done on the system when the potential V(x) changes from V1 to V2, we can find the free energy difference ΔF. According to the Dziarczynski equation, the free energy difference between the (equilibrium) states at the initial and final potentials is: e -βΔF = <e -βW > (178) where 〈 〉 denotes the average over all possible trajectories of the system between time t = 0 and time t = τ, weighted by their respective probabilities. This can be approximated as: [Table 1] where d is the matrix dimension, κ 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 digital SOTA, the complexity of solving symmetric positive definite linear systems, matrix inversion, Lyapunov equations, and determinant problems is for algorithms based on conjugate gradient methods, fast matrix multiplication / inversion, Bartels-Stewart algorithm, and Cholesky decomposition, respectively. ω ≈ 2.3 denotes the matrix multiplication constant. Average over N repeated trials

number

[0490] In summary, the determinant A1 is approximated as follows:

number

[0491] In practice, we may be interested in the logarithmic determinant to avoid computational overflow, which looks like this:

number

[0492] At least P δ With probability δ det To estimate the logarithmic determinant within , the waiting time is:

number

[0493] 21.5. Runtime Analysis Table 1 summarizes the various results obtained in comparison with the best 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 using an integrator-multiplier approach that guarantees convergence of the method.

[0494] 22 Conclusion While various embodiments of the present invention have been described and illustrated herein, those skilled in the art will readily envision various other means and / or structures for performing the functions and / or results and / or obtaining one or more of the advantages described herein, and each such variation and / or modification is deemed 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 herein are exemplary, and that the actual parameters, dimensions, materials, and / or configurations will depend on the particular application or applications 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. Accordingly, the foregoing embodiments are presented by way of example only, and it should be understood that, within the scope of the appended claims and their equivalents, embodiments of the invention may be practiced otherwise than as specifically described and claimed. The inventive embodiments of the present disclosure are directed to each individual feature, system, article, material, kit, and / or method described herein. Furthermore, any combination of two or more such features, systems, articles, materials, kits, and / or methods is included within the inventive scope of the present disclosure, if such features, systems, articles, materials, kits, and / or methods are not mutually inconsistent.

[0495] Also, various inventive concepts may be embodied as one or more methods, examples of which are provided. The acts performed as part of a method may be ordered in any suitable manner. Thus, while the exemplary embodiments are shown as sequential acts, embodiments may be constructed in which acts are performed in a different order than shown, which may include performing some acts simultaneously.

[0496] All definitions defined and used herein should be understood to supersede any dictionary definitions, definitions in documents incorporated herein by reference, and / or ordinary meanings of the defined terms.

[0497] The indefinite articles "a" and "an," as used in the specification and claims, unless expressly indicated to the contrary, should be understood to mean "at least one."

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

[0499] As used herein and in the 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" shall be interpreted as inclusive, i.e., the inclusion of at least one of a number or list of components, and optionally including multiple components and additional unlisted items. Only terms expressly indicated to the contrary, such as "only one" or "exactly one," or, when used in the claims, "consisting of," shall refer to the inclusion of exactly one component of a number or list of components. Generally, the term "or" as used herein shall only be interpreted as indicating exclusive alternatives (i.e., "one or the other, but not both") when preceded by exclusive terms such as "either," "one," "only one," or "exactly one." "Consisting essentially of," when used in the claims, shall have its ordinary meaning as used in the field of patent law.

[0500] As used herein and in the claims, the phrase "at least one" in reference to a list of one or more components should be understood to mean at least one component selected from any one or more components in the list of components, but does not necessarily include at least one of each component specifically listed in the list of components, and does not exclude any combination of components in the list of components. This definition also allows that components other than those specifically identified in the list of components to which the phrase "at least one" refers may optionally be present, whether related or unrelated to the specifically identified components. 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") can refer in one embodiment to at least one, optionally including multiple, A, in the absence of B (and optionally including components other than B); in another embodiment to at least one, optionally including multiple, B, in the absence of A (and optionally including components other than A); in yet another embodiment to at least one, optionally including multiple, A, and at least one, optionally including multiple, B (and optionally including other components), etc.

[0501] In the claims as well as in the foregoing specification, all transitional phrases such as "comprising," "including," "carrying," "having," "containing," "involving," "holding," "consisting of," and the like, shall be understood to be open-ended, i.e., to mean including, but not limited to. Only the transitional phrases "consisting of" and "consisting essentially of," shall be closed or semi-closed transitional phrases, respectively, as defined in the U.S. Patent Office Manual of Patent Examining Procedures § 2111.03.< / x> < / x> < / x> < / x>

Claims

1. 1. A system for sampling a multivariate distribution, comprising: a network of analog electrical circuits for sampling the multivariate distribution, each analog electrical circuit in the network of analog electrical circuits including at least one tunable passive electrical component; a digital controller operatively coupled to the network of analog electrical circuits, for adjusting the tunable passive electrical components according to parameters of the multivariate distribution and for sampling voltages at points within the network of analog electrical circuits; A system comprising:

2. The system of claim 1 , wherein the at least one tunable passive electrical component comprises at least one of a tunable capacitor or a tunable resistor.

3. The system of claim 1 , wherein the at least one tunable passive electrical component comprises a tunable inductor.

4. The system of claim 1 further comprising a voltage-controlled voltage source operatively coupled to the network of analog electrical circuits to modify the multivariate distribution.

5. The system of claim 4 , wherein the voltage-controlled voltage source uses an artificial neural network to relate an input voltage to the voltage-controlled voltage source to an output voltage of the voltage-controlled voltage source.

6. 6. The system of claim 5, wherein the input voltage is based on the voltage sampled by the digital controller, and the voltage-controlled voltage source is configured to apply the output voltage to the network of analog electrical circuits to vary a distribution of potential energy across the network of analog electrical circuits.

7. The multivariate distribution has at least third-order non-zero cumulants, and 2. The system of claim 1, comprising a multiport device operatively coupled to at least three of the analog electrical circuits in the network of analog electrical circuits to generate k-body terms in a potential energy distribution across the network of analog electrical circuits, where k is an integer greater than or equal to 3.

8. The system of claim 7 , wherein the multi-port device includes at least one of a transistor or a voltage-controlled capacitor.

9. The system of claim 1 , wherein each analog electrical circuit further comprises a stochastic noise source.

10. The system of claim 9 , wherein the stochastic noise source is one of a thermal noise source, a shot noise source, or a digital pseudorandom noise source.

11. The system of claim 1 , wherein the analog electrical circuits in the network of analog electrical circuits are coupled to each other via capacitive and / or resistive coupling.

12. The system of claim 1 , wherein the system is configured to generate samples from the multivariate distribution via a Langevin Monte Carlo algorithm and / or a Metropolis adjusted Langevin algorithm.

13. The system of claim 1 , wherein the multivariate distribution is a Gaussian distribution whose covariance matrix is ​​encoded in at least a portion of the parameters of the network of analog electrical circuits.

14. 1. A system for sampling a multivariate distribution, comprising: a network of coupled analog unit cells characterized by a total energy function that maps to the multivariate distribution, each coupled analog unit cell in the network of coupled analog unit cells comprising a digitally adjustable capacitor and a digitally adjustable voltage source and / or a digitally adjustable current source; a digital controller operably coupled to the network of coupled analog unit cells, initializing the digitally adjustable capacitor and the digitally adjustable voltage source and / or the digitally adjustable current source; repeatedly reading voltages and currents at nodes in the network of coupled analog unit cells and updating the digitally adjustable voltage sources and / or the digitally adjustable current sources, the voltages and currents corresponding to position and momentum of the total energy function; converting the voltages and currents into samples of the multivariate distribution; and A system comprising:

15. 15. The system of claim 14, wherein the total energy function is Hamiltonian and the network of coupled analog unit cells has dynamics used to implement a Hamiltonian Monte Carlo protocol.

16. 15. The system of claim 14, wherein the network of coupled analog unit cells has dynamics used to implement a Langevin Monte Carlo protocol.

17. 15. The system of claim 14, wherein the network of coupled analog unit cells includes one analog unit cell for each dimension of the multivariate distribution.

18. 15. The system of claim 14, wherein each analog unit cell in the network of coupled analog unit cells includes a digitally tunable inductor.

19. 15. The system of claim 14, wherein the analog unit cells in the network of coupled analog unit cells are capacitively coupled to one another.

20. 15. The system of claim 14, wherein the analog unit cells in the network of coupled analog unit cells are resistively coupled to one another.

21. 15. The system of claim 14, wherein the analog unit cells in the network of analog unit cells are inductively coupled to each other.

22. The system of claim 14 , wherein the multivariate distribution is a multivariate Gaussian distribution.

23. The system of claim 14 , wherein the multivariate distribution has non-zero cumulants higher than second order.

24. 1. A system for sampling a multivariate distribution, comprising: A gradient calculator, a momentum integrator device operably coupled to said gradient calculator; a position integrator device operably coupled to the momentum integrator device; and A system comprising:

25. 25. The system of claim 24, wherein the gradient calculator comprises an analog neural network.

26. 25. The system of claim 24, wherein the slope calculator comprises a network of resistive circuit elements.

27. 1. A thermodynamic system for sampling a target probability distribution, comprising: an analog dynamic system configured to evolve via the Langevin equations and / or Hamilton's equations of motion; a digital controller operably coupled to the analog dynamic system, the digital controller configured to receive proposed samples of the target probability distribution from the analog dynamic system at each time step of the Langevin equations and / or Hamilton's equations of motion, and to accept or reject the proposed samples; A thermodynamic system comprising:

28. 1. A method for inverting a matrix, comprising: uploading the elements of said matrix to respective values ​​of tunable circuit elements of a network of analog unit cells; allowing the network of analog unit cells to reach thermal equilibrium; calculating a covariance matrix that is the inverse of said matrix based on dynamic variables of said network of analog unit cells; The method comprising:

29. Calculating the covariance matrix extracting samples from the thermal equilibrium distribution of the dynamic variables of the network of analog unit cells; calculating the covariance matrix of the samples; 29. The method of claim 28, comprising:

30. 30. The method of claim 28, wherein calculating the covariance matrix comprises integrating the dynamic variables over time with a plurality of analog integrators.

31. A method for solving a system of linear equations represented by a matrix and a vector, comprising the steps of: uploading the elements of the matrix to respective values ​​of tunable circuit elements of the network of analog unit cells; uploading the elements of the vector to respective current or voltage sources within said analog unit cell; allowing the network of analog unit cells to reach thermal equilibrium; calculating average values ​​of dynamic variables of the network of analog unit cells, the average values ​​representing solutions to the system of linear equations; The method comprising:

32. 32. The method of claim 31, wherein calculating the average value comprises integrating the dynamic variable over time with a plurality of analog integrators.

Images