FIELD

The present techniques generally relate to sampling from statistical distributions and using the samples in computational tasks, such as those performed by a computational system comprising a quantum processor.

Sampling

Throughout this specification and the appended claims, the terms “sample”, “sampling”, “sampling device”, and “sample generator” are used. These terms are used herein in like manner to their corresponding uses in the arts of statistics and statistical analysis, and electrical engineering.

In statistics, a sample is a subset of a population, i.e., a selection of data taken from a statistical population. Sampling is the process of taking the sample, and typically follows a defined procedure. For example, in a population, database, or collection of objects, a sample may refer to an individual datum, data point, object, or subset of data, data points, and/or objects.

In electrical engineering and related disciplines, sampling relates to taking a set of measurements of an analog signal or some other physical system. Sampling may include conversion of a continuous signal to a discrete signal.

In many fields, including simulations of physical systems, and computing, especially analog computing, the foregoing meanings may merge. For example, a hybrid computer can draw samples from an analog computer. The analog computer, as a provider of samples, is an example of a sample generator. The analog computer can be operated to provide samples from a selected probability distribution, the probability distribution assigning a respective probability of being sampled to each data point in the population.

An analog processor, for example a quantum processor and in particular a quantum processor designed to perform quantum annealing and/or adiabatic quantum computation, may be operated as a sample generator. The population can correspond to all possible states of the processor, and each sample can correspond to a respective state of the processor. Using an analog processor as a sample generator may be a preferred mode of operating the processor for certain applications. Operating an analog processor as a sample generator may also enable a broader range of problems to be solved compared to, for example, using an analog processor to find a low energy state of a Hamiltonian that encodes an optimization problem.

Importance Sampling

Importance sampling is a technique for estimating properties of a distribution of interest, by drawing samples from a different distribution, and weighting the samples as needed to recover the distribution of interest. When combined with the normalization constants of both distributions, the resulting estimator can be unbiased. Alternatively, one can normalize the estimator by the sum of the weights to yield an asymptotically unbiased estimator.

Consider the expectation values of a function h(x) over a distribution of interest, π(x), also called the target distribution:

I≡∫h(x)π(x)dx (1)

Another distribution, the proposal distribution q(x), can be introduced into equation (1), provided q(x) is strictly positive with respect to π(x). In this case, the integral becomes:

I
=

∫

h


(
x
)



π


(
x
)

q


(
x
)



q


(
x
)


d

x

(
2
)

The factor w(x)=π(x)/q(x) can be called the importance weight. Since each of the distributions π(x) and q(x) may only be known to within a respective constant factor, a new weight can be defined as follows {tilde over (w)}(x)≡{tilde over (π)}(x)/{tilde over (q)}(x) where the tilde {tilde over ( )} denotes the absence of normalization for the distributions.

The integral in equation (2) can be approximated by a discrete sum over samples from the distribution. For N samples, an unbiased estimate of the integral is given by:

I
^

=

∑

i
=
1

N



h


(

X

(
i
)

)



w
~



(

X

(
i
)

)

∑

i
=
1

N



w
~



(

X

(
i
)

)

(
3
)

A challenge in importance sampling can be finding a suitable proposal distribution q(x). A poor choice of proposal can result in a large variance, which means that an undesirably large number of samples must be drawn from the proposal distribution q(x) before the weighted set is representative of the target distribution π(x). As the number of dimensions in the distribution increases, the problem can become more pronounced.

Markov Chain Monte Carlo

Markov Chain Monte Carlo (MCMC) is a class of computational techniques which include, for example, simulated annealing, parallel tempering, population annealing, and other techniques. A Markov chain may be used, for example when a probability distribution cannot be used. A Markov chain may be described as a sequence of discrete random variables, and/or as a random process where at each time step the state only depends on the previous state. When the chain is long enough, aggregate properties of the chain, such as the mean, can match aggregate properties of a target distribution.

The Markov chain can be obtained by proposing a new point according to a Markovian proposal process (generally referred to as an “update operation”). The new point is either accepted or rejected. If the new point is rejected, then a new proposal is made, and so on. New points that are accepted are ones that make for a probabilistic convergence to the target distribution. Convergence is guaranteed if the proposal and acceptance criteria satisfy detailed balance conditions and the proposal satisfies the ergodicity requirement. Further, the acceptance of a proposal can be done such that the Markov chain is reversible, i.e., the product of transition rates over a closed loop of states in the chain is the same in either direction. A reversible Markov chain is also referred to as having detailed balance. Typically, in many cases, the new point is local to the previous point.

Superconducting Qubits

Superconducting qubits are solid state qubits based on circuits of superconducting materials. Operation of superconducting qubits is based on the underlying principles of magnetic flux quantization, and Josephson tunneling. Superconducting effects can be present in different configurations, and can give rise to different types of superconducting qubits including flux, phase, charge, and hybrid qubits. The different configurations can vary in the topology of the loops, the placement of the Josephson junctions, and the physical parameters of elements of the superconducting circuits, such as inductance, capacitance, and Josephson junction critical current.

Quantum Processor

A quantum processor may take the form of a superconducting quantum processor. A superconducting quantum processor may include a number of superconducting qubits and associated local bias devices. A superconducting quantum processor may also include coupling devices (also known as couplers) that selectively provide communicative coupling between qubits.

In one implementation, the superconducting qubit includes a superconducting loop interrupted by a Josephson junction. The ratio of the inductance of the Josephson junction to the geometric inductance of the superconducting loop can be expressed as 2πL/c/Φ₀(where L is the geometric inductance, Ic is the critical current of the Josephson junction, and Φ₀is the flux quantum). The inductance and the critical current can be selected, adjusted, or tuned, to increase the ratio of the inductance of the Josephson junction to the geometric inductance of the superconducting loop, and to cause the qubit to be operable as a bistable device. In some implementations, the ratio of the inductance of the Josephson junction to the geometric inductance of the superconducting loop of a qubit is approximately equal to three.

In one implementation, the superconducting coupler includes a superconducting loop interrupted by a Josephson junction. The inductance and the critical current can be selected, adjusted, or tuned, to decrease the ratio of the inductance of the Josephson junction to the geometric inductance of the superconducting loop, and to cause the coupler to be operable as a monostable device. In some implementations, the ratio of the inductance of the Josephson junction to the geometric inductance of the superconducting loop of a coupler is approximately equal to, or less than, one.

Further details and embodiments of exemplary quantum processors that may be used in conjunction with the present systems and devices are described in, for example, U.S. Pat. Nos. 7,533,068; 8,008,942; 8,195,596; 8,190,548; and 8,421,053.

Adiabatic Quantum Computation

A Hamiltonian is an operator whose eigenvalues are the allowed energies of the system. Adiabatic quantum computation can include evolving a system from an initial Hamiltonian to a final Hamiltonian by a gradual change. One example of adiabatic evolution is a linear interpolation between the initial Hamiltonian H_iand the final Hamiltonian H_f, as follows:

H_e=(1−s)H_i+sH_f (4)

where H_eis the evolution, or instantaneous, Hamiltonian, and s is an evolution coefficient that can control the rate of evolution.

As the system evolves, the evolution coefficient s changes value from 0 to 1. At the start, the evolution Hamiltonian H_eis equal to the initial Hamiltonian H_i, and, at the end, the evolution Hamiltonian H_eis equal to the final Hamiltonian H_f.

The system is typically initialized in a ground state of the initial Hamiltonian H_i, and the goal of the adiabatic evolution is to evolve the system such that it ends up in a ground state of the final Hamiltonian H_fat the end of the evolution. If the evolution is too fast, then the system can transition to a higher energy state of the system, such as the first excited state.

The process of changing the Hamiltonian in adiabatic quantum computing may be referred to as evolution. An adiabatic evolution is an evolution that satisfies an adiabatic condition such as:

{dot over (s)}| [Image Omitted] 1|dH_e/ds|0|=δg²(s) (5)

where {dot over (s)} is the time derivative of s, g(s) is the difference in energy between the ground state and first excited state of the system (also referred to herein as the gap size) as a function of s, and δ is a coefficient and δ<<1.

If the rate of change (for example, {dot over (s)}), is slow enough that the system is always in the instantaneous ground state of the evolution Hamiltonian, then transitions at anti-crossings (i.e., when the gap size is smallest) are avoided. Equation (4) above is an example of a linear evolution schedule. Other evolution schedules can be used, including non-linear, parametric, and the like. Further details on adiabatic quantum computing systems, methods, and apparatus are described in, for example, U.S. Pat. Nos. 7,135,701; and 7,418,283.

Quantum Annealing

Quantum annealing is a computation method that may be used to find a low-energy state of a system, for example the ground state of the system. Similar in concept to classical simulated annealing, the method relies on an underlying principle that natural systems tend towards lower energy states because lower energy states are more stable. While classical annealing uses classical thermal fluctuations to guide a system to a low-energy state and ideally its global energy minimum, quantum annealing may use quantum effects, such as quantum tunneling, as a source of disordering to reach a global energy minimum more accurately and/or more quickly than classical annealing. In quantum annealing thermal effects and other noise may be present, and the final low-energy state may not be the global energy minimum.

Adiabatic quantum computation may be considered a special case of quantum annealing for which the system, in an ideal implementation, begins and remains in its ground state throughout an adiabatic evolution. Thus, those of skill in the art will appreciate that quantum annealing systems and methods may generally be implemented on an adiabatic quantum computer. Throughout this specification and the appended claims, any reference to quantum annealing is intended to encompass adiabatic quantum computation unless the context requires otherwise.

Quantum annealing can use quantum mechanics as a source of disorder during the annealing process. An objective function, such as an optimization problem, can be encoded in a problem Hamiltonian H_p, and the method can introduce quantum effects by adding a disordering Hamiltonian H_Dthat does not commute with the problem Hamiltonian H_p. An example case is as follows:

H_E=A(t)H_D+B(t)H_p (6)

where A(t) and B(t) are time-dependent envelope functions. For example, A(t) can be a function that changes from a large initial value to substantially zero during the evolution. In this case, H_Emay be thought of as an evolution Hamiltonian similar to H_edescribed above in the context of adiabatic quantum computation.

The disorder is slowly removed by reducing the impact of H_Dthrough a reduction in the value of A(t) over time. In this respect, quantum annealing can be similar to adiabatic quantum computation in that the system starts with an initial Hamiltonian, and evolves through an evolution Hamiltonian to a final problem Hamiltonian H_pwhose ground state encodes a solution to the problem.

If the evolution is slow enough, the system may settle in a global minimum, or in a local minimum close in energy to the global minimum. The performance of the computation may be assessed via a residual energy as a function of evolution time. The computation time can be defined as the time required for the residual energy to fall below a predetermined threshold.

In quantum annealing, H_pmay encode an optimization problem and therefore H_pmay be diagonal in the subspace of the qubits that encode the solution, but the system does not necessarily stay in the ground state at all times. The energy landscape of H_pmay be crafted so that its global minimum is the answer to the problem to be solved, and low-lying local minima are good approximations.

The gradual reduction of disordering Hamiltonian H_Din quantum annealing may follow a defined schedule known as an annealing schedule. Unlike adiabatic quantum computation where the system begins and remains in its ground state throughout the evolution, in quantum annealing the system may not remain in its ground state throughout the entire annealing schedule. As such, quantum annealing may be implemented as a heuristic technique, where low-energy states with energy near that of the ground state may provide approximate solutions to the problem. The method can include adding a disordering Hamiltonian H_Dand subsequently removing it as described above.

The techniques and structures described herein generally relate to sampling from an available probability distribution to create a desirable probability distribution. This resultant distribution can be used for computing values used in computational techniques including: Importance Sampling and Markov chain Monte Carlo systems.

BRIEF SUMMARY

A method of operation of a hybrid computing system to generate approximate equilibrium statistics, the hybrid computing system comprising an analog processor and a digital processor, can be summarized as including: determining an intermediate value of an annealing parameter, the annealing parameter normalized between a value of zero and a value of one; determining a ramp time t_rampwherein t_rampis smaller than the relaxation time of the analog processor and greater than a time needed for non-adiabatic transition to higher energy states; initializing the analog processor in a first state; evolving the analog processor from a value of one for the annealing parameter until the annealing parameter is equal to the intermediate value, wherein evolving the analog processor from a value of one for the annealing parameter until the annealing parameter is equal to the intermediate value occurs in the ramp time; holding the annealing parameter of the analog processor to the intermediate value for a predetermined dwell time; evolving the analog processor from the intermediate value of the annealing parameter until the annealing parameter is equal to one, wherein evolving the analog processor from the intermediate value of the annealing parameter until the annealing parameter is equal to one occurs in the ramp time and after the ramp time the analog processor is in a second state; and reading out the second state of the analog processor.

Determining an intermediate value of an annealing parameter may include determining a value of the annealing parameter where the evolution of the analog processor starts deviating toward a distribution that is not a Boltzmann distribution.

Initializing the analog processor in a first state may include initializing the state of the analog processor via the digital processor to correspond to a distribution of interest.

Initializing, the analog processor in a first state may include evolving the analog processor from an end state of a previous annealing operation.

The dwell time may be a time greater than a time needed for the evolution of the analog processor to recover from any distortion caused by evolving the analog processor by the ramp time.

Reading out the second state of the analog processor may include reading out the second state of the analog processor with the digital processor.

The second state may not be equal the first state.

The second state may be a sample from a distribution obtained from an evolution in t_rampfrom a thermal equilibrium distribution

P
β



(
z
)

=

e

-
β



E


(
z
)

z
β

.

The method may further include: after reading out the second state of the analog processor, post-processing the second state via the digital processor to obtain a third state of the analog processor.

Post-processing the second state of the quantum processor may include applying one or more passes of quantum Monte Carlo simulations to the second state.

The method may further include: assigning importance weights to the third state of the analog processor via the digital processor.

A hybrid computing system to generate approximate equilibrium statistics, the hybrid computing system comprising an analog processor and a digital processor, may execute any of the acts described above.

A method of operation of a hybrid computing system to generate approximate equilibrium statistics, the hybrid computing system comprising an analog processor and a digital processor, may be summarized as including: initializing a counter i to i=0; determining an intermediate value s_i* of an annealing parameter, the annealing parameter normalized between a value of zero and a value of one; determining a ramp time t_rampwherein t_rampis smaller than the relaxation time of the analog processor and greater than the time needed for non-adiabatic transition to higher energy states; initializing, the analog processor to a state i; Iterating until an exit condition is met: evolving the analog processor from a value of one for the annealing parameter until the annealing parameter is equal to s_i*, wherein evolving the analog processor from a value of one for the annealing parameter until the annealing parameter is equal to the intermediate value occurs in the ramp time; holding the annealing parameter of the analog processor to s_i*for a defined dwell time; evolving the analog processor from s_i* of the annealing parameter until the annealing parameter is equal to one, wherein evolving the analog processor from the intermediate value of the annealing parameter until the annealing parameter is equal to one occurs in the ramp time and after the ramp time the analog processor is in a state i+1; and reading out the state i+1 of the analog processor; if the exit condition is not met, determining a new value s_i+1*>s_i* of the annealing parameter; and updating the counter i to i=i+1.

Determining an intermediate value of an annealing parameter, the annealing parameter normalized between a value of zero and a value of one may include determining a value of the annealing parameter where the evolution of the analog processor starts deviating toward a distribution that is not a Boltzmann distribution.

The dwell time may be a time greater than the time needed for the evolution of the analog processor to recover from any distortion caused by evolving the analog processor by the ramp time.

Reading out the state i=1 of the analog processor may include reading out the state i+1 of the analog processor with the digital processor.

The state i+1 may not be equal the state i.

The state i+1 may be a sample from a distribution obtained from an evolution in t_rampfrom a thermal equilibrium distribution

P
β



(
z
)

=

e

-
β



E


(
z
)

Z
β

.

The method may further include: after reading out the state i+1 of the analog processor, post-processing the state i+1 via the digital processor to obtain a post-processed state i+1 of the analog processor.

Post-processing the state i+1 of the quantum processor may include applying one or more passes of quantum Monte Carlo simulations to the state i+1.

The method may further include: assigning importance weights to the post-processed state i+1 of the analog processor via the digital processor; and initializing, by the digital processor, the analog processor to the post-processed state i+1.

The method may further include: after reading out the state i+1 of the analog processor storing the state i+1 via the digital processor; and after meeting the exit condition post processing the states of the analog processor via the digital processor to obtain importance weights.

Post processing the state of the analog processor via the digital processor may include performing one or more passes of quantum Monte Carlo post processing.

A hybrid computing system to generate approximate equilibrium statistics, the hybrid computing system comprising an analog processor and a digital processor, may execute any of the acts described above.

BRIEF DESCRIPTION OF THE SEVERAL VIEWS OF THE DRAWING(S)

In the drawings, identical reference numbers identify similar elements or acts. The sizes and relative positions of elements in the drawings are not necessarily drawn to scale. For example, the shapes of various elements and angles are not necessarily drawn to scale, and some of these elements may be arbitrarily enlarged and positioned to improve drawing legibility. Further, the particular shapes of the elements as drawn, are not necessarily intended to convey any information regarding the actual shape of the particular elements, and may have been solely selected for ease of recognition in the drawings.

FIG. 1 is a schematic diagram of an exemplary hybrid computing system including a digital processor and an analog processor useful in performing the methods described in the present disclosure.

FIG. 2 is a flow diagram of an exemplary method of operation of a hybrid computing system to generate statistics.

FIG. 3 is a flow diagram of an exemplary method of operation of a hybrid computing system to generate equilibrium statistics via post processing and calculating importance weights.

FIG. 4 is a flow diagram of an exemplary method of operation of a hybrid computing system to generate equilibrium statistics and calculating importance weights by sequential iteration of the analog processor.

FIG. 5A is a chart showing an exemplary set of evolutions performed during an example implementation of the method of FIG. 2.

FIG. 5B is a chart showing an exemplary set of evolutions performed during an example implementation of the method of FIG. 6.

FIG. 6 is a flow diagram of an exemplary method of operation of a hybrid computing system to implement an update operation for a classical MCMC method using re-equilibration.

DETAILED DESCRIPTION

In the following description, some specific details are included to provide a thorough understanding of various disclosed embodiments. One skilled in the relevant art, however, will recognize that embodiments may be practiced without one or more of these specific details, or with other methods, components, materials, etc. In other instances, well-known structures associated with quantum processors, such as quantum devices, coupling devices, and control systems including microprocessors and drive circuitry have not been shown or described in detail to avoid unnecessarily obscuring descriptions of the embodiments of the present methods. Throughout this specification and the appended claims, the words “element” and “elements” are used to encompass, but are not limited to, all such structures, systems, and devices associated with quantum processors, as well as their related programmable parameters.

Unless the context requires otherwise, throughout the specification and claims which follow, the word “comprise” and variations thereof, such as, “comprises” and “comprising” are to be construed in an open, inclusive sense, that is as “including, but not limited to.”

Reference throughout this specification to “one embodiment” “an embodiment”, “another embodiment”, “one example”, “an example”, or “another example” means that a particular referent feature, structure, or characteristic described in connection with the embodiment or example is included in at least one embodiment or example. Thus, the appearances of the phrases “in one embodiment”, “in an embodiment”, “another embodiment” or the like in various places throughout this specification are not necessarily all referring to the same embodiment or example. Furthermore, the particular features, structures, or characteristics may be combined in any suitable manner in one or more embodiments or examples.

It should be noted that, as used in this specification and the appended claims, the singular forms “a,” “an,” and “the” include plural referents unless the content clearly dictates otherwise. Thus, for example, reference to a problem-solving system including “a quantum processor” includes a single quantum processor, or two or more quantum processors. It should also be noted that the term “or” is generally employed in its sense including “and/or” unless the content clearly dictates otherwise.

The headings provided herein are for convenience only and do not interpret the scope or meaning of the embodiments.

Overview

Sampling from a Boltzmann distribution of a classical or quantum Hamiltonian can be challenging. The statistical nature of an open physical quantum annealer, such as a quantum annealing processor or an adiabatic quantum computer, can be exploited for sampling purposes.

Typically, it is a challenge for a quantum annealer to be able to equilibrate quickly enough at points during the annealing process to provide equilibrium samples efficiently. One approach to mitigating this challenge is to combine the use of a quantum processor performing quantum annealing with the use of a classical or digital processor performing classical methods.

One way to combine the use of a quantum processor and a classical or digital processor is to use post-processing techniques (such as described in International PCT Application WO2016029172A1 “Systems and Methods for Problem Solving, Useful For Example In Quantum Computing”). Typically, the amount of post-processing needed to achieve a desired accuracy can become large enough that it at least partially negates the benefit of obtaining initial samples from the quantum processor.

The present systems, devices, articles, and methods relate to sampling from an available probability distribution and making use of the samples. Sampling can be performed via a digital computer on its own, or performed via a digital computer in cooperation with an analog computer, such as, a quantum computer. FIG. 1 shows a digital computer communicatively coupled to a quantum computer. Use of a quantum computer for sampling may advantageously allow samples to be taken from many disparate low energy states.

FIG. 1 illustrates a hybrid computing system 100 including a digital computer 102 coupled to an analog computer 104. In some implementations the analog computer 104 is a quantum computer. The exemplary digital computer 102 includes a digital processor 106 that may be used to perform classical digital processing tasks.

Digital computer 102 may include at least one digital processor 106 with one or more cores, at least one system memory 108, and at least one system bus 110 that couples various system components, including system memory 108 to digital processor 106.

The digital processor 106 may be any circuitry that forms a logic processing unit, such as one or more central processing units (“CPUs”), graphics processing units (“GPUs”), digital signal processors (“DSPs”), application-specific integrated circuits (“ASICs”), programmable gate arrays (“FPGAs”), programmable logic controllers (PLCs), etc.

Digital computer 102 may include a user input/output subsystem 112. In some implementations, the user input/output subsystem includes one or more user input/output components such as a display 114, mouse 116, and/or keyboard 118.

System bus 110 can employ any known bus structures or architectures, including a memory bus with a memory controller, a peripheral bus, and a local bus. System memory 108 may include non-volatile memory, such as read-only memory (“ROM”), static random access memory (“SRAM”), Flash NAND; and volatile memory such as random access memory (“RAM”) (not shown).

Digital computer 102 may also include other non-transitory computer- or processor-readable storage media or non-volatile memory 120. Non-volatile memory 120 may take a variety of forms, including: spinning media for instance a hard disk drive (HDD) for reading from and writing to a magnetic hard disk and/or an optical disk drive for reading from and writing to removable optical disks, and/or non-spinning media for instance a solid state drive (SSD) for reading from and writing to solid state memory. The optical disk can be a CD-ROM or DVD, while the magnetic disk can be a magnetic floppy disk or diskette or one or more magnetic disc platters. Non-volatile memory 120 may communicate with digital processor via system bus 110 and may include appropriate interfaces or controllers 122 coupled to system bus 110. Non-volatile memory 120 may serve as long-term storage for processor- or computer-readable instructions, data structures, or other data (sometimes called program modules) for digital computer 102.

Although digital computer 102 has been described as employing hard disks, optical disks and/or magnetic disks, those skilled in the relevant art will appreciate that other types of non-volatile computer-readable media may be employed, such magnetic cassettes, flash memory cards, Flash, ROMs, smart cards, etc. Those skilled in the relevant art will appreciate that some computer architectures employ volatile memory and non-volatile memory. For example, data in volatile memory can be cached to non-volatile memory or in a solid-state drive that employs integrated circuits to provide non-volatile memory.

Various processor- or computer-readable instructions, data structures, or other data can be stored in system memory 108. For example, system memory 108 may store instruction for communicating with remote clients and scheduling use of resources including resources on the digital computer 102 and analog computer 104.

In some implementations system memory 108 may store processor- or computer-readable calculation instructions to perform pre-processing, co-processing, and post-processing to analog computer 104. System memory 108 may store at set of analog computer interface instructions to interact with the analog computer 104.

Analog computer 104 may include an analog processor such as quantum processor 124. The analog computer 104 can be provided in an isolated environment, for example, in an isolated environment that shields the internal elements of the quantum computer from heat, magnetic field, and other external noise (not shown) and/or which cools the analog processor to temperatures (i.e., critical temperature) at or below which the circuitry of the analog processor becomes superconductive. In contrast, the digital computer 102 will typically operate at much higher temperatures (e.g., room temperature) at which superconductivity does not occur and/or may employ materials that do not superconduct even at or below the critical temperature.

The systems and methods described in this document can generate weighted samples from a quantum Boltzmann distribution with a moderate amount of post-processing. The methods comprise generating samples using a quantum annealing processor and performing annealed importance sampling, and is referred to herein as quantum importance sampling.

Sampling from a Classical Distribution

A thermal equilibrium probability distribution can be expressed as follows:

P
β



(
z
)

=

e

-
β



E


(
z
)

Z
β

(
7
)

where β=1/k_BT, k_Bis the Boltzmann constant, T is temperature, E(z) is the classical energy corresponding to a state z∈{−1, +1}^N, and Z_β=Σ_ze^−βE(z)is the partition function.

This distribution can be reached by stochastic processes such as Metropolis algorithm or Markov chain Monte Carlo (MCMC), which can be represented by a transition matrix P_β(z|z′). A very large number of applications of such algorithms, starting from any random state, might be able to give P_β(z) with an acceptable degree of accuracy. For that to hold, P_β(z) should be a stationary state of such a process, or an Eigen function of P_β(z|z′) with eigenvalue 1:

P
β



(
z
)

=

∑

z
′



P
β



(

z
|

z
′

)



P
β



(

z
′

)

(
8
)

The detailed balance condition, i.e. P_β(z|z′)P_β(z′)=P_β(z′|z)P_β(z), would imply give equation (8); however, the opposite is not true. Equations (7) and (8) allow to describe the equilibrium expectation value of an observable 0 at β in terms of another distribution P_β′(z) at a different temperature β′:

〈
O
〉

β

=

∑

z
,

z
′



O


(
z
)



w


(

z
′

)



P
β



(

z
|

z
′

)



P

β
′



(

z
′

)




where

(
9
)

w


(
z
)

=

P
β



(
z
)

P

β
′



(
z
)

=

Z

β
′

Z
β



e

-

(

β
-

β
′

)



E


(
z
)

(
10
)

Is a weight assigned to the state z′. If the new temperature is larger (β′>β), then P_β, is an easier distribution to construct. It is therefore possible to use a simpler distribution to generate statistics from a harder distribution.

The above sum (9) can be approximated by weighted sampling as follows:

[Image Omitted] O_β=_z,z′,[O(z)w(z′)] (11)

Where [Image Omitted] _z,z′ is an expectation value estimated using sample z and z′.

Sampling can be performed as follows: for the k^thsample, generate the first state z′ by sampling from the distribution P_β′, generate the second state z using a MCMC at temperature β starting from state z′, and weight the sample by:

w
k

=

Z

β
′

Z
β



e

-

(

β
-

β
′

)



E


(
z
)

(
12
)

The process can be repeated many times to generate a large number of samples, each weighted by a w_k. The pre-factor is the same in each case, so it can be eliminated by normalizing the weights.

Typically, in simulated annealing, the process starts at a high temperature, and gradually reduces the temperature, evolving the state using an MCMC, without introducing weights. If the goal is to obtain a true Boltzmann distribution at a low temperature, then the process should keep the system in equilibrium at intermediate temperatures. To do this, the process can apply many MCMC after each temperature change, and the process can be time-consuming.

By introducing weights, simulated annealing can anneal the system more quickly while producing the desired equilibrium statistics. The trade-off is that some samples can be lost because of the exponentially small weights accumulated during sampling.

Sampling from a Quantum Boltzmann Distribution

In one implementation of the present systems and methods, annealed importance sampling is used with quantum Monte Carlo simulation to generate samples from a quantum Boltzmann distribution, i.e., a Boltzmann distribution of a quantum Hamiltonian.

Annealing can be performed by either changing the temperature or other Hamiltonian parameters such as the transverse field. The notation θ may be used to represent all such parameters, including temperature.

In quantum Monte Carlo (QMC) simulation, the quantum statistics are reproduced by Monte Carlo sampling from trajectories (paths) x(r) along the imaginary time dimension τ∈[0,β]. Each trajectory can appear in the sampling with probability

P
θ



[

x


(
τ
)

]

=

e

-

S
θ



[

x


(
τ
)

]

Z
θ

(
13
)

Monte Carlo updates change the trajectory x→x′, which can be represented by a transition matrix P_θ[x(τ)|x′(τ)] that satisfied the stationary condition

P_θ[x(τ)]=∫ [Image Omitted] [x′(τ)]P_θ[x(τ)|x′(τ)]P_θ[x′(τ)] (14)

Equation (13) and (14) can replace equations (7) and (8) as the basis for importance sampling. Equilibrium statistics can therefore the generated corresponding to a Hamiltonian H at the point θ using weighted equilibrium samples obtained at another point θ′.

Quantum Annealing

Quantum annealing can be implemented in a physical quantum annealer such as an adiabatic quantum computer. In some implementations, the physical quantum annealer implements the following Hamiltonian:

H


(
s
)

=

-

A


(
s
)



∑
i



σ
i
x

+

B


(
s
)

[

∑
i



h
i



σ
i
z

+

∑

i
,
j
,



J

i
,
j



σ
i
z



σ
j
z

]

(
15
)

where h_iand J_i,jare tuneable or programmable parameters, and A(s) and B(s) are monotonic functions of an annealing parameter s, where 0≤s≤1, A(0)>>B(0)≈0 and B(1)>>A(1)≈0.

In some implementations, the annealing parameters is varied linearly with time t: s(t)=t/t_α, where t_α is the annealing time. Alternatively, or in addition, the annealing parameter s may be varied non-linearly with time.

For an infinitely slow evolution, a quantum annealer will follow equilibrium distribution of its instantaneous Hamiltonian H(s). However, in many applications, a system comprising a quantum annealer will evolve close to equilibrium up to some point and then gradually deviate from equilibrium distribution. This can present challenges in certain applications. For example, if one wishes to sample from a Boltzmann distribution, this de-equilibration can lead to a final distribution that is not exactly a Boltzmann distribution. A Boltzmann distribution may be obtained, for example, by reading out the state of a quantum annealer at an intermediate point during the anneal, before the deviation starts.

Re-Equilibrated Quantum Sampling

In some implementations, the physical quantum annealer performs a re-equilibration act. This may be in response to an act which pushes the state of the annealer away from equilibrium, and/or as a countermeasure to the gradual deviation from equilibrium described above.

In some implementations, the physical quantum annealer may de-equilibrate by executing a fast ramp operation (a steep linear variation of the annealing parameters) from an intermediate point part way through annealing to the end of the annealing:

s


(
t
)

=

s
*

+

(

t
-

t
1

)

×

(

1
-

s
*

)

(

t
a

-

t
1

)

(
16
)

where t₁is the start of the ramp, and s* is the annealing parameter at an intermediate point during the anneal.

Examples of fast ramp operations, including reverse annealing to an intermediate value of s, are found in International Patent Application No. US2016/059,169 and U.S. Provisional Application No. 62/417,940.

The ramp can be selected to be faster than the relaxations of the system, while not so fast as to cause non-adiabatic transitions to higher energy states. At least some distortion caused by non-adiabaticity can be recovered by a sufficiently quick thermalization step, provided the ramp is slow compared to the time scale determined by the temperature. A suitable choice of ramp time t_rampis given by:

β≲t_ramp≲t_relax (17)

where t_relaxis the shortest relaxation time at intermediate point s*.

The distribution obtained with a fast ramp operation from the equilibrium distribution P_θ(z) at point θ can be defined as:

{tilde over (P)}_θ({tilde over (z)})=P_θ(z)+fast ramp (18)

As mentioned before, in general {tilde over (P)}_θ≠P_θ because of the non-equilibrium evolution during the fast ramp operation. In some cases, where t_rampis smaller than all the relaxation times at point s*, the difference between the distributions may be negligible and {tilde over (P)}_θ may be used to represent an equilibrium distribution with some confidence. In other cases, however t_rampmay be larger than some local relaxations, for example relaxations that flip only one or two qubits in the quantum processor. In such case, it may be possible to recover the equilibrium distribution with post-processing.

A method for using a quantum processor with fast ramp may be summarized by a transition matrix P_θ^QA({tilde over (z)},{tilde over (z)}′) as follows:

P
θ

Q

A



(

z
˜

,

z
˜

′

)

=

{

1.



initialize



qubits



in



state




z
˜

′

2.



reverse



anneal



to



s
*


within




t
ramp

3.



hold



at



s
*
for



time




t
hold

4.



ramp



to



s

=

1



within




t
ramp

5.



measure



the



state




z
˜

(
19
)

Equation (19) will be described in FIG. 2.

FIG. 2 is a flow chart illustrating a method 200 for operating a hybrid computing system including at least one digital processor and at least one analog processor with a fast ramp operation and a reverse annealing operation for generating equilibrium statistics. The analog processor may be a quantum annealing processor.

Method 200 comprises acts 202 through 226. For method 200, as with others methods taught herein, the various acts may be performed in a different order than that illustrated and described. Additionally, the methods can omit some acts, and/or employ additional acts. One or more of these acts may be performed by or via one or more circuits, for instance one or more processors (e.g., digital processors such as microprocessors, analog processor such as quantum processors, a hybrid computer including a digital process and an analog processor). A sampling device, such as, a device including a hybrid computer, may perform one or more of the acts in method 200.

Method 200 starts at 202, for example in response to another routine.

At 204 the at least one digital processor initializes a counter i to a starting value, e.g., i=0.

At 206, the at least one digital processor determines an intermediate value s* of the annealing parameter s. The intermediate value s_i* may be supplied to the at least one digital processor as part of a set of input or determined by the at least one digital processor.

At 208, the at least one digital processor determines a time interval t_ramp. The value of t_rampmay be supplied to the at least one digital processor as part of a set of input or determined by the at least one digital processor. t_rampis chosen so as to satisfy equation (17).

At 210, the at least digital processor programs the at least one analog processor so that at least a portion of the qubits in the at least one analog processor is in an initial state z at an initial time corresponding to counter i. The initial state z may be an equilibrium state corresponding to a distribution P_θ(z). In one implementation, distribution P_θ(z) is a Boltzmann distribution. The at least one analog processor may be programmed so that after an annealing operation the samples obtained corresponds to samples from distribution P_θ(z).

From 212 to 222, method 200 iterates over linear increments of counter i.

At 212, the at least one digital processor programs the at least one analog processor so that the at least analog processor evolves from a value s=1 to a value s_i* of the annealing parameter s within a time t_ramp, where s_i* and t_rampwhere determined at 206 and 208, respectively. FIG. 5A shows an exemplary graph 500A illustrating this evolution as line 504 occurring over time interval t_ramp. (Axes 522 and 524 correspond to time t and annealing parameter s, respectively).

At 214, the at least one digital processor clamps the state of the at least one analog processor for a time t_hold. As mentioned before, the time t_holdshould be long enough for the system to re-equilibrate and recover from any distortion that may have caused the state of the at least one analog processor to deviate from equilibrium. A distortion may be caused by the fast ramp operation at 212. The time t_holdmay be determined by the at least digital processor at 214 or supplied to the at least one digital processor as part of a set of inputs. An example of this clamping is shown as line 506 occurring over time interval t_holdin FIG. 5A.

At 216, the at least one digital processor programs the at least one analog processor so that the at least one analog processor evolves from s_i* to s=1 in time t_ramp. The time t_rampmay be the same time interval applied at 212 or different (e.g., a longer time interval), but still satisfying condition (17). An example of this evolution is shown as line 508 occurring over time interval t_rampin FIG. 5A.

At 218, while the at least analog processor is in a classical state, the at least one digital processor reads out the state of the at least one analog processor. The state read out at 218 is sampled from the same distribution P_θ(z).

Quantum Importance Sampling

Annealed importance sampling (AIS) cannot be easily directly applied to a quantum annealing system because in order to calculate the weights, the energy of the system needs to be measured, but in a quantum annealing system measurements are commonly done in a computation basis, not an energy basis. In addition, a measurement of the state of a quantum annealer, when the quantum annealer is in a quantum state, will destroy the quantum state of the system, so that the state read is a classical state, which is not desirable in annealed importance sampling.

The below method can be applied with a quantum processor and quantum Monte Carlo simulation for obtaining weighted samples from an equilibrium distribution.

Optionally, at 220, the at least one digital processor determines whether an exit condition has been met. In one implementation the exit condition may be the completion of a predetermined numbers of iterations, while in other implementations, the exit conditions may be the measure of a performance factor reaching a pre-determined threshold. If the exit condition is not met at 220, control passes to 222, otherwise to 226. In some implementations, the method proceeds to end at 226 without performing a determination of an exit condition.

At 222, the at least one digital processor determines a new intermediate value s*₁₊₁for the annealing parameter s, where s*_i+1>s*_i.

At 224, the at least one digital processor increases the value of the counter i to i=i+1, before passing control to 212.

At 226, method 200 terminates, until for example it is invoked again. At each iteration of method 200 the at least one analog processor starts a new evolution cycle from an equilibrium state determined at the previous iteration. That is, when the at least one analog processor recovers from any distortion, if the initial state {tilde over (z)}′ is sampled from {tilde over (P)}_θ({tilde over (z)}′), the state read out at 218 will also be sampled from the same distribution, so that

P
˜

θ



(

z
˜

)

=

∑

z
˜

′



P
θ

Q

A



(

z
˜

|

z
˜

′

)



P
˜

θ



(

z
˜

′

)

(
20
)

Equation (20) may be used instead of (8) for importance sampling. However, {tilde over (P)}_θ({tilde over (z)}) has no closed form that can be used as an equivalent of equation (7), preventing a direct application of annealed importance sampling via method 200.

Quantum Importance Sampling with Quantum Monte Carlo Post-Processed Samples

The distribution obtained after a fast ramp from an equilibrium distribution may deviate locally from the desired Boltzmann distribution. In one implementation, the system uses quantum Monte Carlo (QMC) post-processing to recover the desired distribution.

As mentioned before, for small enough t_rampdeviation of {tilde over (P)}_θ from P_θ is local, i.e. it involves only short hamming distance relaxations. Therefore, it may be possible to recover, or approximate with a degree of certainty, P_θ with QMC post-processing.

After an annealing operation, the read-out state {tilde over (z)} is post-processed to obtain a state z, for example with quantum Monte Carlo post-processing. The post-processing operation can be defined by the probability distribution P_θ^QMC(z|{tilde over (z)}), which is the probability distribution of a state z obtained after a few passes of QMC simulation with parameter θ, staring from state {tilde over (z)}.

P_θ^QA({tilde over (z)},{tilde over (z)}′) as defined in (19) can be used with P_σ^QMC(z|{tilde over (z)}) as a hybrid operation P_θ^HYBwhich take an initial state z′ to generate {tilde over (z)}, then runs a few iterations of quantum Monte Carlo simulation at parameter θ to generate state z:

P
θ

H

Y

B



(

z
|

z
′

)

=

∑

z
˜



P
θ

Q

M

C



(

z
|

z
˜

)



P
θ

Q

A



(

z
˜

|

z
′

)

(
21
)

Post-processing, for example QMC post-processing, may restore equilibrium from local deviations due to a fast ramp operation. Thus, starting from an equilibrium state it is possible to get back to an equilibrium state:

P
θ



(
z
)

=

∑
z



P
θ

H

Y

B



(

z
|

z
′

)



P
θ



(

z
′

)

(
22
)

Similar to (10), importance weights can be defined as

w


(
z
)

=

P
θ



(
z
)

P

θ
′



(
z
)

(
23
)

To estimate weights (23) it is possible to define a dimensionless local free-energy [Image Omitted] _θ(z) as a quantum analog to the local (free) energy to be used to calculate the weights:

[Image Omitted] (24)

Where [Image Omitted] is a normalization factor, and only local equilibrations of QMC are needed to calculate _θ.

The above method can be generalized by taking multiple samples z₁, . . . z_nwith a quantum annealing processor and a digital processor, where each sample, z_i, is obtained from z_i−1using P_θ^HYB(z|z′) as defined in (21), then calculating normalized importance weights:

w(z₁, . . . ,z_n-1)

[Image Omitted] (25)

FIG. 3 is a flow diagram of an exemplary method 300 of operation of a hybrid computing system comprising at least one digital processor and at least one analog processor to generate equilibrium statistics via post processing and calculating importance weights. The analog processor may be a quantum annealing processor.

Method 300 comprises acts 302 through 332. For method 300, as with others methods taught herein, the various acts may be performed in a different order than that illustrated and described. Additionally, the methods can omit some acts, and/or employ additional acts. One or more of these acts may be performed by or via one or more circuits, for instance one or more processors (e.g., digital processors such as microprocessors, analog processor such as quantum processors, a hybrid computer including a digital process and an analog processor). A sampling device, such as, a device including a hybrid computer, may perform one or more of the acts in method 300.

Method 300 starts at 302, for example is response to a call from another routine.

At 304, the at least one digital processor initializes a counter i to a starting value, e.g., i=0.

At 306, the at least one digital processor determines an intermediate value s_i* of the annealing parameter s. The intermediate value s_i* may be supplied to the at least one digital processor as part of a set of input or determined by the at least one digital processor.

At 308, the at least one digital processor determines a time interval t_ramp. The value of t_rampmay be supplied to the at least one digital processor as part of a set of input or determined by the at least one digital processor. t_rampis chosen so as to satisfy equation (17).

At 310, the at least digital processor programs the at least one analog processor so that at least a portion of the qubits in the at least one analog processor is in an initial state at an initial time corresponding to counter i. The initial state is an equilibrium state corresponding to a distribution P_θ(z). In one implementation, distribution P_θ(z) is a Boltzmann distribution.

From 312 to 328, method 300 iterates over linear increments of counter i.

At 312, the at least one digital processor programs the at least one analog processor so that the at least analog processor evolves from a value s=1 to a value s_i* of the annealing parameter s within a time t_ramp, where s_i* and t_rampwhere determined at 306 and 308, respectively.

At 314, the at least one digital processor clamps the state of the at least one analog processor for a time t_hold. As mentioned before, the time t_holdshould be long enough for the system to re-equilibrate and recover from any distortion that may have caused the state of the at least one analog processor to deviate from equilibrium. A distortion may be caused by the fast ramp operation at 312. The time t_holdmay be determined by the at least digital processor at 314 or supplied to the at least one digital processor as part of a set of inputs.

At 316, the at least one digital processor programs the at least one analog processor so that the at least one analog processor evolves from s_i* to s=1 in time t_ramp. The time t_rampmay be the same time interval applied at 312 or different, but still satisfying condition (17).

At 318, while the at least analog processor is in a classical state, the at least one digital processor reads out the state of the at least one analog processor. The state of at least a portion of the qubits in the at least one analog processor may be not corresponding to the initial distribution P_θ; at least a portion of the qubits in the at least one analog processor may be in a state {tilde over (z)} corresponding to a distribution {tilde over (P)}_θ that is a deviation from the initial distribution P_θ.

At 320, the at least one digital processor performs post-processing on the state {tilde over (z)} read out at 318. As mentioned before, in some implementations (e.g., small t_ramp), it is possible to recover the initial distribution P_θ from the distribution {tilde over (P)}_θ measured at 318, therefore the state z can be recovered from {tilde over (z)}. In at least one implementation, the at least one digital processor performs QMC post-processing.

At 322, the at least one digital processor calculates importance weights using equation (25).

At 324, the at least one digital processor determines whether an exit condition has been met, as described with reference to 220 of method 200. If the exit condition is not met, control passes to 326, otherwise to 332.

At 326, the at least one digital processor programs the at least one analog processor so that at least a portion of the qubits in the at least one analog processor is in state z, obtained at after post-processing at 320.

At 328, the at least one digital processor determines a new intermediate value s*_i+1for the annealing parameter s, where s*_i+1>s*_i.

At 330, the at least one digital processor increases the value of the counter i to i=i+1, before passing control to 312.

At 332, method 300 terminates, until it is, for example, invoked again.

In method 300, all the samples are generated from post-processing, for example QMC post-processing. Therefore, every time the quantum annealing processor starts from a state z_ithat was the outcome of the last post-processing operation, its output {tilde over (z)}_i+1will be post-processed to generate z_i+1. Method 300 may not be efficient in generating samples due to the need to reinitialize the parameters (e.g. qubits states) of the quantum annealing processor at each iteration.

Quantum Importance Sampling with Sequential Hardware Samples

In one implementation, the quantum annealer is initialized using its own previous state rather than a state generated by quantum Monte Carlo post-processing. One benefit of initializing the quantum annealer using its own previous state is that it can save hardware re-initialization time, and can increase the speed of quantum importance sampling.

Sequential operation can be described by the following:

P
˜

θ



(

z
˜

)

≈

∑

z
˜

′



w


(

z
˜


′

)



P
θ

Q

A



(

z
˜

|

z
˜

′

)



P
˜

θ
′



(

z
˜

′

)

(
26
)

where the weights are given by:

w


(

z
˜

′

)

=

∑

z
˜

′



P
θ



(

z
′

)

P

θ
′



(

z
′

)



P

θ
′

Q

M

C



(

z
′

|

z
˜

′

)

(
27
)

In implementations where θ is not very different from θ′ and t_holdis sufficiently long, the difference in output of P_θ^QAwhen starting from z′ or {tilde over (z)}′ is negligible, i.e. P_σ^QA({tilde over (z)}|z′)≠P_θ^QA({tilde over (z)}|{tilde over (z)}′). The two distributions are not equal (z′ is post-processed at θ′ and P_θ^QAoperates at θ), however, increasing t_holdwill reduce the error in the approximation.

In one implementation, statistical error in estimating the weights w({tilde over (z)}′) can be reduced by generating many quantum Monte Carlo samples of z′ for each {tilde over (z)}′ generated by the quantum annealer at parameter θ′.

Using (26) eliminates the use of post-processing between successive runs of the quantum annealing operations, so that post-processing is only used for calculating the importance weights, which can be done at the end, after generating all the samples. During the quantum annealing operation, each output of the quantum annealing processor determines the initial state of the next P_θ^QAoperations. Thus the quantum annealing processor is run in a sequential way each time starting with its own previous results without re-initialization.

FIG. 4 is a flow diagram of an exemplary method 400 of operation of a hybrid computing system comprising at least one digital processor and at least one analog processor to generate equilibrium statistics and calculating importance weights by sequential iteration of the analog processor. The analog processor may be a quantum annealing processor.

Method 400 comprises acts 402 through 430. For method 400, as with others methods taught herein, the various acts may be performed in a different order than that illustrated and described. Additionally, the methods can omit some acts, and/or employ additional acts. One or more of these acts may be performed by or via one or more circuits, for instance one or more processors (e.g., digital processors such as microprocessors, analog processor such as quantum processors, a hybrid computer including a digital process and an analog processor). A sampling device, such as, a device including a hybrid computer, may perform one or more of the acts in method 400.

Method 400 start at 402, for example is response to a call from another routine.

At 404, the at least one digital processor initializes a counter i to a starting value, e.g., i=0.

At 406, the at least one digital processor determines an intermediate value s_i* of the annealing parameter s, as described with reference to 306 of method 300.

At 408, the at least one digital processor determines a time interval t_ramp, as described with reference to 308 of method 300.

At 410, the at least digital processor programs the at least one analog processor so that at least a portion of the qubits in the at least one analog processor is in an initial state at an initial time corresponding to counter I, as described with reference to 310 of method 300.

From 412 to 424, method 400 iterates over linear increments of counter i.

At 412, the at least one digital processor programs the at least one analog processor so that the at least analog processor evolves from a value s=1 to a value s_i* of the annealing parameter s within a time t_ramp, where s_i* and t_rampwhere determined at 406 and 408, respectively.

At 414, the at least one digital processor clamps the state of the at least one analog processor for a time t_hold, as described with reference to 314 of method 300.

At 416, the at least one digital processor programs the at least one analog processor so that the at least one analog processor evolves from s_i* to s=1 in time t_ramp, as described with reference to 316 of method 300.

At 418, while the at least analog processor is in a classical state, the at least one digital processor reads out the state of the at least one analog processor, as described with reference to 318 of method 300.

At 420, the at least one digital processor stores the state of the at least one analog processor. The at least one digital processor may employ a memory storage system as described in FIG. 1.

At 422, the at least one digital processor determines whether an exit condition has been met, as described with reference to 220 of method 200. If the exit condition is not met, control passes to 424, otherwise to 428.

At 424, the at least one digital processor determines a new intermediate value s*_i+1for the annealing parameter s, where s*_i+1>s*_i.

At 426, the at least one digital processor increases the value of the counter i to i=i+1, before passing control to 412.

At 428, the at least one digital processor calculates importance weights based on equation (27) using the values stored at 420.

At 430, method 400 terminates, until it is, for example, invoked again.

Quantum Update Operator

The efficiency of Markov Chain Monte Carlo (MCMC) methods tends to be limited by the update operation used. For instance, update operations such as Metropolis-Hastings or Gibbs sampling tend to be highly local, which may impede exploration of the state space and thus reduce the quality of samples. Some MCMC methods attempt to mitigate this locality effect (e.g. as in parallel tempering), but the efficiency of such methods still tends to be strongly dependent on the power of the update operator.

In some implementations, a re-equilibration by a physical quantum annealer is used to provide a broad-scale update operator in an MCMC method sampling from classical and/or quantum Boltzmann distributions. The re-equilibration comprises an operation of the physical quantum annealer which uses thermal and/or quantum effects to evolve (or at least tend to evolve) the state of the physical quantum annealer towards equilibrium. This operation may be considered to be an update operation and may be denoted U_RE: Ψ→Ψ, where Ψ is the space of all quantum states representable by the physical quantum annealer. Thus, ψ′=U_REψ denotes the quantum state ψ being transformed to updated quantum state ψ′ according to operator U_RE. Further operators introduced herein use similar notation.

In some implementations the re-equilibration operation U_REcomprises holding the parameters of a physical quantum annealer at a particular point in the anneal s* for time t_hold, during which time the quantum and/or thermal dynamics of the system will tend to evolve the quantum state toward an equilibrium quantum Boltzmann distribution with parameters A(s), B(s), and temperature T (corresponding to the physical temperature of the physical quantum annealer). This operation may be denoted by the operator U_s*→s*^PQA(t_hold): Ψ→Ψ.

The re-equilibrating behavior of U_s→s^PQA(t_hold) may be at least partially explained by the observation that open quantum systems have a state described by a density matrix p(t) and, in a Markovian environment in particular, this state evolves with time according to:

{dot over (ρ)}(t)= [Image Omitted] (t)*ρ(t) (28)

where [Image Omitted] is a superoperator in the Hilbert space of the density matrices. satisfies the detailed balance condition, which results in the behavior of U_s*→s*^PQA(t_hold) described above.

One may extend U_RE(e.g. U_s*→s*^PQA(t_hold)) by composing it with other operators without sacrificing the equilibration dynamics described above, so long as those compositions satisfy detailed balance conditions. For instance, in some implementations, U_REis composed with an invertible operation U_Ito generate a composed update operation U_c=U_I⁻¹U_REU_I. This pattern of composition may be repeated with further invertible operations. Note that, since MCMC methods and quantum annealing are probabilistic methods, invertible operation U_Imay be invertible in a probabilistic sense but not necessarily invertible in an exact sense (i.e. U_I⁻¹U_Iwill not necessarily yield an identity output for every input).

In some implementations, the physical quantum annealer is initialized in a classical state z corresponding to s=1, evolved to a quantum state ψ, clamped for a period of time t_holdyielding a re-equilibrated quantum state ψ′, and then annealed to a classical state z′ corresponding to s=1. The resulting classical state z′ may be considered to be a re-equilibrated update to classical state z. An example method implementing such a process is shown in FIG. 2 and illustrated by FIG. 5A, both of which are described above.

This may be more formally expressed by denoting the evolution from s=1 to s=s* over time t_rampas U_1→s*^PQA(t_ramp):Z→Ψ and denoting the evolution from s=s* to s=1 over time t_rampas U_s*→1^PQA(t_hold):Ψ→Z, where Z is the space of classical states. These operations are inverses of each other and, in essence, allow for the conversion of classical states to quantum states and vice-versa. Accordingly, the re-equilibration operator U_s*→s*^PQA(t_hold) may be extended from the quantum space to the classical space by performing the following operation:

z′=U_1→s*^PQA(t_ramp)U_s*→s*^PQA(t_hold)U_s*→1^PQA(t_ramp)z (29)

The above operation may be further extended to accommodate classically-implemented MCMC annealing methods where the state that must be updated is a classically-represented intermediate state x corresponding to a point s=s* in the evolution of the MCMC method. Although such a state x may notionally correspond to a quantum state, it may not necessarily be directly representable on the physical quantum annealer (e.g. because the physical quantum annealer cannot be initialized in a quantum state, because the state x comprises some internal representation specific to the MCMC method, and/or for other reasons).

In some implementations, the intermediate state x is converted to a classical state z via a classical update operation U_s*→1^MCMC:X→Z, updated to an updated classical state z′ by applying a suitable update operation (e.g. the operation described in equation (29)), and then converted to an updated intermediate state x′ via an inverse classical update operation U_1→s*^MCMC:Z→X, where X is the space of intermediate states representable by the MCMC method. This yields the following update operation:

x′=U_1→s*^MCMCU_s*→1^PQA(t_ramp)U_s*→s*^PQA(t_hold)U_1→s*^PQA(t_ramp)U_s*→1^MCMCx (30)

As described above, this operation preserves detailed balance conditions and so may be used as an update operation for at least some classical MCMC methods, including those directed to sampling from a Boltzmann distribution.

In some implementations, U_s*→1^MCMC:X→Z comprises annealing quickly to the end of the evolution using classical updates (e.g. Metropolis-Hastings, Gibbs sampling, or other classical updates employed by the relevant MCMC method). The specific speed meant by “quickly” will tend to vary between MCMC methods, update operations, and intermediate states s*, but will generally be faster than the ordinary annealing rate of the method. Preferably, the annealing rate is selected to be sufficiently fast to limit changes introduced over the course of an operation z=U_s*→1^MCMCx to local changes. Similarly, the operator U_1→s*^MCMCmay comprise reverse annealing quickly back to point s*, although in some embodiments the annealing rate of U_1→s*^MCMCis slower than the annealing rate of U_s*→1^MCMCto provide a greater opportunity for UM to invert changes introduced by U_sn^mj. The annealing rate and other characteristics of operation U_s*→1^MCMCand U_1→s*^MCMCmay be pre-selected, provided by a user, and/or otherwise obtained. Accelerating the annealing rate may involve, for example, relaxing one or more proposal acceptance criteria.

In some implementations, update operation U_s*→1^MCMC:X→Z comprises finding a projection from an intermediate state x to a classical state z without necessarily annealing. For example, in an example MCMC method simulating quantum annealing an intermediate state x may comprise a set of i spins {s_qⁱ} for each simulated qubit q. A projection may be determined by randomly selecting, for each simulated qubit, one of its associated spins in state x. That is, U_s*→1^MCMCX→Z may be implemented according to:

z=U_1→s*^MCMCx={s_q^j|j=random{1, . . . ,i}∀q} (31)

In some implementations, X is a higher-dimensional space than Z, so the inverse operation U_1→s*^MCMC:Z→X may comprise a non-projection operation which maps from Z to X and inverts the operation of U_s*→1^MCMC:X→Z. For example, in the example MCMC method described above, U_1→s*^MCMC:Z→X may comprise generating an initial intermediate state x_i′ where each spin s_qⁱfor a simulated qubit q is set based on the corresponding classical value s_qfor that qubit in classical state z′. U_1→s*^MCMC:Z→X may further comprise generating an output intermediate state x_iby reverse-annealing from x_i′ as described above. This further reverse annealing step, although not necessarily required, is likely to improve the performance of the MCMC method, since the example initial intermediate state x_i′ is relatively homogenous.

In addition, or alternatively, operations U_s*→1^MCMC:X→Z and/or U_1→s*^MCMC:Z→X may comprise any other suitable operation for converting an intermediate state to a classical state.

A chart 500B illustrating an example execution of the update operation of equation (30) is illustrated in FIG. 5B. Like-numbered elements of FIGS. 5A and 5B are similar and are described above with reference to FIG. 5A. Lines 502 and 510 correspond to evolutions according to operations U_s*→1^MCMCand U_1→s*^MCMCover time periods 512 and 520, respectively. (For consistency with FIG. 5A, FIG. 5B refers to point s*_iin the anneal, although this is simplified to s* in the following discussion for convenience). In some implementations, the evolutions of lines 502 and 510 are performed by a classical computer using a classical update operator associated with an MCMC method, as described above. The intervening evolutions of lines 504, 506, and 508 may be performed by a physical quantum annealer.

FIG. 6 shows a flowchart illustrating an example method 600 for implementing an update operation for a classical MCMC method using re-equilibration, for example via a physical quantum annealer. At 602, an intermediate state x is obtained. This may comprise, for example, receiving a communication containing intermediate state x from another computing system, determining intermediate state x in the course of executing an MCMC method, and/or otherwise receiving intermediate state x.

At 604, intermediate state x is converted to a classical state z. In some implementations, this comprises forward-annealing x from s=s* to s=1 to yield classical state z. As described above, this forward annealing may be done quickly so as to limit changes in global state introduced during annealing. The forward-annealing may be performed, for example, by a classical computer implementing the same or a similar method to the MCMC method which generated intermediate state x (subject to any changes to the method involved in accelerating the annealing rate). In some implementations, this conversion act comprises finding a projection of x to yield classical state z, as described above.

At 606, classical state z is reverse annealed (e.g. via a ramp operation over time t_ramp) from s=1 to s=s* in the evolution to yield quantum state ψ. This reverse annealing may be executed by a physical quantum annealer. For instance, classical state z may be read out from a classical processor, transmitted to a physical quantum annealer (e.g. a physical quantum annealer in a hybrid computer which also comprises the classical processor), and reverse-annealed by the physical quantum annealer.

At 608, the evolution of quantum state ψ is paused for a period of time (e.g. for a period t_hold), during which thermal and/or quantum effects may induce re-equilibration. This operation yields re-equilibrated quantum state ψ′. This re-equilibration may be executed by a physical quantum annealer. Note that quantum states ψ and ψ′ are not necessarily observed or read out during the course of method 600, and are named merely for convenience of explanation.

At 610, re-equilibrated quantum state ψ′ is forward-annealed (e.g. via a ramp operation over time t_ramp, which may be the same as or different than the time period of reverse-annealing act 606) s=s* to s=1. This operation yields re-equilibrated classical state z′. This forward annealing may be executed by a physical quantum annealer.

At 612, re-equilibrated classical state z′ is reverse-annealed from s=1 to s=s* to yield re-equilibrated intermediate state x′. This reverse-annealing may be performed, for example, by a classical computer implementing the same or a similar method to the MCMC method of act 604 (subject to any changes involved in reversing the annealing and, optionally, using a different annealing rate). For instance, re-equilibrated classical state z′ may be read out from the physical quantum annealer, transmitted to a classical processor (e.g. a classical processor in a hybrid computer which also comprises the physical quantum annealer), and reverse-annealed by the classical processor.

At 614, re-equilibrated classical state x′ is used as a basis for an update proposal for the MCMC method which generated intermediate state x. For instance, the MCMC method may use re-equilibrated intermediate state x′ as an update proposal with or without further post-processing and/or other modification. The resulting update proposal may be accepted or rejected according to the criteria of MCMC method.

The above method provides an exemplary implementation of a quantum update operator for a classical MCMC method and illustrates a use of physical quantum annealers to refine classical methods.

The above description of illustrated embodiments, including what is described in the Abstract, is not intended to be exhaustive or to limit the embodiments to the precise forms disclosed. Although specific embodiments of and examples are described herein for illustrative purposes, various equivalent modifications can be made without departing from the spirit and scope of the disclosure, as will be recognized by those skilled in the relevant art. The teachings provided herein of the various embodiments can be applied to other methods of quantum computation, not necessarily the exemplary methods for quantum computation generally described above.

The various embodiments described above can be combined to provide further embodiments. All of the U.S. patents, U.S. patent application publications, U.S. patent applications, foreign patents, foreign patent applications and non-patent publications referred to in this specification and/or listed in the Application Data Sheet including U.S. patent application Ser. No. 15/870,411, U.S. Patent Publication No. US20150161524, U.S. Patent Publication No. US20150269124, International Patent Application No. US2016/059,169, U.S. Provisional Application No. 62/445,496, U.S. Provisional Application No. 62/447,904, and U.S. Provisional Application No. 62/417,940 are incorporated herein by reference, in their entirety. Aspects of the embodiments can be modified, if necessary, to employ systems, circuits, and concepts of the various patents, applications, and publications to provide yet further embodiments.

These and other changes can be made to the embodiments in light of the above-detailed description. In general, in the following claims, the terms used should not be construed to limit the claims to the specific embodiments disclosed in the specification and the claims, but should be construed to include all possible embodiments along with the full scope of equivalents to which such claims are entitled. Accordingly, the claims are not limited by the disclosure.