Systems and methods for improving the computational efficiency of processor-based devices when solving constrained quadratic models

Abstract

A system and method for optimization algorithms, updating samples, and penalizing constraint violations are discussed. The method for updating samples includes receiving a problem definition having an objective function and a constraint function, an initial sample, and values ​​of a progress parameter. For each variable, a total energy change is determined based on the sample value of the variable and an objective energy change based on one or more terms of the objective function that includes the variable, and a constraint energy change based on the sample value of the variable and each of the constraint functions defined by the variable. A sampling distribution is selected based on the variable type, and update values ​​are sampled based on the total energy change and the progress parameter. An updated sample is returned having an updated value for each variable in the set of variables. This can improve the operation of a processor-based system.

Description

[Technical Field]

[0001] Field The present disclosure relates generally to systems and methods for improving the efficiency of processor-based devices in solving constrained quadratic models, for example, employing a hybrid approach that utilizes digital and analog processors. [Background technology]

[0002] background quantum devices Quantum devices are structures in which quantum mechanical effects can be observed. Quantum devices include circuits in which current transport is governed by quantum mechanical effects. Such devices include spintronics and superconducting circuits. Both spin and superconductivity are quantum mechanical phenomena. Quantum devices can be used in measuring instruments, computers, etc.

[0003] quantum computing A quantum computer is a system that directly uses at least one quantum mechanical phenomenon, such as superposition, tunneling, and entanglement, to perform operations on data. The elements of a quantum computer are qubits. Quantum computers can speed up certain classes of computational problems, such as those that simulate quantum physics.

[0004] quantum processor The quantum processor may take the form of a superconducting quantum processor, which may include a number of superconducting qubits and associated local bias devices, and may also include coupling devices (also known as couplers or qubit couplers) that selectively communicatively couple the qubits.

[0005] A quantum processor is any computer processor designed to exploit at least one quantum mechanical phenomenon (e.g., superposition, entanglement, tunneling, etc.) in processing quantum information. Regardless of the specific hardware implementation, all quantum processors encode and manipulate quantum information in quantum mechanical objects or devices called quantum bits or "qubits," all quantum processors employ structures or devices for communicating information between quantum bits, and all quantum processors employ structures or devices for reading out the state of at least one quantum bit. A quantum processor may include a large number (e.g., hundreds, thousands, millions, etc.) of programmable elements, including, but not limited to, qubits, couplers, readout devices, latching devices (e.g., quantum flux parametron latching circuits), shift registers, digital-to-analog converters, and / or demultiplexer trees, as well as programmable subcomponents of these elements, such as programmable subcomponents for correcting device variations (e.g., inductance tuners, capacitance tuners, etc.), compensating for undesired signal drift, etc.

[0006] Further details and embodiments of exemplary quantum processors that can be used with the systems and devices of the present invention are described, for example, in U.S. Patent Nos. 7,533,068, 8,008,942, 8,195,596, 8,190,548, and 8,421,053.

[0007] Hybrid computing systems including quantum processors A hybrid computing system can include a digital computer communicatively coupled to an analog computer, hi some implementations, the analog computer is a quantum computer and the digital computer is a classical computer.

[0008] A digital computer may include a digital processor that can be used to perform the classical digital processing tasks described in the present systems and methods. A digital computer may include at least one system memory that can be used to store various collections of computer or processor readable instructions, application programs, and / or data.

[0009] A quantum computer may include a quantum processor that includes programmable elements such as qubits, couplers, and other devices. The qubits can be read out by a readout system, and the results are communicated to a digital computer. The qubits and couplers can be controlled by a qubit control system and a coupler control system, respectively. In some implementations, the qubit control system and coupler control system can be used to implement quantum annealing on an analog computer.

[0010] quantum annealing Quantum annealing is a computational technique that can be used to find a low-energy state of a system, typically and preferably the system's ground state. This method relies on the fundamental principle that natural systems tend toward low-energy states because they are more stable. Quantum annealing can use quantum effects, such as quantum tunneling, as a source of delocalization to reach an energy minimum.

[0011] Quantum processors can be designed to perform quantum annealing and / or adiabatic quantum computation. An evolution Hamiltonian proportional to the sum of a first term proportional to the problem Hamiltonian and a second term proportional to the delocalized Hamiltonian can be constructed as follows: H E ∝A(t)H P +B(t)H D However, H E is the evolutionary Hamiltonian, and H P is the problem Hamiltonian, and H Dis the nonlocalized Hamiltonian, and A(t), B(t) are coefficients that can control the rate of evolution, typically in the range [0,1].

[0012] In some implementations, a time-varying envelope function may be placed on the problem Hamiltonian. The appropriate delocalized Hamiltonian is

number

number

number

[0013] A general problem Hamiltonian includes a first component proportional to the diagonal single-qubit term and a second component proportional to the diagonal multi-qubit term, and may be of the form:

number

number

[0014] where

number

[0015] Throughout this specification, the terms "problem Hamiltonian" and "final Hamiltonian" are used interchangeably unless the context dictates otherwise. Particular states of a quantum processor are energetically favored, or simply preferred, by the problem Hamiltonian. These include ground states and may also include excited states.

[0016] H in the above two equations D and H P Hamiltonians such as θ can be physically realized in a variety of different ways. A concrete example is realized by the implementation of superconducting qubits.

[0017] sampling Throughout this specification and the appended claims, the terms "sample," "sampling," "sampling device," and "sample generator" are used.

[0018] In statistics, a sample is a subset of a population, i.e., a selection of data taken from a statistical population. In electrical engineering and related fields, sampling involves taking a set of measurements of an analog signal or some other physical system.

[0019] In many fields involving simulation and computation of physical systems, particularly analog computing, the above meanings can blend. For example, a hybrid computer can derive samples from an analog computer. An analog computer, as a source of samples, is an example of a sample generator. The analog computer can be operated to provide samples from a selected probability distribution, which assigns each data point in the population a respective probability of being sampled. The population can correspond to all possible states of the processor, and each sample can correspond to a respective state of the processor.

[0020] Markov Chain Monte Carlo Markov Chain Monte Carlo (MCMC) is a class of computational techniques that includes, for example, simulated annealing, parallel tempering, population annealing, and other techniques. Markov chains can be used, for example, when probability distributions are not available. Markov chains can be described as sequences of discrete random variables and / or as random processes whose states at each time increment depend only on the previous state. If the chain is long enough, the aggregate properties of the chain, such as the mean, can match the aggregate properties of the target distribution.

[0021] A Markov chain can be obtained by proposing new points according to a Markov proposal process (commonly called an "update operation"). The new point is either accepted or rejected. If the new point is rejected, a new proposal is made, and so on. An accepted new point is one that produces probabilistic convergence to the target distribution. Convergence is guaranteed if the proposal and acceptance criteria satisfy detailed equilibrium conditions and the proposal satisfies ergodicity requirements. Furthermore, proposal acceptance can be done so 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 said to have detailed equilibrium. Typically, new points are often local to the previous point.

[0022] The above examples of the prior art and its associated limitations are intended to be illustrative, not exhaustive. Other limitations of the prior art will become apparent to those skilled in the art upon reading this specification and studying the drawings. Summary of the Invention [Means for solving the problem]

[0023] Quick Overview It is generally desirable to improve the computational efficiency and / or accuracy of the operation of processor-based devices, and particularly desirable when using processor-based devices as solvers for solving constrained quadratic models.

[0024] One way to address constraints when solving an optimization problem is to transform the constraints into part of the objective function using slack variables and minimize them over the new objective function. However, in some implementations, this may incur high computational costs, such as long processing time, a large number of slack variables required, and / or a significantly increased solution complexity. This disclosure describes systems and methods useful for improving computational efficiency when solving constrained quadratic problems. In particular, the systems and methods described herein may allow constraints to be specified without the need to transform the constraints and include them in the objective function.

[0025] To provide some flexibility in solving while also ensuring that constraints are satisfied, constraints can be assigned penalty values, and the weights given to constraints can be varied so that constraints can be violated in some situations (such as at the beginning of simulated annealing). When adding constraints to the objective function, it may be necessary to select penalty values ​​for the constraints without any guidance regarding appropriate penalty values ​​for a given problem. This may result in the need to solve the problem multiple times with different penalty values ​​to find a better solution. The systems and methods described herein may also enable the selection of penalty values ​​without solving the problem multiple times with different penalty values, and the dynamic adjustment of penalty values ​​during the solving process.

[0026] The methods and systems described herein advantageously enable computationally efficient solving of optimization problems with constraints by directly processing the constraints, as opposed to converting the constraints into part of the objective function. This may advantageously enable a processor to return a feasible solution to an optimization problem in a time-efficient and scalable manner. This may also enable the processor to generate multiple solutions in parallel. The contribution of the constraints may be implicitly handled by an update process that considers each constraint. Processing the constraint functions directly and deriving their energy values ​​computationally may increase the efficiency of the solution. The methods and systems described herein may also advantageously enable automatic and dynamic adjustment of penalty weights for constraint functions, allowing the search to efficiently guide towards feasibility while returning better solutions. Many real-world problems, such as those solved in industry, are problems with constraints on feasible solutions. The methods and systems described herein may also enable the solution of problems with different types of variable sets, such as a mixture of binary and discrete variables. Sampling may be performed based on a determination of the variable type and taking into account the current energy values ​​and penalties. In some implementations, the methods and systems described herein can be used in combination with hybrid quantum computing techniques to yield improved solutions.

[0027] According to one aspect, there is provided a method of operating a computing system for updating samples in an optimization algorithm to improve convergence to feasibility, the method being executed by a processor and including receiving a problem definition including a set of variables, an objective function defined over the set of variables, and one or more constraint functions, each of the constraint functions being defined by at least one variable of the variable set; receiving sample values ​​of the set of variables and values ​​of a progress parameter; for each variable of the variable set, determining a variable type for the variable; selecting a sampling distribution based on the variable type; determining an objective energy bias based on the sample values ​​of the variable and one or more terms of the objective function including the variable; determining one or more constraint energy biases based on the sample values ​​of the variable and each of the constraint functions defined by the variable; and sampling updated values ​​for the variables from the sampling distribution based on the objective energy bias, the one or more constraint energy biases, and the progress parameter; and returning the updated sample, wherein the updated sample includes an updated value for each variable of the variable set.

[0028] According to another aspect, the method may further include receiving a value of a penalty parameter; sampling an update value of the variable from the sampling distribution may further include sampling an update value of the variable from the sampling distribution based on the value of the penalty parameter; receiving the value of the penalty parameter may include receiving a value of a Lagrangian parameter that depends on the value of the progress parameter; determining a variable type of the variable may include determining that the variable type is one of binary, discrete, integer, or continuous; determining that the variable type is one of binary, discrete, integer, or continuous may include determining that the variable type is binary; and selecting the sampling distribution based on the variable type may include selecting a Bernoulli distribution. selecting a distribution, determining that the variable type is one of binary, discrete, integer, or continuous may include determining that the variable type is discrete, selecting a sampling distribution based on the variable type may include selecting a softmax distribution, determining that the variable type is one of binary, discrete, integer, or continuous may include determining that the variable type is one of integer or continuous, selecting a sampling distribution based on the variable type may include selecting a conditional probability distribution, sampling update values ​​from the sampling distribution may include slice sampling from the conditional probability distribution, and receiving a value for the progression parameter may include receiving an inverse temperature.

[0029] According to one aspect, a system for updating samples in an optimization algorithm is provided, the system including at least one non-transitory processor-readable medium that stores at least one of processor-executable instructions and data, and at least one processor communicatively coupled to the at least one non-transitory processor-readable medium, the at least one processor performing a method described herein in response to executing the at least one of the processor-executable instructions and data.

[0030] According to one aspect, a method of operating a computing system is provided, the computing system including one or more processors, the method being executed by at least one of the one or more processors and including receiving a problem definition including a variable set, an objective function defined over the variable set, and one or more constraint functions, each of the constraint functions being defined by at least one variable in the variable set; initializing a sample solution and progress parameters for the objective function; iteratively incrementing a stage of an optimization algorithm until a termination criterion is met; for each variable in the variable set, selecting an i-th variable from the variable set, the i-th variable having a current value; calculating an objective energy bias for the objective function based on the current value of the i-th variable; calculating a constraint energy bias for each constraint function defined by the i-th variable based on the current value of the i-th variable; and sampling an updated value for the i-th variable based on the objective energy bias and the constraint energy bias, the updated value replacing the current value; sampling; incrementing the progress parameters; evaluating the termination criterion; and outputting a solution including the current values ​​of the variable set.

[0031] According to another aspect, selecting an i-th variable from the set of variables may include selecting a binary variable, the binary variable having a current value and an alternative value; calculating an objective energy bias for the objective function based on the current value of the i-th variable may include calculating a difference in energy of the objective function between the current value and the alternative value of the binary variable; calculating a constraint energy bias based on each constraint function defined by the i-th variable based on the current value of the i-th variable may include calculating a difference in energy of each constraint function defined by the binary variable between the current value of the binary variable and the alternative value of the binary variable; and sampling an update value for the i-th variable based on the objective energy bias and the constraint energy bias may include sampling an update value for the i-th variable based on a difference in energy values ​​of the objective function and each constraint function defined by the i-th variable.

[0032] According to another aspect, the method may further include initializing a penalty parameter, adjusting the penalty parameter for each constraint function based on an energy difference and a progress parameter for the constraint function defined by the i-th variable, wherein calculating the energy difference for each constraint function defined by the i-th variable includes penalizing each constraint function by the penalty parameter, and incrementing a stage of the optimization algorithm with respect to the objective function may include incrementing one of a simulated annealing algorithm or a parallel tempering algorithm, and wherein the problem definition includes the objective function and one or more constraint functions. receiving the problem definition may include receiving a problem definition including a quadratic objective function and one or more quadratic equality or inequality constraint functions; receiving the problem definition including the set of variables may include receiving a problem definition including one or more sets of binary variables, integer variables, or discrete variables; receiving the problem definition including the set of variables may include receiving a problem definition including one or more integer variables; sampling the updated value of the ith variable may include sampling from a conditional probability distribution, the ith variable includes an integer variable from the one or more integer variables; and the termination criteria may include one of a number of iterations, an amount of time, an average change in value limit, or a value of a progress parameter.

[0033] According to one aspect, a system for use in optimization is provided, the system including at least one non-transitory processor-readable medium that stores at least one of processor-executable instructions and data, and at least one processor communicatively coupled to the at least one non-transitory processor-readable medium, the at least one processor performing a method described herein in response to executing the at least one of the processor-executable instructions and data.

[0034] According to another aspect, the system may further include a quantum processor, and after performing the method, the at least one processor may instruct the quantum processor to perform quantum annealing based on the output solution.

[0035] According to one aspect, there is provided a method of operating a hybrid computing system, the hybrid computing system including a quantum processor and a classical processor, the method being executed by the classical processor, comprising receiving a constrained quadratic optimization problem, the constrained quadratic optimization problem including a set of variables, an objective function defined over the set of variables, one or more constraint functions, each of the constraint functions defined by at least one variable of the set of variables, and progress parameters for the optimization, the progress parameters including a set of values ​​that increment between an initial value and a final value; sampling a sample set of values ​​of the set of variables from an optimization algorithm iteratively until a final value of the progress parameters is reached; updating the sample set of values ​​using an update algorithm, wherein for each variable of the variable set, determining a variable type of the variable; updating the sample set of values, including selecting a sampling distribution based on the variable type, determining an objective energy bias based on sample values ​​of the variables from the sample set of values ​​and one or more terms of the objective function including the variables, determining one or more constraint energy biases based on each of the sample values ​​of the variables and the constraint functions defined by the variables, and sampling updated values ​​for the variables from the sampling distribution based on the objective energy bias, the one or more constraint energy biases, and a progress parameter; and returning the updated samples, the updated samples including updated values ​​for each variable in the variable set; incrementing the progress parameter; transmitting the one or more final samples to a quantum processor; instructing the quantum processor to refine the samples; and outputting the solution.

[0036] According to other aspects, transmitting one or more final samples to a quantum processor may include transmitting pairs of samples to the quantum processor, and instructing the quantum processor to refine the samples may include instructing the quantum processor to perform quantum annealing to select between the samples, and the method may further include returning the output solution as a sample set of values ​​for the set of variables as input to an optimization algorithm.

[0037] According to one aspect, a hybrid computing system is provided that includes a quantum processor and a classical processor, at least one non-transitory processor-readable medium that stores at least one of processor-executable instructions and data, and at least one processor communicatively coupled to the at least one non-transitory processor-readable medium, the at least one processor performing a method described herein in response to executing the at least one of the processor-executable instructions and data.

[0038] According to one aspect, there is provided a method of operating a computing system for guiding a search space toward feasibility to improve performance of the computing system, the computing system including one or more processors, the method being executed by at least one of the one or more processors and including receiving a sample from an optimization algorithm, determining an energy value of one or more constraint functions, evaluating the feasibility of the sample, increasing a penalty value if the sample is not feasible, decreasing the penalty value if the sample is feasible, and returning the penalty value to the optimization algorithm.

[0039] According to another aspect, the method may further include determining whether violations have decreased compared to a previous sample and increasing the initial adjuster value if the violations have not decreased, and determining whether a current best solution has improved compared to a previous sample and increasing the initial adjuster value if the current best solution has not improved.

[0040] According to one aspect, a method of operating a computing system is provided, the computing system including one or more processors, the method being executed by at least one of the one or more processors, the method comprising: receiving a problem definition including a set of variables, an objective function defined over the set of variables, and one or more constraint functions, each of the constraint functions being defined by at least one variable of the variable set, the variable set including a first subset of variables having at least one variable that is one of binary, integer, discrete, or continuous, and a second subset of variables having at least one variable that is continuous; initializing a sample solution and progress parameters for the objective function; initializing a continuous problem defined over the second subset of variables, the continuous problem comprising a linear programming model; iteratively incrementing stages of an optimization algorithm until a termination criterion is met; for each variable in the first subset of variables, determining a variable type for the variable; selecting a sampling distribution based on the variable type; determining an objective energy bias based on one or more terms of an objective function including the values ​​and variables; determining one or more constraint energy biases based on the sample values ​​of the variables and each of the constraint functions defined by the variables; sampling updated values ​​of the variables from the sampling distribution based on the objective energy bias, the one or more constraint energy biases, and a progress parameter; returning the updated samples, where the updated samples include updated values ​​for each variable in the set of variables; solving a linear programming model for a second subset of variables while holding the values ​​of each variable in the first subset of variables fixed to the updated values; sampling updated values ​​for each variable in the second subset of variables based on the solved linear programming model, where the updated values ​​replace the current values; updating the objective energy bias and one or more constraint energy biases based on the updated values ​​of the second subset of variables; incrementing the progress parameter; evaluating a termination criterion; and outputting a solution including the updated values ​​of the set of variables.

[0041] According to another aspect, the method may further include initializing one or more penalty parameters; sampling update values ​​of the variables from the sampling distribution further includes sampling update values ​​of the variables from the sampling distribution based on at least one value of the one or more penalty parameters; sampling update values ​​of each variable in the second subset of variables includes sampling update values ​​of each variable in the second subset of variables based on the solved linear programming model and the at least one value of the one or more penalty parameters; initializing the one or more penalty parameters may include initializing at least one Lagrangian parameter that depends on the value of the progress parameter; determining a variable type of the variables may include determining that the variable type is one of binary, discrete, integer, or continuous; determining that the variable type is one of binary, discrete, integer, or continuous may include determining that the variable type is binary; and selecting the sampling distribution based on the variable type. may include selecting a Bernoulli distribution, determining that the variable type is one of binary, discrete, integer, or continuous may include determining that the variable type is discrete, selecting a sampling distribution based on the variable type may include selecting a softmax distribution, determining that the variable type is one of binary, discrete, integer, or continuous may include determining that the variable type is one of integer and continuous, selecting a sampling distribution based on the variable type may include selecting a conditional probability distribution, sampling updated values ​​for the variable from the sampling distribution may include slice sampling from the conditional probability distribution, receiving values ​​for the progression parameter may include receiving an inverse temperature, incrementing a stage of the optimization algorithm for the objective function may include incrementing one of a simulated annealing algorithm or a parallel tempering algorithm, receiving a problem definition including the objective function and one or more constraint functionsThe method may include receiving a problem definition including a quadratic objective function and one or more quadratic equality or inequality constraint functions, and the termination criteria may include one of a number of iterations, an amount of time, an average change in a value limit, or a value of a progress parameter.

[0042] According to one aspect, a system for use in optimization is provided, the system including at least one non-transitory processor-readable medium that stores at least one of processor-executable instructions and data, and at least one processor communicatively coupled to the at least one non-transitory processor-readable medium, the at least one processor performing a method described herein in response to executing the at least one of the processor-executable instructions and data.

[0043] According to another aspect, the system may further include a quantum processor, and after performing the methods described herein, the at least one processor instructs the quantum processor to perform quantum annealing based on the output solution.

[0044] In other aspects, the above-described features can be combined together in any reasonable combination, as will be appreciated by one of skill in the art.

[0045] BRIEF DESCRIPTION OF THE DRAWINGS In the figures, identical reference numbers identify similar elements or acts. The sizes and relative positions of elements in the figures are not necessarily drawn to scale. For example, the shapes and angles of various elements are not necessarily drawn to scale, and some of these elements may be arbitrarily enlarged and positioned for clarity of the drawings. Furthermore, the particular shapes of elements as drawn are not intended to convey any information regarding the actual shape of the particular elements, but may have been selected merely to facilitate ease of understanding the drawings. [Brief explanation of the drawings]

[0046] [Figure 1]1 is a schematic diagram of a hybrid computing system including a digital computer coupled to an analog computer in accordance with the present systems, devices, and methods. [Figure 2] 1 is a schematic diagram of a portion of an exemplary superconducting quantum processor. [Figure 3] 1 is a flow diagram of an example of a method for performing sampling in an optimization algorithm. [Figure 4] 1 is a flow diagram of an example of a method of operation of a computing system for finding a solution to a constrained problem. [Figure 5] 1 is a flow diagram of another example method of operation of a computing system for finding a solution to a constrained problem. [Figure 6] 1 is a flow diagram of an example method of operating a hybrid computing system. [Figure 7] FIG. 1 is an illustration of functions of constraints on equalities and inequalities. [Figure 8] 1 is a flow diagram of an example of a method for adjusting penalty parameters to steer a search space toward feasibility. [Figure 9] 10 is a flow diagram of another example of a method for adjusting penalty parameters to steer a search space toward feasibility. [Figure 10] 1 is a diagram of the domain of the variable xk used in calculating the partition function. [Figure 11] 1 is a diagram of the general case of the partition function of segment s. [Figure 12] FIG. 10 is an example of a probability distribution of a binary variable based on energy bias and progress parameters. [Figure 13] FIG. 1 is a diagram of an example of a probability distribution function for a continuous variable. [Figure 14] 1 is a flow diagram of an example of a method of operation of a computing system for finding a solution to a constrained problem with continuous variables. DETAILED DESCRIPTION OF THE INVENTION

[0047] Detailed Description In the following description, certain specific details are set forth to provide a thorough understanding of the various disclosed implementations. However, those skilled in the art will recognize that implementations can be practiced without one or more of these specific details, or with other methods, components, materials, etc. In other instances, well-known structures related to computer systems, server computers, and / or communication networks are not shown or described to avoid unnecessarily obscuring the description of implementations.

[0048] Unless the context requires otherwise, throughout this specification and the appended claims, the word "comprising" is synonymous with "including" and is inclusive or open-ended (i.e., does not exclude additional, unrecited elements or method acts).

[0049] References throughout this specification to "one implementation" or "an implementation" mean that a particular feature, structure, or characteristic described in connection with that implementation is included in at least one implementation. Thus, the appearances of the phrase "in one implementation" or "in an implementation" in various places throughout this specification do not necessarily all refer to the same implementation. Furthermore, particular features, structures, or characteristics may be combined in any suitable manner in one or more implementations.

[0050] As used in this specification and the appended claims, the singular forms "a," "an," and "the" include plural referents unless the context clearly dictates otherwise. It should also be noted that the term "or" is generally used in its sense to include "and / or" unless the context clearly dictates otherwise.

[0051] The headings provided herein and the Abstract of this disclosure are for convenience only and do not interpret the scope or meaning of the implementations.

[0052] As an illustrative example, a superconducting quantum processor designed to perform adiabatic quantum computing and / or quantum annealing will be used in the following description. However, as previously mentioned, those skilled in the art will understand that the present systems and methods may be applied to any form of quantum processor hardware (e.g., superconducting, photonic, ion trap, quantum dot, topological, etc.) that implements any form of quantum algorithm (e.g., adiabatic quantum computing, quantum annealing, gate / circuit-based quantum computing, etc.). A quantum processor may be used in combination with one or more classical or digital processors. The methods described herein may be performed by a classical or digital processor in communication with a quantum processor that implements a quantum algorithm.

[0053] Exemplary Computing System 1 illustrates a computing system 100 that includes a digital computer 102. This example digital computer 102 includes one or more digital processors 106 that can be used to perform classical digital processing tasks. The digital computer 102 can further include at least one system memory 122 and at least one system bus 120 that couples various system components, including the system memory 122, to the digital processor 106. The system memory 122 can store one or more sets of processor-executable instructions, sometimes referred to as modules 124.

[0054] The digital processor 106 may be any logic processing device or circuit (e.g., integrated circuit), 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"), and / or combinations thereof.

[0055] In some implementations, computing system 100 includes analog computer 104, which may include one or more quantum processors 126. Quantum processor 126 may include at least one superconducting integrated circuit. Digital computer 102 may communicate with analog computer 104, for example, via controller 118. As described in more detail herein, certain calculations may be performed by analog computer 104 at the direction of digital computer 102.

[0056] The digital computer 102 may include a user input / output subsystem 108. In some implementations, the user input / output subsystem includes one or more user input / output components, such as a display 110, a mouse 112, and / or a keyboard 114.

[0057] The system bus 120 may employ any known bus structure or architecture, including a memory bus with a memory controller, a peripheral bus, and a local bus. The system memory 122 may include read-only memory ("ROM"), static random access memory ("SRAM"), non-volatile memory such as flash NAND, and volatile memory (not shown) such as random access memory ("RAM").

[0058] The digital computer 102 may also include other non-transitory computer-readable or processor-readable storage media or non-volatile memory 116. The non-volatile memory 116 may take various forms, including a hard disk drive for reading from and writing to a hard disk (e.g., a magnetic disk), an optical disk drive for reading from and writing to a removable optical disk, and / or a solid-state drive (SSD) for reading from and writing to a solid-state medium (e.g., a NAND-based flash memory). The non-volatile memory 116 may communicate with the digital processor via a system bus 120 and may include an appropriate interface or controller 118 coupled to the system bus 120. The non-volatile memory 116 may serve as long-term storage of processor-readable or computer-readable instructions, data structures, or other data (sometimes referred to as program modules or modules 124) for the digital computer 102.

[0059] While the digital computer 102 has been described as using hard disks, optical disks, and / or solid-state storage media, those skilled in the art will appreciate that other types of non-transitory and non-volatile computer-readable media may be used. Those skilled in the art will appreciate that some computer architectures use non-transitory volatile memory and non-transitory non-volatile memory. For example, data in volatile memory may be cached in non-volatile memory, or solid-state disks that use integrated circuits to provide the non-volatile memory.

[0060] Various processor or computer readable and / or executable instructions, data structures, or other data may be stored in system memory 122. For example, system memory 122 may store instructions for communicating with remote clients and scheduling the use of resources, including resources on digital computer 102 and analog computer 104. Further, for example, system memory 122 may store at least one of processor-executable instructions or data that, when executed by at least one processor, cause the at least one processor to perform various algorithms for executing instructions. In some implementations, system memory 122 may store processor or computer-readable computational instructions and / or data for performing pre-processing, co-processing, and post-processing for analog computer 104. System memory 122 may store a set of analog computer interface instructions for interacting with analog computer 104. For example, system memory 122 may store processor or computer readable instructions, data structures, or other data that, when executed by a processor or computer, cause the processor or computer to perform one, more, or all of the actions of methods 300 (FIG. 3) through 900 (FIG. 9) and 1400 (FIG. 14).

[0061] Analog computer 104 may include at least one analog processor, such as quantum processor 126. Analog computer 104 may be provided in an isolated environment, for example, an isolated environment that shields the internal elements of the quantum computer from heat, magnetic fields, and other external noise. The isolated environment may include a refrigerator, for example, a dilution refrigerator, operable to cryogenically cool the analog processor, for example, to a temperature below about 1 K.

[0062] Analog computer 104 may include programmable elements such as qubits, couplers, and other devices (also referred to herein as controllable devices). Qubits can be read out by readout system 128. The readout results can be transmitted to other computer-readable instructions or processor-readable instructions of digital computer 102. Qubits can be controlled by qubit control system 130. Qubit control system 130 can include on-chip digital-to-analog converters (DACs) and analog lines operable to apply biases to target devices. Couplers coupling qubits can be controlled by coupler control system 132. Coupler control system 132 can include tuning elements such as on-chip DACs and analog lines. Qubit control system 130 and coupler control system 132 can be used to implement the quantum annealing schedules described herein on analog processor 104. Programmable elements can be included in quantum processor 126 in the form of integrated circuits. Qubits and couplers can be disposed within a layer of an integrated circuit comprising a first material. Other devices, such as readout control system 128, may be located in other layers of the integrated circuit that include the second material. In accordance with the present disclosure, quantum processors, such as quantum processor 126, may be designed to perform quantum annealing and / or adiabatic quantum computing. Examples of quantum processors are described in U.S. Patent No. 7,533,068.

[0063] Exemplary Superconducting Quantum Processor Figure 2 is a schematic diagram of a portion of an exemplary superconducting quantum processor 200 according to at least one implementation. Superconducting quantum processor 200 may be implemented within computing system 100 of Figure 1 to form all or part of quantum processor 126. The portion of superconducting quantum processor 200 shown in Figure 2 includes two superconducting qubits 201 and 202. Also shown is tunable coupling (diagonal coupling) between qubits 201 and 202 by coupler 210 (i.e., providing a two-local interaction). Although the portion of quantum processor 200 shown in Figure 2 includes only two qubits 201, 202 and one coupler 210, those skilled in the art will understand that quantum processor 200 may include any number of qubits and any number of couplers coupling information therebetween.

[0064] Quantum processor 200 includes a number of interfaces 221-225 that are used to configure and control the state of quantum processor 200. Each of interfaces 221-225 may be implemented by a respective inductively coupled structure as shown, as part of a programming subsystem and / or evolution subsystem. Alternatively or additionally, interfaces 221-225 may be implemented by a galvanically coupled structure. In some implementations, one or more of interfaces 221-225 may be driven by one or more DACs. Such programming and / or evolution subsystems may be separate from quantum processor 200, or may be included locally (i.e., on-chip with quantum processor 200).

[0065] In operation of quantum processor 200, interfaces 221 and 224 can be used to couple flux signals into compound Josephson junctions 231 and 232 of qubits 201 and 202, respectively, thereby introducing a tunable tunneling term (Δ i This coupling realizes the off-diagonal σ xterms, and these flux signals are examples of “delocalized signals.” Examples of Hamiltonians (and their terms) used in quantum computing are described in more detail, for example, in U.S. Patent Application Publication No. 2014 / 0344322.

[0066] Similarly, interfaces 222 and 223 can be used to apply flux signals into the qubit loops of qubits 201 and 202, respectively, thereby changing the h in the system Hamiltonian i This coupling realizes the diagonal σ z Additionally, interface 225 can be used to couple the flux signal into coupler 210, thereby providing the J ij (a dimensionless local field for the coupler). This coupling is

number

[0067] In Figure 2, the contribution of each of the interfaces 221-225 to the system Hamiltonian is shown within boxes 221a-225a, respectively. As shown, in the example of Figure 2, boxes 221a-225a are elements of a time-varying Hamiltonian for quantum annealing and / or adiabatic quantum computation.

[0068] Although Figure 2 shows only two physical qubits 201, 202, one coupler 210, and two readout devices 251, 252, a quantum processor (e.g., processor 200) can use any number of qubits, couplers, and / or readout devices, including a larger number (e.g., hundreds, thousands, or more) of qubits, couplers, and / or readout devices. Examples of superconducting qubits include superconducting flux qubits, superconducting charge qubits, etc. In superconducting flux qubits, the Josephson energy dominates or is equal to the charging energy. In charge qubits, this is reversed. Examples of flux qubits that can be used include radio frequency superconducting quantum interference devices, which include a superconducting loop interrupted by one Josephson junction, persistent current qubits, which include a superconducting loop interrupted by three Josephson junctions, etc.

[0069] Constrained quadratic problem A quadratic function refers to a polynomial function with one or more variables that interact with each other at most quadratically. Many real-world problems can be expressed as a quadratic function to be optimized combined with some constraints placed on the feasible outcomes of the variables. In other words, the quadratic function (also referred to herein as the objective function) defines the interactions between the variables, which are affected by a constraint or set of constraints. To obtain an optimal or near-optimal solution to a given problem, the objective function can be minimized or maximized subject to the constraints on the feasible outcomes. In the implementation described below, the problem is structured so that minimization or low energy corresponds to an improved solution. However, it will be understood that a minimization problem can alternatively be structured as a maximization problem (e.g., by reversing the sign of the function), and the following principles generally apply to situations where the objective function is extremized. While minimization is discussed below for clarity, it will be understood that similar principles apply to functions that are maximized, the term "extremization" can be substituted for "minimization," and the energy penalties described below apply to penalize deviations from improved solutions.

[0070] One way to deal with constraints when solving optimization problems is to transform the constraints into part of the objective function using penalty terms and, in the case of inequalities, additional "slack" variables and attempt to minimize over the new objective function. However, in some implementations, this may incur high computational costs, such as long processing times, a large number of slack variables required, and / or a significantly increased solution complexity. This disclosure describes systems and methods useful for improving computational efficiency when solving constrained quadratic problems. In particular, the systems and methods described herein may allow constraints to be specified without the need to transform the constraints and include them in the objective function.

[0071] To provide some flexibility in solving, which may beneficially increase the probability of finding a global optimum while also ensuring that the constraints are satisfied by the final solution, constraints may be assigned penalty values ​​and weights given to constraints may be varied so that constraints may be violated in some situations (such as at the beginning of simulated annealing). When adding constraints to the objective function, it may be necessary to select penalty values ​​for the constraints without any guidance regarding appropriate penalty values ​​for a given problem. This may result in the need to solve the problem multiple times with different penalty values ​​to find a better solution. The systems and methods described herein may also enable the selection of penalty values ​​without solving the problem multiple times with different penalty values ​​and the dynamic adjustment of penalty values ​​during the solving process.

[0072] The desired result of solving the constrained optimization problem is g c The goal of a quadratic optimization problem is to optimize some function, referred to herein as f(x), subject to some inequality or equality constraint, such as (x)≦0. In a quadratic optimization problem, we consider f(x) and g c (x) is a function with polynomial interactions of up to second order. The optimization is

number

number

number

[0073] In a binary implementation with linear equality constraints, the linear constraints can be expressed as functions:

number

[0074] To treat this constraint as a quadratic term, we can square the left-hand side (violation from the desired value) (n=2). If the violation is small, it may be beneficial to use a linear term (n=1) or a constant term (n=0) instead, especially if the violation is less than 1. As discussed in more detail below, the variable x i For each proposed update to E is determined. Since the desired result is to minimize the energy of the system, increasing energy values ​​penalize solutions that yield those increasing values.

[0075] In implementations with linear inequality constraints, functions are treated similarly, but in this case violations are penalized only if the function is positive, and not if the function is negative. So, for example, Σ i A c , i x i +b c ≦0 In the inequality expressed as, ,violation is penalized only if the error function is positive, otherwise no penalty is applied.

[0076] In an implementation with quadratic constraints, the function can similarly be expressed as:

number

[0077] As mentioned above, the value of n can be selected based on the value of the violation. Inequality constraints can be penalized only for a given range of values, as discussed in more detail below.

[0078] 3 is a flow diagram of an example method 300 for performing sampling in an optimization algorithm executed on a processor-based system. Method 300 may be performed on a hybrid computing system including at least one digital or classical processor and a quantum processor, such as hybrid computing system 100 of FIG. 1, or may be performed by a classical computing system including at least one digital or classical processor. Example implementations of method 300 are discussed in more detail below with respect to constrained binary problems.

[0079] Although method 300 includes acts 302-322, those skilled in the art will understand that the number of acts shown is exemplary and that in some implementations, certain acts may be omitted, additional acts may be added, and / or the order of the acts may be changed.

[0080] The method 300 begins at 302, for example, in response to a call or invocation from another routine.

[0081] At 302, a processor receives a problem definition having a set of variables, an objective function defined over the set of variables, and one or more constraint functions, each of the constraint functions defined by at least one variable of the set of variables.

[0082] The problem definition can be expressed as follows: Minimize:Σ i a i x i +Σ i≦j b ij x i x j +c conditions:

number

[0083] At 304, the processor receives sample values ​​of the set of variables and values ​​of the progression parameter. In some implementations, the progression parameter can be an inverse temperature. The sample values ​​can be predetermined starting values ​​based on known characteristics of the problem, can be provided by another routine, can be randomly selected, or can be provided using other techniques known in the art.

[0084] At 306, the processor selects the ith variable from the variable set, which may be, for example, the first variable in the variable set in the first iteration, the second variable in the variable set in the second iteration, etc., until all of the variables in the variable set have been selected. In other implementations, the order of the variables may be selected according to randomization or other metrics known in the art.

[0085] At 308, the processor determines the variable type of the variable. The variable type may be, for example, one of binary, discrete, integer, and continuous. It will be appreciated that for a given problem definition, all of the variables in the variable set may be of a single type, or the problem may have a mix of different variable types.

[0086] At 310, the processor selects a sampling distribution based on the variable type. In some implementations where the variable type is binary, the selected sampling distribution may be a Bernoulli distribution. In other implementations where the variable type is discrete, the selected sampling distribution may be a softmax distribution. Alternatively, in some implementations, the discrete variables may be included as binary variables using a one-hot constraint, and the sampling distribution may be a Bernoulli distribution. The one-hot constraint refers to converting a categorical variable to a binary variable by assigning a dummy binary variable of 0 / 1 or true / false to each category. For example, a discrete variable may be Σi∈D x i =1. As discussed in more detail below with respect to Figures 10 and 11, in other implementations where the variable type is integer, sampling may be performed by Gibbs sampling, i.e., sampling based on conditional probability given the current state. Alternatively, Gibbs sampling can be performed for all variable types, with the conditional probability function determined by the variable type. It will be appreciated that in some implementations, integer variables can also be converted to binary variables, and sampling from a Bernoulli distribution can be performed. As discussed in more detail below, in some implementations where the variable type is continuous, sampling based on conditional probability may be performed by slice sampling, or the sampling distribution may be provided by a linear programming model.

[0087] At 312, the processor determines a target energy bias to act on the variable under consideration given the current values ​​of all other variables. As described above, the optimization problem can be structured with an objective function that defines energy, and during optimization, the processor returns a solution with the objective of reducing this energy. In some implementations, such as when the variable under consideration is a binary variable, the target energy change

number

[0088] At 314, the processor determines the constraint energy bias acting on the variable under consideration given the current values ​​of all other variables and each of the constraint functions that include that variable. As discussed above, the constraint energy bias is Δ E =Σ c (│δ c,1 │ n -│δ c,0 │ n ) The penalty applied to each constraint can be determined by the magnitude of the violation.

[0089] In some implementations, the processor determines the total energy change based on the target energy change and the constraint energy change. As discussed above, for binary variables, the function

number

[0090] At 318, the processor samples an update value for the variable from the sampling distribution based on the objective energy bias and constraint energy bias and the progress parameter. For example, in some implementations, the processor may sample from the distribution such that if the sample value is an improvement in energy and does not violate any constraint functions, the sample value is accepted, whereas if the sample value is not an improvement in energy and / or violates a constraint function, the sample value is accepted only with some probability depending on the progress parameter; otherwise, the previous value is accepted. The progress parameter may tolerate more violations in early stages and less violations in later stages. In other implementations, the sampled update value may be sampled from a weighted distribution based on the total energy change and the progress parameter. For example, if the variable is a binary variable, the sampling distribution may be a Bernoulli distribution, and the variable is

number

number

[0091] Another implementation of Gibbs sampling for binary variables uses a conditional partition function. k The conditional probability that =1 is

number

number

number

number

number

[0092] For discrete variables, Gibbs sampling can be performed based on the use of one-hot constraints. D k is the set of cases for variable k

number

number

[0093] For integer variables, the partition function is

number

number

number

number

number

[0094] At 322, the processor returns updated samples with updated values ​​for each variable in the variable set. Method 300 can then end, for example, until called again, or it can repeat with the new sample values ​​from act 304. The sample output at 322 can be passed to another algorithm for further processing or returned as a solution to the problem. In some implementations, method 300 may be applied after a simulated annealing algorithm returns a solution, where method 300 is applied at zero temperature, causing the system to relax to the nearest local minimum.

[0095] In some implementations, the method 300 can include receiving a value for a penalty parameter, where the total energy change can also depend on the value of the penalty parameter. The penalty parameter is discussed in more detail below. In some implementations, the penalty parameter can be a Lagrangian parameter that depends on the value of the progress parameter.

[0096] An example of pseudocode for method 300 in an implementation with only binary variables is as follows: initialization;

number

[0097] As mentioned above, the progression parameter can be the reverse temperature β. In some implementations, the reverse temperature can be incremented over a set range, such as starting from 0.1 and incrementing to 1. In other implementations, the initial temperature (the highest temperature in the simulated annealing algorithm, and therefore the smallest reverse temperature) β min and the final temperature (the lowest temperature in the simulated annealing algorithm) β maxIt may be beneficial to set based on the probability of moving away from a local minimum. That is, the starting and ending temperatures should be chosen high enough so that there is a significant probability of moving toward a less desirable solution. In some implementations, the solution returned by the processor may be sensitive to the initial temperature, and determining the initial temperature based on the particular problem may beneficially return an improved solution.

[0098] In some implementations, the energy bias contributed by the objective function and any constraints can be used to determine an initial inverse temperature parameter that provides a significant probability (e.g., a 50% probability) of sampling a locally non-optimal value. At the high temperature at the start of the algorithm, this probability can be selected to ensure that even variables with the strongest possible energy bias have a relatively large probability of sampling the least favorable value. At the low temperature at the end of the algorithm, this probability can be selected to ensure that even variables with the smallest possible energy bias have at least a relatively small probability of sampling the least favorable value. In other words, throughout the algorithm, the probability of selecting an unfavorable value starts at a given relatively high probability and then decreases, but is always non-zero.

[0099] 4 is a flow diagram of an example method 400 of operation of a computing system for finding a solution to an optimization problem or other problem that can be expressed as a constrained quadratic model, or for generating sample solutions that can be used as input to other algorithms. Method 400 may be performed on a hybrid computing system that includes at least one digital or classical processor and a quantum processor, such as hybrid computing system 100 of FIG. 1, or may be performed by a classical computing system that includes at least one digital or classical processor. Method 400 may be used in implementations with binary variables, or may include one or more additional acts for representing other variable types as binary variables.

[0100] Although method 400 includes acts 402-430, one skilled in the art will understand that the number of acts shown is exemplary and that in some implementations, certain acts may be omitted, additional acts may be added, and / or the order of the acts may be changed.

[0101] The method 400 begins at 402, for example, in response to a call or invocation from another routine or in response to input by a user.

[0102] At 402, the processor receives a problem definition. The problem definition includes a set of variables, an objective function defined over the set of variables, and one or more constraint functions, each of the constraint functions defined by at least one variable in the set of variables. The problem definition may be received from a user (e.g., entered by an input device), transmitted from another processor, retrieved from memory, or provided as the output of another process executed by the processor. The objective function may be a quadratic function, and the problem definition may define a quadratic optimization problem. The set of variables may include one or more of continuous variables, discrete variables, binary variables, and integer variables. The constraint functions may include one or more quadratic equality or inequality constraint functions.

[0103] At 404, the processor initializes a sample solution to the objective function. The sample solution may be a random solution to the objective function. The random solution may be selected randomly from the entire variable space or may be selected within a limited range of the variable space based on known characteristics of the problem definition or other information. The sample solution may be generated by another algorithm or provided as input by a user.

[0104] At 406, the processor initializes a progression parameter. The progression parameter can be a set of incrementally changing values ​​that define the optimization algorithm. For example, the progression parameter can be an inverse temperature, which can increment from an initial high temperature to a final low temperature. In some implementations, the inverse temperature is provided as a progression parameter for a simulated annealing algorithm. The selection of the inverse temperature is described in more detail above.

[0105] At 408, the processor optionally initializes penalty parameters. As discussed in more detail below, in some implementations, the penalty parameters may be Lagrange multipliers. In other implementations, the penalty parameters may be selected as constant values ​​or values ​​that depend on one or more other variables.

[0106] At 410, the processor optionally calculates the current value of each constraint function at the sample value. For example, the constraint function may be calculated as ΣA c,i x i +b c =0 or as an inequality.

[0107] At 412, the processor increments the stage of the optimization algorithm, such as by providing progress parameters to the optimization algorithm. The optimization algorithm may include simulated annealing, parallel tempering, Markov chain Monte Carlo techniques, branch-and-bound algorithms, and greedy algorithms, which may be executed by a classical computer. The optimization algorithm may also include algorithms executed by a quantum computer, such as quantum annealing, quantum-approximate optimization algorithms (QAOA) or other noisy intermediate-scale quantum (NISQ) algorithms, quantum-implemented fault-tolerant optimization methods, or other quantum optimization algorithms. The quantum computer may include a quantum annealing processor or a gate-model-based processor. Some example implementations of the optimization algorithm are described in U.S. Provisional Patent Application No. 62 / 951,749 and U.S. Patent Application Publication No. 2020 / 0234172. Over successive iterations, the incremented optimization algorithm may provide samples.

[0108] At 414, the processor selects an i-th variable from the variable set. This may be, for example, the first variable in the variable set in the first iteration, the second variable in the variable set in the second iteration, and so on, until all of the variables in the variable set have been selected. The i-th variable has a current value from the sample initialized in act 404. In some implementations, the i-th variable may also have an alternative value. The alternative value may be provided in act 412 by an optimization algorithm or may be generated based on known characteristics of the i-th variable. For example, if the i-th variable is binary and the current value of the i-th variable is 0, the alternative value of the i-th variable will be 1.

[0109] At 416, the processor calculates the energy bias provided by the objective function based on the current value of the i-th variable. In implementations where the variables have alternative values ​​that are binary, this may include calculating the difference in the energy of the objective function between the current value of the i-th variable and the alternative value of the i-th variable. As described above, the objective energy change for a binary variable

number

[0110] At 418, the processor calculates the energy bias provided by the constraint including the i-th variable. There may be one or more constraint functions that contribute to the constraint energy bias. In an implementation with alternative values, this may include calculating the difference in energy of each constraint function defined by the i-th variable between the current value of the i-th variable and the alternative value of the i-th variable. Calculating the difference in energy of each constraint function defined by the i-th variable may include penalizing each constraint function by a respective penalty parameter. For example, for a binary variable with linear equality constraints, the difference in energy of the constraint functions is δ c,0 =δ c -x i A c,i and δ c,1 =δ c +(1-x i )A c,i and the energy difference is given by Δ E =Σ c (│δ c,1 │ n -│δ c,0 │ n ) is given by

[0111] At 420, the processor samples a value for the i-th variable based on the objective energy bias and the constraint energy bias. In implementations with alternative values, this may include sampling based on the difference in energy values ​​of the objective function and each constraint function defined by the i-th variable, as calculated in acts 416 and 418. As described above, the sampled value may be sampled from a weighted distribution based on total energy change and a progress parameter, such as a Bernoulli distribution for binary variables, a softmax distribution for discrete variables, a conditional probability distribution for integer variables, or a linear programming model or slice sampling model for continuous variables.

[0112] At 422, the processor evaluates whether all of the variables in the variable set have been considered. If all of the variables have not been considered, control returns to act 414, where the next variable is selected. Once all of the variables have been considered, control passes to act 424. In some implementations, the processor may incrementally consider each variable in the variable set in turn, and evaluate all of the variables in the variable set as having been considered when the last variable in the variable set has been considered.

[0113] At 424, the processor increments the progression parameter. For example, if the progression parameter is an inverse temperature and was initialized at t0 in act 406, the progression parameter may be incremented to t1 during the first iteration at 424. In some implementations, t1 may be a temperature lower than t0. Over successive iterations, the temperature may be further decreased until a final temperature is reached.

[0114] At 426, penalty parameters for each constraint may optionally be adjusted. In some implementations, the penalty parameters may be adjusted based on a change in energy of the constraint function defined by the i-th variable and the progress parameter. In other implementations, the penalty parameters may be adjusted as described in methods 800 and 900, discussed in more detail below.

[0115] At 428, the processor evaluates one or more termination criteria. In some implementations, the termination criteria may be the value of a progress parameter. In other implementations, the termination criteria may include a number of iterations, an amount of time, an average energy change threshold between updates, a measure of the quality of the current value of a variable, or other metrics known in the art. If the termination criteria are not met, the method continues with act 412. Acts 412-428 are iteratively repeated until the termination criteria are met. In some implementations, the set of variables sampled in act 420 may be received as input by an optimization algorithm in act 412, and the samples may be modified by the optimization algorithm before passing to act 414. As described above, incrementing a stage of the optimization algorithm with respect to the objective function may include incrementing a simulated annealing algorithm or a parallel tempering algorithm. In other implementations, the optimization algorithm may be an MCMC algorithm or a greedy algorithm. Other optimization algorithms include branch and bound algorithms, quantum annealing, quantum approximate optimization algorithms (QAOA) or other noisy intermediate-scale quantum (NISQ) algorithms, quantum-implemented fault-tolerant optimization methods, or other quantum optimization algorithms.

[0116] If one or more termination criteria are met, control passes to 430 where the solution is output. At 430, method 400 ends, e.g., until called again. The solution output in act 430 can be passed to other algorithms, such as a quantum annealing algorithm.

[0117] In some implementations, after outputting the solution in act 430, method 400 can begin again with new samples that were initialized in act 404. In some implementations, method 400 can be run multiple times in parallel starting from different initialized samples or randomly generated sets of samples. In some implementations, solutions can be paired and a binary problem can be constructed to evaluate the set of solutions using cross-Boltzmann updates as described in U.S. Provisional Patent Application No. 62 / 951,749.

[0118] The methods described herein use an optimization algorithm to provide candidate solutions, and in some implementations, simulated annealing may be used. Simulated annealing probabilistically determines whether to accept a candidate solution based on both the change in the energy of the candidate solution and the stage of simulated annealing. In early stages of simulated annealing, candidate solutions that do not provide lower energy are more likely to be accepted, while in later stages of simulated annealing, candidate solutions that do not provide lower energy are more likely to be rejected. This can be thought of as evolving a system from a high temperature to a low temperature. The probability of accepting a candidate solution depends on a probability acceptance function that depends on both the energy of the current solution and the candidate solution, and a time-varying parameter often represented or referred to as temperature. Another alternative is parallel tempering (also known as replica-exchange Markov chain Monte Carlo). Similar to simulated annealing, parallel tempering starts with a random initial solution and exchanges candidate solutions based on temperature and energy. Other optimization algorithms include MCMC algorithms, greedy algorithms, branch and bound algorithms, quantum annealing, quantum approximate optimization algorithms (QAOA) or other noisy intermediate-scale quantum (NISQ) algorithms, quantum-implemented fault-tolerant optimization methods, or other quantum optimization algorithms.

[0119] The method described herein converts the constraint function into part of the objective function, allowing the solver to implicitly consider the constraint function and penalize violations without having to determine appropriate penalty values. The method considers the interaction of each variable with all connected variables, as well as the interaction of each variable with each constraint. The method can also adjust the strength of the penalty within the function until each constraint is satisfied. For each variable, the method can perform variable type-specific updates to allow for the inclusion of different variable types.

[0120] 5 is a flow diagram of an example method 500 of operation of a computing system for finding a solution to an optimization problem or other problem having binary variables, which may be expressed as a constrained quadratic model. It will be understood that similar principles can be applied to other types of variables discussed herein by converting other variable types to binary variables or by addressing energy bias and sampling distributions as discussed elsewhere herein. Method 500 may be performed on a hybrid computing system including at least one digital or classical processor and a quantum processor, such as hybrid computing system 100 of FIG. 1, or may be performed by a classical computing system including at least one digital or classical processor.

[0121] Although method 500 includes acts 502-524, one skilled in the art will understand that the number of acts shown is exemplary and that in some implementations certain acts may be omitted, additional acts may be added, and / or the order of acts may be changed.

[0122] The method 500 begins at 502, for example, in response to a call or invocation from another routine.

[0123] At 502, the processor receives a problem definition. The problem definition includes a set of variables, an objective function defined over the set of variables, and one or more constraint functions, each of the constraint functions defined by at least one variable in the set of variables. The problem definition may be received from a user (e.g., by an input device), transmitted from another processor, retrieved from memory, or provided as the output of another process executed by the processor. The objective function may be a quadratic function, and the problem definition may define a quadratic optimization problem. The set of variables may include one or more of continuous variables, discrete variables, binary variables, and integer variables.

[0124] At 504, the processor generates two candidate values ​​for each variable in the objective function. Examples of methods for generating candidate values ​​are described in U.S. Provisional Patent Application No. 62 / 951,749. The candidate values ​​may be generated by an optimization algorithm, may be generated randomly, or a random solution may be selected and a replacement value based on the random solution generated. A solution may be selected randomly from the entire variable space or within a range of the variable space based on known characteristics of the problem definition or other information. For example, if the i-th variable is binary and the current value of the i-th variable is 0, the replacement value for the i-th variable is 1. For non-binary variables, replacement values ​​may be generated by sampling from a known distribution. In some implementations, the processor may use a native sampler, such as a Gibbs sampler, in the space of integer or continuous variables to select replacement values. It will be appreciated that replacement values ​​may be generated by a variety of algorithms and may depend on the variable type and known parameters of the problem.

[0125] At 506, the processor selects the i-th variable from the variable set, which may be, for example, the first variable in the variable set in the first iteration, the second variable in the variable set in the second iteration, etc., until all of the variables in the variable set have been selected.

[0126] At 508, the processor calculates the difference in the energy of the objective function based on the two candidate values ​​of the i-th variable generated at 504.

number

[0127] At 510, the processor calculates the difference in energy for each constraint function between two candidate values ​​of the i-th variable. In some implementations, the energy difference is calculated as Δ E =Σ c (│δ c,1 │ n -│δ c,0 │n ) can be given by

[0128] At 512, the processor samples the value of the i-th variable based on the difference in energy values ​​of the objective function and each constraint function defined by the i-th variable. As described above, the sampled value may be sampled from a weighted distribution based on the total energy change and a progress parameter, such as a Bernoulli distribution for binary variables.

[0129] At 514, the processor evaluates whether all of the variables in the variable set have been considered. In some implementations, the processor may incrementally consider each variable in the variable set in turn, and evaluate all of the variables in the variable set as having been considered when the last variable in the variable set has been considered. If all of the variables have not been considered, control returns to act 506, where the next variable is selected. Once all of the variables have been considered, control passes to act 516.

[0130] At 516, the processor stores the energy value for each constraint function based on the sampled values.

[0131] At 518, the processor optionally adjusts the penalty value applied to the energy calculation with respect to changes in the energy values ​​of the constraint functions. As discussed in more detail below, the penalty value can be increased if the constraint is not met, while the penalty value can be decreased if the constraint is met. Implementations of automatic adjustment of penalty values ​​are discussed below with respect to methods 800 and 900. In other implementations, the penalty value can be adjusted by a fixed value in response to a constraint being met or violated.

[0132] At 520, the processor evaluates one or more termination criteria. The termination criteria may include the value of a progress parameter, the number of iterations, the amount of time, an average energy change threshold between updates, a measure of the quality of the current value of a variable, or other metrics as known in the art. If the termination criteria are not met, the method continues with act 522. Acts 504-522 are repeated iteratively until the termination criteria are met. If the termination criteria are met, control passes to 524, where the solution is output. At 524, the method 500 ends, for example, until called again.

[0133] At 522, the processor increments the optimization algorithm for the objective function to provide a new set of acceptable values ​​for the set of variables. The optimization algorithm can start from the sampled values ​​of act 512 and generate alternative solutions to provide two candidate values ​​for each variable. The optimization algorithm can include simulated annealing, parallel tempering, Markov chain Monte Carlo techniques, branch-and-bound algorithms, and greedy algorithms, which can be executed by a classical computer. The optimization algorithm can also include algorithms executed by a quantum computer, such as quantum annealing, quantum-approximate optimization algorithms (QAOA) or other noisy intermediate-scale quantum (NISQ) algorithms, quantum-implemented fault-tolerant optimization methods, or other quantum optimization algorithms. The quantum computer can include a quantum annealing processor or a gate-model-based processor. Examples of optimization algorithms are described in U.S. Provisional Patent Application No. 62 / 951,749 and U.S. Patent Application Publication No. 2020 / 0234172.

[0134] 6 is a flow diagram of an example method 600 of operation of a hybrid computing system for finding a solution to an optimization problem or other problem that can be expressed as a constrained quadratic model. Method 600 can be performed on a hybrid computing system that includes at least one digital or classical processor and a quantum processor, such as hybrid computing system 100 of FIG. 1.

[0135] Although method 600 includes acts 602-608, one skilled in the art will understand that the number of acts shown is exemplary and that in some implementations, certain acts may be omitted, additional acts may be added, and / or the order of acts may be changed.

[0136] The method 600 begins at 602, for example, in response to a call from another routine.

[0137] At 602, the processor receives a constrained quadratic model (CQM) problem (also referred to herein as a constrained quadratic optimization problem), the constrained quadratic optimization problem having a set of variables, an objective function defined over the set of variables, one or more constraint functions, each of the constraint functions defined by at least one variable of the set of variables, and progression parameters for optimization, the progression parameters including a set of values ​​that increment between an initial value and a final value.

[0138] At 604, the processor generates a sample solution to the CQM problem using an optimization algorithm 604a and a sampling algorithm 604b. These methods may be similar to methods 300, 400, and 500 described herein. In some implementations, act 604 may include iteratively sampling a sample set of values ​​of the set of variables from the optimization algorithm until a final value of the progress parameter is reached; updating the sample set of values ​​using an update algorithm, where for each variable of the set of variables, determining an objective energy change based on the sample values ​​of the variable and one or more terms of an objective function that includes the variable; determining a constraint energy change based on the sample values ​​of the variable and each of the constraint functions defined by the variable; determining a total energy change based on the objective energy change and the constraint energy change; determining a variable type for the variable; selecting a sampling distribution based on the variable type; and updating the sample set of values ​​using the update algorithm by sampling updated values ​​of the variable from the sampling distribution based on the total energy change and the progress parameter; then returning the updated samples, where the updated samples include updated values ​​for each variable of the set of variables; and incrementing the progress parameter. In other implementations, act 604 may include any other method described herein, such as method 1400 of FIG.

[0139] At 606, one or more final samples generated after reaching a final value of the progress parameter or another termination criterion described above may be transmitted to a quantum processor, and the processor may instruct the quantum processor to refine the samples. In some implementations, transmitting one or more final samples to the quantum processor includes transmitting (at 606a) pairs of samples to the quantum processor, and instructing (at 606b and 606c) the quantum processor to refine the samples by performing quantum annealing to select between the samples.

[0140] At 608, the processor outputs solutions to the CQM problem. In some implementations, these solutions may be returned to a user (e.g., via an output device). In other implementations, the solutions may be passed to another algorithm for refinement, or may be returned to act 604 as a sample set of values ​​as input to an optimization algorithm. At 608, method 600 ends unless it is repeated, or for example, called again.

[0141] As described above with respect to method 600, in some implementations, the obtained samples may be further refined by a quantum computer, for example, as part of a hybrid algorithm using cluster shrinkage. The cluster shrinkage hybrid algorithm is disclosed in more detail in U.S. Patent Application Publication No. 2020 / 0234172. In some implementations, solutions may be paired, and a binary problem may be constructed to evaluate a sequence of solutions using cross-Boltzmann updates, as described in U.S. Provisional Patent Application No. 62 / 951,749. The initial solution may be improved using a hybrid computing system including at least one classical or digital computer and a quantum computer, such as hybrid computing system 100 of FIG. 1.

[0142] Cross Boltzmann updates can find or yield two possible solutions along with two different values ​​for a given variable. A new binary variable is defined to choose from two different values. This can then be reduced to an optimizable quadratic unconstrained binary optimization (QUBO) problem. In some implementations, the resulting QUBO problem can be sent to a quantum processor, and quantum annealing can be used to generate a solution. A classical or digital computer can use cross Boltzmann updates to construct a binary problem that a quantum computer can solve. Cross Boltzmann updates use the energy associated with updating two variables individually and dependently to yield a Hamiltonian that expresses the decision of whether to update those variables. The digital processor can select two candidate values ​​and send this update to the quantum processor as an optimization problem.

[0143] Define the QUBO problem using a super-diagonal matrix Q, an NxN upper triangular matrix of real weights, and a vector x of binary variables to minimize the function:

number

[0144] This is more succinctly

number

[0145] In scalar notation, the objective function expressed as a QUBO is:

number

[0146] inequality As mentioned above, the contribution of the constraint to the change in energy can be expressed as:

number

[0147] For example, in the inequality expressed as:

number

[0148] δ c Θ(δ c ), where Θ(δ) is the function shown in graph 704, penalizing penalizes violations only if the value of the constraint is positive.

[0149] Lagrangian parameters As mentioned above, the change in energy contributed by both the objective function and the constraints is

number

[0150] The Lagrangian multipliers can be adjusted during optimization by evaluating violations of the constraint functions: if a given constraint is violated, the associated Lagrangian multiplier is increased; if a given constraint is satisfied, the associated Lagrangian multiplier is decreased.

[0151] As per the implementation above, the CQM problem can be expressed as:

number

number

[0152] In simulated annealing, the magnitude of the change in the objective value is evaluated for each variable where the simulated state of the variable changes from its current value.

number

number

number

number

number

number

[0153]

number

number

number

number

[0154]

number

[0155] After finding the first feasible solution, a global decay is used to stabilize the penalty value. The parameters are the decay coefficient df:

number

[0156]

number

[0157] 8 and 9 are flow diagrams of example methods 800 and 900 for adjusting penalty parameters to guide the search space toward feasibility, which can advantageously improve the performance of a processor-based system executing the methods 800 and 900. The methods 800 and 900 can be implemented as a subprocess of simulated annealing or as part of an optimization algorithm. As described below, the method 900 begins with initializing Lagrangian parameters and multipliers, specifically the Lagrangian initial values ​​μ, α, and γ. In addition, counters tf and tf are initialized to count the number of iterations during which the search fails to reduce the infeasibility (total violations) and develop the best solution, respectively. tf is the number of iterations during which the search remains in the infeasible region without reducing the infeasibility (total violations), while tf is the number of iterations during which the search remains in the feasible region without improving the best feasible solution. As described above, t is the number of iterations the algorithm waits before increasing α and γ.

[0158] During the simulated annealing update method, given a state, its current objective, and the energy of the left-hand side of the constraints, the system checks the feasibility of the solution, its distance from the constraint bindings, and their violations. The system then adjusts the search towards feasibility or infeasibility according to the improvement of the best solution and the sum of violations. Finally, the Lagrangian multiplier is adjusted based on the new values ​​of α, γ, violation / constraint bindings, and the Lagrangian multiplier from the previous iteration.

[0159] In some implementations, the Lagrangian parameters may be varied directly as the simulated annealing progresses, along with a progression parameter such as inverse temperature. At the beginning of the simulated annealing, the Lagrangian parameters may be small and the penalties assigned to violations may be small. Towards the end of the simulated annealing, the penalties assigned to violations may be increased so that the constraints may be satisfied in the final solution.

[0160] 8 is a flow diagram of an example of a method 800 for guiding a search space to feasibility in an optimization algorithm, which may advantageously improve the performance of a processor-based system when performing the optimization. Method 800 may be performed on a hybrid computing system including at least one digital or classical processor and a quantum processor, such as hybrid computing system 100 of FIG. 1, or may be performed by a classical computing system including at least one digital or classical processor.

[0161] Although method 800 includes acts 802-822, those skilled in the art will understand that the number of acts shown is exemplary and that in some implementations, certain acts may be omitted, additional acts may be added, and / or the order of acts may be changed.

[0162] The method 800 begins at 802, for example, in response to a call or invocation from another routine.

[0163] At 802, a processor receives samples from an optimization algorithm, such as a simulated annealing algorithm as described above. In some implementations, method 800 may be invoked by any one of methods 300, 400, 500, and 600 described above, or method 1400 discussed below. Method 800 may be used to adjust penalty parameters, such as in act 426 of method 400 or act 518 of method 500.

[0164] At 804, the processor determines energy values ​​or biases of one or more constraint functions. In some implementations, the processor may evaluate the constraint functions for a given sample. In other implementations, such as the methods described above, the energy values ​​of one or more constraint functions may be calculated by a separate function and passed to act 804, such as by act 418 of method 400 or act 516 of method 500.

[0165] At 806, the processor evaluates the feasibility of the sample. This may be determined, for example, by received values ​​for the energy values ​​of one or more constraint functions. In some implementations, a function may be provided that zeros out the energy values ​​of feasible samples, as discussed above with respect to FIG. 7. If a constraint is violated, the sample is deemed infeasible.

[0166] At 808, if the sample is not executable, control passes to 810; if the sample is executable, control passes to 812.

[0167] At 810, the processor increases the penalty value of the constraint that is violated for the infeasible sample. In some implementations, the penalty value may be incremented by a fixed value for each iteration that returns a infeasible sample. Because only feasible samples are desired, it is beneficial to increment the penalty value until the constraint functions are sufficiently weighted so that the algorithm returns only feasible solutions.

[0168] At 812, the processor decreases the penalty value for feasible samples. The penalty value may be decreased if a feasible sample is returned because overweighting a constraint can lead to a feasible solution that is less than optimal. As discussed below with respect to method 900, a decay factor may be used to cause the penalty value to converge to a lower value that returns a feasible solution. In some implementations, the penalty value may be decreased by a fixed value for each iteration that returns a feasible sample until a non-feasible sample is returned. The penalty value may then be increased and maintained for the remaining iterations until a feasible sample is returned.

[0169] At 814, the processor optionally determines whether the violation has decreased compared to the previous sample. For example, this can be done by comparing constraint violation functions as described above. If the total energy penalty contributed by the constraint function increases, the violation is increasing. If the total energy penalty contributed by the constraint function decreases, the violation is decreasing. In some implementations, the violation may remain the same and therefore not decrease. If the violation has decreased, control passes to 818; otherwise, to 822.

[0170] At 818, the processor increases the initial adjuster value. For example, as discussed above, α can be increased exponentially by a constant factor (e.g., 1.2).

[0171] At 816, the processor optionally determines whether the current best solution has improved compared to the previous sample. If the current best solution has improved, control passes to 822; otherwise, control passes to 820. For example, this can be done by comparing the energy of the objective function in the previous iteration to the energy of the objective function.

[0172] At 820, the processor increases the initial adjuster value. For example, as discussed above, γ can be increased exponentially by a constant factor (e.g., 1.2).

[0173] At 822, the processor outputs the penalty value to return to the optimization algorithm. The method 800 may then end, for example, until called again.

[0174] Figure 9 is a flow diagram of an example implementation of method 800 for guiding a search space for feasibility in an optimization algorithm in the form of a method 900 that can advantageously improve the performance of a processor-based system when performing the optimization. Figure 9 illustrates how the Lagrangian multipliers can be automatically adjusted during a simulated annealing procedure for solving a CQM problem. Figure 9 shows the main simulated annealing procedure 902, along with acts 904-922 that illustrate the automatic adjustment procedure of the multipliers.

[0175] At 904, the system initializes the Lagrangian parameters. In the implementation of FIG. 9, the initial adjuster α 0 and γ 0 are set to 2 and 0.05, respectively, and the Lagrangian parameters

number

[0176] At 906, the system receives the energy of the constraint functions and the current states of the variables from an optimization or update algorithm. For example, sample values ​​generated by one of methods 300, 400, 500, 600, or 1400 may be passed to the system along with the energy for the constraint functions found by one of those methods. In some implementations, the system may receive the energy of the constraint functions and the current states of the variables from a simulated annealing algorithm.

[0177] At 908, the system evaluates the damping factors, violations, and bindings. The damping factors may be provided as initial parameters or may be adjusted depending on the variation between feasible and infeasible values ​​in previous iterations. The violations may be evaluated from constraint functions such as those described above. For example, the constraints may be evaluated and multiplied by the functions discussed with respect to FIG. 7. The bindings may be calculated as the distance from the binding for each constraint where the constraint is not violated and / or the total distance from the binding across all constraints.

[0178] At 910, the system determines whether the current sample is feasible, i.e., whether any of the constraints are violated. This may be determined based on the energy of the constraints, as described above. If one or more of the constraints are violated, the method proceeds to act 912, whereas if the sample is feasible (i.e., no constraints are violated), the method proceeds to act 914.

[0179] At 912, the system determines whether the violation has been reduced. As described above, this can be determined based on the energy of the constraint. If the violation has been reduced, control passes to 916 and the counter values ​​are unchanged (tinf and tf remain zero). If the violation has not been reduced, control passes to 918.

[0180] In step 918, the counter tinf is incremented by 1 and compared with t. If tinf is equal to or greater than t, the initial adjuster α 0Increase the value of the initial adjuster α (by some constant factor, such as an exponent of 1.2). 0 If the value of is incremented, tinf is reset to zero. Control then passes to 916.

[0181] At 914, the system determines whether the best solution has improved, i.e., whether the energy of the objective function has decreased or whether a more optimal solution has been found. If the solution has improved, control proceeds to 916 and the counter values ​​are unchanged (tinf and tf remain zero). If the best solution has not improved, control proceeds to 920.

[0182] In step 920, the counter tf is incremented by 1 and compared with t. If tf is equal to or greater than t, the initial adjuster γ 0 Increase the value of the initial adjuster γ (by some constant factor, such as an exponent of 1.2). 0 If the value of is incremented, tf is reset to zero. Control is then passed to 916.

[0183] At 916, the system is

number

number

number

number

[0184] At 922, the system calculates the updated Lagrange multipliers for each constraint based on whether the constraint is satisfied, as well as the Lagrange multipliers discussed above.

number

number

number

[0185] The updated Lagrange multipliers are returned to 902 and can be used in other algorithms, such as being returned to a stage of the simulated annealing algorithm.

[0186] Pseudocode example In some implementations, the problem to be solved may involve only binary variables. One implementation of such a binary problem with binary variables x1 and x2 is discussed below. It will be understood that this example problem represents only one possible implementation, and other implementations may surpass what is discussed below.

[0187] The processor receives a constrained quadratic model (CQM) problem having an objective function, such as f(x1, x2), and constraint functions, such as g(x1, x2) and h(x1). The constraint functions may include equalities or inequalities.

[0188] The processor receives or initializes samples. The samples may be generated by another algorithm, received as input, sampled from a distribution, or generated randomly. In this example, a random sample is generated.

number

[0189] The processor receives or initializes a progression parameter. In some implementations, the progression parameter can be an inverse temperature. The processor can initialize the inverse temperature as t, for example. For example, in some implementations, simulated annealing can be used to perform the optimization. The temperature of the system can be initially set to a high temperature and incremented to a low temperature. This can have the effect of accepting high-energy solutions early in the optimization, reducing the probability of getting stuck in a local minimum, and preventing high-energy solutions from being accepted near the end of the optimization.

[0190] The processor may optionally receive or initialize a penalty parameter, e.g., a function of the progress parameter and the magnitude of the energy violation, e.g., λ c (t,δ c ) and the Lagrangian parameters can be initialized as discussed in more detail above.

[0191] The processor calculates the current value of each constraint based on the received or initialized samples.

number

[0192] The processor then considers each variable in turn. The variables have current values ​​given by the initialized samples, and the processor can determine alternative values ​​based on the properties of the variables or constraints on the variables, or can receive alternative values ​​as input or from another algorithm. For example, if x1 is a binary variable,

number

[0193] For each constraint that includes x1, the processor determines an energy value for the constraint based on the pairwise terms of the initially calculated constraint values ​​g(0,1), h(0), and g(1,1), the alternative values ​​of x1 that include x1. In the given example, the constraint energy is δ x1,1 = g(1,1) + h(1) and δ x1,0 = g(0,1) + h(0). These terms can then be calculated by the processor as Δ E =Δ E,x1 +λ t0,δx1,1 δ x1,1 2 Θ(δ x1,1 )-λ t0,δx1,0 δ x1,0 2 Θ(δ x1,0 ) to the current energy value, where Θ(δ) is a function due to the inequality constraints, as discussed in more detail above with respect to Figure 7. If no inequality constraints are violated, Θ(δ) is set to zero, so that the term does not contribute to the energy.

[0194] Next, Δ EA sample value of the variable x1 is generated by the processor based on t0 and t0. The sampling can vary depending on the variable type, for example, sampling from a distribution or generating the sample in another way. In some implementations, if the sampled value does not violate any constraints, the sampled value is accepted, while if the sampled value does violate a constraint, the sampled value is accepted with a probability that depends on t0. This value is then calculated as x 1,t0 , e.g. x 1,t0 =1.

[0195] The processor then considers x2 and repeats the process to find x 2,t0 Sample the value of

[0196] Once all of the variables have been considered and a new sample value has been generated, the processor increments the progress parameter.

number

number

number

number

number

[0197] At the final inversion temperature, the resulting sample

number

number

number

[0198] Here is an example of pseudocode for the simple binary problem discussed above: Receive the objective function f(x1,x2), constraints g(x1,x2), and h(x1) Sample initialization

number

number

[0199] The methods and systems described herein advantageously enable computationally efficient solving of optimization problems with constraints by directly processing the constraints, as opposed to converting the constraints into part of the objective function. This may advantageously enable a processor to return a feasible solution to an optimization problem in a time-efficient and scalable manner. This may also enable the processor to generate multiple solutions in parallel. The contribution of the constraints may be implicitly handled by an update process that considers each constraint. Processing the constraint functions directly and deriving their energy values ​​through computations may increase the efficiency of the solution. The methods and systems described herein may also advantageously enable automatic and dynamic adjustment of penalty weights for constraint functions, allowing the search to efficiently guide the search toward feasibility while returning better solutions. Many real-world problems, such as those solved in industry, are problems with constraints on feasible solutions. The methods and systems described herein may also enable the solution of problems with different types of variable sets, such as a mixture of binary and discrete variables. Sampling may be performed based on a determination of the variable type and taking into account the current energy values ​​and penalties. In some implementations, the methods and systems described herein can be used in combination with hybrid quantum computing techniques to yield improved solutions.

[0200] Integer variables As mentioned above, constrained problems can be defined over a variety of variable types, such as binary, discrete, integer, and continuous. The variable types can determine the sampling distributions selected by the processor to sample update values ​​of those variables. In some implementations where the variable types are binary, the selected sampling distribution can be a Bernoulli distribution. In other implementations where the variable types are discrete, the selected sampling distribution can be a softmax distribution.

[0201] When the variable type is integer (or continuous), there may not be a single well-defined distribution to choose from. Many optimization problems involve integer variables and constraints on the values ​​those integer variables can take. To support integer variables in a constrained problem solver as described above, it may be beneficial to provide "native" support for integer variables. Updating integer variables taking into account their constraints and Lagrangian parameters is discussed below.

[0202] One technique for implementing integer updates involves temporarily relaxing integer variables to continuous variables, computing the objective and gradient via Lagrangian relaxation of constraints involving that variable, proposing an update via that gradient, and probabilistically accepting or rejecting it, similar to the process described above. However, this type of update for integer variables can result in high algorithmic complexity (on the order of the number of adjacent constraints) when the constraints are complex and / or there are many overlapping constraints. Alternatively, integer variables defined over a bounded range of size R can be implemented similarly to binary variables by converting the integer variable to a binary variable of approximately log(R). This can then be updated by selecting a Bernoulli distribution as the sampling distribution and proceeding with the update as described above for binary variables. However, this type of update requires choosing the best way to convert the integer variable to a binary variable. Some examples include base-2 representation or Gray code. For a binary variable of log(R), sampling can take exponential time, resulting in a running time on the order of R.

[0203] It may be beneficial to perform integer updates by examining the neighborhood constraints of a given integer variable and dividing the range of the integer variable into a series of segments, each segment having the same form of constraint penalty. A segment may be selected randomly based on the integral of the objective and constraint penalty within that segment. A new value for the integer may be sampled from within that segment. The following method is discussed in the context of integer variables, but can also be extended for use with continuous variables, as discussed below.

[0204] This process takes O(M log(M)+Mt) time, where M is the number of constraints and t is the execution time required to compute the partition function or samples on the segment. If the integer variable has only linear terms (or multiple linear terms) in the objective function, then t = O(1), and the update written below is expected to take O(M log(M)) time. If the integer variable has quadratic terms in the objective, then t is expected to be O(log(r)), where r is the range of the segment, which may be O(1) in some cases. This is beneficially expected to be faster than the other methods mentioned above, such as relaxing or converting integer variables to continuous variables.

[0205] When sampling values ​​of integer variables, Gibbs sampling can be performed, i.e., sampling with a conditional probability given the current state, as will be explained in more detail below. As mentioned above, in the context of solving CQM problems, both the objective function and the constraint functions are considered when sampling new values. Therefore, to update a given integer variable, sampling can be performed with a conditional probability given the current state or currently accepted value of the variable. The conditional probability of the objective function can be calculated analytically as described below.

[0206] Although quadratic functions involve squared variables, it will be understood that some of the equations in the implementations described herein may not work for squared integer variables. In those implementations, a new variable can be defined. For example, if an integer variable x1 is squared in the objective function or constraints, a new variable x2 can be defined and the square can be replaced with x1x2, along with the additional constraint x2=x1. It will be understood that this substitution may be required in equations used above, such as when variables are written with i and j subscripts, where i=j. It will be understood that other substitutions can be made as well to accommodate squared variables. In other implementations, the following equations can be modified. An alternative implementation for accommodating squared (quadratic) variables without adding an extra variable is also discussed below using slice sampling.

[0207] Constraints are considered by calculating a binding value. Assuming the effective linear bias of a variable with respect to its adjacent constraints is greater than zero, there are three possible outcomes: the binding value may be less than the variable's lower bound, the binding value may be between the variable's upper and lower bounds, or the binding value may be greater than the variable's upper bound. If the binding value is greater than the variable's upper bound, the constraint is satisfied for all possible values ​​of the variable, and the partition function is given by a conditional partition function that includes only the conditional objective function. If the binding value is less than the variable's lower bound, all possible values ​​of the variable violate the constraint and a penalty must be imposed. The opposite is true if the variable's effective linear bias is less than zero. As mentioned above, the penalty can be applied as a constant or can vary depending on the magnitude of the violation and the stage of the optimization algorithm.

[0208] When the binding value is between the upper and lower bounds of the variable, the domain can be divided into two segments, as shown in Figure 10. As above, Figure 10 assumes that the effective linear bias is greater than zero. When the effective linear bias is less than zero, the result is inverted. The first segment is defined by the lower bound of the variable and an integer less than or equal to the binding value, and the second segment is defined by the upper bound and an integer greater than or equal to the binding value. The partition function can be calculated as the sum of the partition functions of these two segments. The first segment can use the same partition function when the binding value is greater than the upper bound of the variable, i.e., the partition function conditional on the objective function. The second segment can use the same partition function when the binding value is less than the lower bound of the variable, with a penalty applied. When the effective linear bias of the variable is less than zero, the situation is the same as above, but the cases are swapped. This can be used to generalize to multiple constraints.

[0209] Given a function of the objective and constraints, new values ​​can be sampled by Gibbs sampling based on the conditional probability given the current state.

[0210] As mentioned above, a constrained quadratic model (CQM) is defined by an objective function and a set of constraints. A CQM is a set of variables x, which can be, for example, a combination of binary variables, integer variables, discrete variables, and continuous variables, or other types of variables. i In some implementations, one or more of the variables are integer variables.

[0211] In one implementation, the objective function, e.g., f(x)=Σ i a i x i +Σ i<j w i,j x i x j is defined. The form g c (x)=Σ i a i,c x i +Σ i<j w i,j,c xi x j +C c A set of constraints on M is provided: ≦, ≧, or =0. i,c or w for any j i,j,c If is non-zero, then the constraint g c (x) is the given variable x i Integer variables in a CQM are considered to be "adjacent" with respect to the variable x k These bounds can alternatively be expressed as constraints.

[0212] Given an integer variable x k To update x, Gibbs sampling is performed based on the conditional probability given the current state of all other variables. Consider the objective function, and the integer variable x k When updating, the conditional partition function is

number

[0213] When inserted into the partition function, the first term on the right-hand side of the equation (exp[-βA]) becomes an irrelevant prefactor, and the partition function can be taken as:

number

[0214] δ eff >0, the partition function is

number

number

[0215] δ eff <0, the partition function can be calculated analytically using the same series identity as above with the substitution m=-n, yielding

number

[0216] Considering the constraints, all constraints are expressed as inequality constraints, and equality constraints are implemented as a combination of two inequalities, with the variable x k Any constraint c adjacent to

number

number

number

[0217]

number

number

number

number

number

number

[0218] Considering case 3, the constraint is x k is satisfied for all possible values ​​of , and there is no need to constrain values. The conditional partition function is given by equation (1).

[0219] Considering case 1, the variable x k All possible values ​​of violate the constraint and should therefore be penalized. The partition function is

number

number

number

[0220] Considering case 2, the variable x is k The domain of is divided into two segments:

number

number

number

[0221] The partition function can be calculated as the sum of the partition functions in the two segments. In the first segment, we can use the results from Case 3, and in the second segment, we can use the results from Case 1, both adapted to the domain of the segments. Then the partition function becomes:

number

[0222]

number

[0223] Generalizing the partition function to the case of multiple constraints, an ordered set of binding values

number

number

number

[0224] This constraint set is called C s and the partition function of the segment s is

number

[0225] This general case is illustrated in FIG.

[0226] For the L0 penalty with α=0, the partition function becomes:

number

[0227] The L1 penalty allows us to incorporate the linear dependence of the violation on the objective function bias, giving us an analytical expression for the partition function, which gives us:

number

[0228] Change the prefactor,

number

[0229] Therefore, the partition functions of L0 and L1 are given by the linear penalty LP with an offset ω and a slope σ for each segment s. s of

number

number

[0230] Inverse Transform Sampling Sampling from within the selected segment can be done using various techniques, for example, using inverse transform sampling. Inverse transform sampling is described below for integer and continuous variables. Let x be the segment x∈[x L ,x U ]. See the discussion above for x L and x U is the segment boundary x s , x s+1 The partition function of a segment is generally

number

number

number

[0231] Segment [x] with distribution exp[-λx] L ,x U To sample a continuous variable in [0,1], we can randomly sample a variable y that is uniformly distributed in [0,1] and then use equation (9).

[0232] The integer case can be sampled in a similar way, and the cumulative distribution function is

number

number

[0233] The above techniques for integer variables will now be discussed with reference to Figure 3 and method 300. As noted above, method 300 is one implementation of how a computing system operates to update samples in an optimization algorithm to improve convergence to feasibility. It will be appreciated that the above techniques for integer variables are equally applicable to any of the methods discussed herein.

[0234] At 302, a processor receives a problem definition having a variable set, an objective function defined over the variable set, and one or more constraint functions, each of the constraint functions defined by at least one variable of the variable set, in this implementation, at least one variable of the variable set is an integer variable.

[0235] At 304, sample values ​​for the set of variables are received by the processor, which may be based on some known characteristic of the system, may be provided by another algorithm, may be randomly selected, or may be provided by other techniques recognized by those skilled in the art. Values ​​for process parameters, such as temperature in the case of simulated annealing, are also received.

[0236] At 306, the processor selects a variable from the variable set. In this example, only integer variables are discussed; however, this method can also be applied to continuous variables. As discussed in more detail above, other variable types can also be updated. It will be appreciated that the order of acts 308-318 may be changed from the order shown in FIG. 3, as discussed in more detail below.

[0237] The variable type is determined at 308. In this example, the variable type is determined to be integer.

[0238] At 310, a sampling distribution for the integer variable is selected. As noted above, there is no single, well-defined distribution that can be selected for an integer variable. Instead, Gibbs sampling, which refers to a technique for sampling from a conditional distribution, is performed. The sampling distribution is selected to be the conditional distribution given by the conditional partition function.

[0239] At 312, a target energy change bias is determined based on the sample values ​​of the variables and one or more terms of the objective function that includes the variables. As mentioned above, for integer variables this is the variable x k is calculated as the effective linear bias of eff =a k +Σ i≠k w ik x i is given by

[0240] At 314, a constraint energy change based on each of the sample values ​​of the variables and the constraint functions defined by the variables is determined.

number

[0241] At 318, update values ​​for the variables are sampled based on the objective energy change and the constraint energy change, as described above. Update values ​​may also be selected based on the current value of each constraint and the type of each constraint, as described above. Specifically, for a given variable

number

number

[0242] At 320, it is determined whether all of the variables have been considered. If there are other integer variables, the above process for sampling is repeated. If there are other variable types, sampling is performed as above.

[0243] At 322, an updated sample is returned, the updated sample including an updated value for each variable in the variable set.

[0244] continuous variables As discussed above, constrained problems can be defined over a variety of variable types, such as binary, discrete, integer, and continuous. The variable type can determine the sampling distribution selected by the processor. In the discussion above for integer variables, there may not be a single, well-defined sampling distribution. Similarly, continuous variables may not have a single, well-defined sampling distribution. A continuous variable can take on any two specific real values, and can also take on all real values ​​between those two values. For example, a continuous variable can take on any value across a range of real numbers. Many types of optimization problems involve decision variables represented by continuous variables. Some examples of such continuous variables include weight, volume, pressure, temperature, speed, flow rate, and elevation.

[0245] One way to address continuous variables in optimization problems is to use discretization of continuous variables to represent them using a set of integer or binary variables and use the techniques described above for the redefined variable type. However, discretization of continuous variables may not provide a computationally efficient solution. For example, converting continuous variables to integer variables may require an infinite upper bound to model the continuous variable with arbitrary precision. Integer variables with infinite bounds may be impractical to implement in many applications. Using binary variables may similarly significantly increase the problem size and lead to inefficient solutions.

[0246] Inverse transform sampling of continuous variables As mentioned above, the inverse transform sampling method for integer variables can be easily adapted to continuous variables, the main difference being how the partition function is calculated and how the segments are divided. The partition function for a continuous variable can be written as:

number

[0247] The prefactor exp{-βA} can be neglected, and δ eff >0 or δ eff< 0, we can analytically calculate the partition function without having to distinguish between the two cases.

[0248] variable x k Constraints adjacent to g c If you consider this, it will be the binding value

number

number

number

number

[0249] Comparing the above equation with its integer equivalent demonstrates how sums can be replaced by integrals and floor and ceiling functions for the binding values ​​eliminated. When considering multiple constraints, the partition function becomes

number

[0250] Even in this case, a uniform approach can be used to analytically calculate the partition functions for the L0 and L1 penalties.

number

number

[0251] The methods described above use inverse transform sampling or similar sampling methods across segments and cannot be used for variables (integer or continuous) with quadratic cases.

[0252] The cumulative distribution function can be analytically inverted only if it has a linear (x) term or a bilinear (x × y) term. 2 For terms like ∑(x) ...

[0253] Similar to the implementation described above, consider the case where we have a constrained quadratic model (CQM) with binary, integer, and continuous variables. Let B ⊆ {1,...,n} denote the set of indices of the binary variables, i.e., if i∈B then x i ∈{0,1}. Similarly, C⊆{1,...,n} denotes the set of indices of continuous variables, i.e., if i∈C, then

number

number

[0254] CQM assumes the following format: Minimize:Σ i a i x i +Σ i≦j b ij xi x j Constraints:

number

[0255] term Σ i≦j b ij x i x j is a pure quadratic term

number

number

[0256] Variable x with k∈C or k∈I k Consider the unconstrained model where we update {x}. In one implementation using the annealing solver and Gibbs sampling, we i} i≠k Given the value of variable x k The conditional probability distribution function of

number

[0257] b kkIf =0, there is no quadratic term, so sampling of integer or continuous variables can be performed using the method previously described for integer variables. The cumulative distribution function of p(x) = exp(-λx) can be inverted for both integer and continuous variables. However, b kk If ≠ 0, the function can no longer be easily inverted. One implementation of sampling for quadratic terms and mixed variable interactions is described below.

[0258] In the case of quadratic terms and mixed variable interactions, sampling of a given variable can be performed by slice sampling. A general description of slice sampling can be found in Neal, R., Slice Sampling, The Annals of Statistics, 2003, Vol. 31, No. 3, 705-767. Slice sampling is generally a Markov chain Monte Carlo algorithm for drawing random samples by uniformly sampling from the area under the graph of the density function of the variable. For a given value of variable x0, the auxiliary variable y is uniformly sampled between 0 and f(x0). A point (x, y) can be sampled from a line (or slice) in y under the curve of f(x0). The value of x at the sampled point becomes the updated value of variable x. While slice sampling is described below, it will be understood that other implementations can use inverse transform sampling with the aid of Newton-Raphson, bisection, or other methods known in the art.

[0259] In a CQM problem,

number

number

[0260] Let u~RAND(0,1) be a random number between 0 and 1.

number

number

number

number

[0261] There are two cases: x1≧lb[k] and x2≦ub[k] hold, x k In the connected domain x ∈[x1,x2] k The new values ​​of can be uniformly sampled from [x1,x2]. x k ∈[lb[k],x1]∨x k In the disconnected region ∈[x2,ub[k]], one of the two disconnected regions is selected with a probability proportional to x1-lb[k] and [x2,ub[k]]. k The new values ​​of are sampled uniformly from the selected region.

[0262] variable x k are integers, then the two values ​​x1, x2 are replaced by x1, x2 → [x1], [x2] in the connected region, and by x1, x2 → [x1], [x2] in the disconnected region.

[0263] Moving on to constraints, we will now discuss "less than" inequalities. It will be understood that "greater than" equalities can be treated similarly, and that an equality constraint can be treated as a combination of two inequality constraints.

[0264] Consider the case with one constraint:

number

[0265] As above, the variable x (which can be an integer or continuous variable) k When updating, the other variables are fixed. The constraints can be written as:

number

number

number

number

number

[0266] b kk >0, the two solutions to the inequality are

number

number

number

number

number

number

[0267] As discussed in more detail above, particularly with respect to methods 800 and 900, a Lagrangian coefficient (penalty parameter) can be provided to penalize infeasible solutions. The inhibited case can also be written as: P(x k =n)∝exp{-β[an 2 +bn+c]} , where the following equation holds:

number

[0268] variable x k The next value of the variable

number

Number

Number

[0269] Each S i can be made empty (when all x within the segment have P(x)<y), or can be truncated as above. Each of these lower regions has a (non-normalized) probability p which is the width of the region (for continuous variables) or the number of integers within the region (for integer variables) i accompanied by. The lower region is selected with probability

Number

[0270] These principles can be extended when there is no solution to the above equation with all ranges appropriately suppressed, and when there is only one solution

Number

Number

Number

[0271] If there are multiple constraints, the variable x k The domain of x can be divided into a set of segments where some constraints are violated. After sampling slice y, the method can proceed similarly to the implementation above, with each segment divided into subregions and P(x) ≥ y. One of the subregions is then sampled with a probability depending on the full width (or number of integers) of the subregion, and the variable x k The new values ​​of are sampled uniformly from within the selected subregion.

[0272] L2 norm penalty In the above discussion of integer variables, the sampling of constraint functions is described in terms of L0 and L1 norms. The general partition function of a segment s is

number

[0273] The L2 norm for the violated constraints where α=2 does not allow for inverse-transform sampling because the above equation contains quadratic terms. However, as noted above, slice sampling can be used to accommodate the use of the L2 norm when constraints have linear or bilinear terms. For constraints with quadratic terms, the L2 penalty introduces up to fourth-order terms, which is beyond the scope of the following discussion. Additional modifications of the methods described herein may be required to accommodate fourth-order cases, or the L0 and L1 norms may be used in these cases.

[0274] Constraints without quadratic terms, i.e.

number

number

number

number

[0275] The square in the exponential function can be expanded as above: x k To sample new values ​​of {tilde over (x)}, the slice sampling procedure described above can be used.

[0276] Pseudocode: Next, an example implementation of the above method is described.

[0277] Input: CQM, state x on the CQM, current Lagrange multiplier λ for each constraint c c , the current energy value E[c] of each constraint c, the inverse temperature β, and the integer variable x to be updated. k .

[0278] The effective bias function is defined as follows, given the states, variables, and constraints:

number

[0279] In one implementation, the following structure is defined: "QuadraticPenalty" has the properties "quadratic", "linear", and "offset" (all real numbers) that encode the quantities introduced in the linear penalty equation defined above. It is used to describe a quadratic penalty with quadratic and linear coefficients, and "offset" is x. k =0 is the penalty value.

[0280] The addition operation can also be defined as follows: Given two QuadraticPenalty objects QP1 and QP2, define QP=QP1+QP2 as follows: QP.quadratic=QP1.quadratic+QP2.quadratic QP.linear=QP1.linear+QP2.linear QP.offset=QP1.offset+QP2.offset.

[0281] A ConstraintPoint has the properties "value" (an integer), "left" and "right" (both "QuadraticPenalty"), which means that one or more penalties vary as different quadratic functions of x. i Describes a point along the domain of x. k , and "left" and "right" refer to different quadratic penalties (or sums of quadratic penalties) that go to 0 closest to "value". "left" refers to x with a negative "slope". k "right" refers to the penalty coming from the negative direction of the domain, and "right" refers to the penalty going in the positive direction of the domain with a positive "slope". For two "ConstraintPoint" objects CP1 and CP2, addition is defined if and only if CP1.value == CP2.value. If this is true, then set CP=CP1+CP2, CP.value = CP1.value, CP.left = CP1.left + CP2.left, and CP.right = CP1.right + CP2.right.

[0282] A "PartialSegment" has the properties "lower_bound" (an integer) and "quadratic_penalty" ("QuadraticPenalty"), which correspond to the sum of the quadratic penalties applied to the variables at the start of each segment (as defined in the mathematical formulation section above) and in the part of the segment following the lower bound.

[0283] "Segment" has the properties "lower_bound", "upper_bound", and "quadratic_penalty".

[0284] In one implementation, the algorithm proceeds as follows: Subroutine 1: sample_integer_continuous_variable input: Inverse Temperature Beta Current State x The variable to be updated xk Lagrangian List CQM Model Initialize two empty arrays of ConstraintPoint objects called tmp_constraints and constraints. Initialize an empty array of PartialSegment objects called partial_segments Initialize an empty array of Segment objects called segments Obtain quadratic_bias and linear_bias from the objective function IF xk is an integer variable THEN: update_integer=1 ELSE: update_integer=0 END IF num_active_constraints=0 FOR each constraint c adjacent to variable xk DO: lagrangian is the Lagrange parameter tmp_constraints,num_active_constraints=PARSE_CONSTRAINT(c,lagrangian,update_integer,tmp_constraints,xk_lower,xk_upper,num_active_constraints) END FOR SORT tmp_constraints by value, then by left and right penalties with quadratic first and offset last constraints,num_constraints=COLLAPSE_CONSTRAINTS(tmp_constraints,num_active_constraints) partial_segments,num_partial_segments=MAKE_PARTIAL_SEGMENTS(quadratic_bias,linear_bias,lower_bound,upper_bound,collapsed_constraints,num_collapsed_constraints,partial_segments,update_integer) segments,num_segments=MAKE_SEGMENTS(partial_segments,num_partial_segments,lower_bound,upper_bound,segments,update_integer) new_value=SAMPLE_VALUE_FROM_SEGMENTS(segments,num_segments,beta) RETURN new_value

[0285] The next subroutine parses the constraint c, i.e., the binding value depending on whether the constraint is linear or quadratic, respectively.

number

number

[0286] Subroutine 2: PARSE_CONSTRAINT input: constraint c lagrangian update_integer tmp_constraints num_active_constraints lower_bound upper_bound From the constraint c, we obtain the quadratic coefficients b_kk, the effective linear bias delta_k, and the offset A_k. IF b_kk=0 THEN #Linear constraints IF delta_k==0 THEN Skip this constraint END IF xc=-delta_k / A_k if update_integer THEN IF delta_k>0 THEN xc=CEIL(xc) ELSE xc=FLOOR(xc) END IF END IF IF L0 penalty is used THEN QP=QUADRATIC_PENALTY(0,0,lagrangian * abs(delta_k)) ELSE IF L1 penalty is used THEN QP=QUADRATIC_PENALTY(0,lagrangian * delta_k,-lagrangian * delta_k * xc) ELSE IF L2 penalty is used THEN QP=QUADRATIC_PENALTY(lagrangian * delta_k ** 2,-lagrangian * delta_k * A_k,lagrangian * A_k ** 2) END IF CP=ConstraintPoint(xc) IF delta_k>0 THEN CP.right=QP ELSE CP.left=QP END IF tmp_constraints[num_active_constraints]=CP num_active_constraints++ ELSE #This constraint is quadratic. #Here x1 and x2 are sorted if they are real numbers x1,x2=SOLVE_QUADRATIC_EQUATION(b_kk,delta_k,A_k) IF (x1<=x2) AND (x1, x2 are real numbers) THEN IF b_kk>0 THEN IF update_integer THEN x1=CEIL(x1) x2=FLOOR(x2) END IF IF L0 penalty is used THEN QP=QUADRATIC_PENALTY(0,0,lagrangian) ELSE IF L1 penalty is used THEN QP=QUADRATIC_PENALTY(lagrangian * b_kk,lagrangian * delta_k,lagrangian * A_k) END IF #In this case, L2 penalty cannot be used CP=ConstraintPoint(x1) CP.left=QP tmp_constraints[num_active_constraints]=CP num_active_constraints++ CP=ConstraintPoint(x2) CP.right=QP tmp_constraints[num_active_constraints]=CP num_active_constraints++ ELSE IF b_kk<0 THEN IF update_integer THEN x1=FLOOR(x1) x2=CEIL(x2) END IF IF L0 penalty is used THEN QP=QUADRATIC_PENALTY(0,0,lagrangian) ELSE IF L1 penalty is used THEN QP=QUADRATIC_PENALTY(lagrangian * b_kk,lagrangian * delta_k,lagrangian * A_k) END IF CP=ConstraintPoint(x1) CP.right=QP tmp_constraints[num_active_constraints]=CP num_active_constraints++ CP=ConstraintPoint(x2) CP.left=-QP tmp_constraints[num_active_constraints]=CP num_active_constraints++ END IF ELSE IF x1, x2 are complex numbers THEN IF b_kk>0 THEN #Constraints are always violated IF L0 penalty is used THEN QP=QUADRATIC_PENALTY(0,0,lagrangian) ELSE IF L1 penalty is used THEN QP=QUADRATIC_PENALTY(lagrangian * b_kk,lagrangian * delta_k,lagrangian * A_k) END IF CP=ConstraintPoint(lower_bound) CP.right=QP tmp_constraints[num_active_constraints]=CP num_active_constraints++ END IF END IF END IF RETURN tmp_constraints,num_active_constraints

[0287] The following subroutine scans tmp_constraints and if there are two or more constraint points within the same value, it ensures that only one remains, which is their union.

[0288] Subroutine 3: COLLAPSE_CONSTRAINTS input: tmp_constraints num_active_constraints constraints Set num_collapsed_constraints=0 current_constraint=tmp_constraints[0] FOR c_i=1,num_active_constraints-1 DO next_constraint=tmp_constraints[c_i] IF current_constraint.value==next_constraint.value THEN current_constraint=current_constraint+next_constraint ELSE constraints[num_collapsed_constraints]=current_constraint current_constraint=next_constraint num_collapsed_constraint++ END IF END FOR constraints[num_collapsed_constraints]=current_constr num_collapsed_constraint++ RETURN constraints,num_collapsed_constraints

[0289] The following subroutine creates the subsegments, ie, the points after which the conditional probability distribution function changes coefficients.

[0290] Subroutine 4: MAKE_PARTIAL_SEGMENTS input: quadratic_bias linear_bias lower_bound upper_bound collapsed_constraints num_collapsed_constraints partial_segments update_integer / / For L1 penalty_sum=QUADRATIC_PENALTY(quadratic_bias,linear_bias,0); FOR i = 0, num_collapsed_constraints - 1 DO IF collapsed_constraints[i].value < lower_bound THEN penalty_sum += collapsed_constraints[i].right ELSE penalty_sum += collapsed_constraints[i].left END IF END FOR current_value = lower_bound num_partial_segments = 0 FOR i = 0, num_collapsed_constraints - 1 DO IF collapsed_constraints[i].value < lower_bound THEN CONTINUE END IF IF collapsed_constraints[i].value > lower_bound + update_integer THEN BREAK END IF IF collapsed_constraints[i].right exists THEN IF current_value!= collapsed_constraints[i].value THEN partial_segments[num_partial_segments].lower_bound = current_value partial_segments[num_partial_segments].quadratic_penalty = penalty_sum; num_partial_segments++; END IF It should be noted that the text seems to have some Japanese or other non - standard notations like "collapsed_constraints[i].rightが存在する" which is not a proper English expression. I've tried to translate it as literally as possible while making sense of the code - like structure. If this is from a specific technical or programming context, it might need more domain - specific knowledge for a more accurate translation.current_value=collapsed_constraints[i].value penalty_sum+=collapsed_constraints[i].right END IF IF collapsed_constraints[i].left exists THEN partial_segments[num_partial_segments].lower_bound=current_value partial_segments[num_partial_segments].quadratic_penalty=penalty_sum; current_value=constraints[i].value+update_integer; penalty_sum-=constraints[i].left; num_partial_segments++; END IF END FOR partial_segments[num_partial_segments].lower_bound=current_value partial_segments[num_partial_segments].quadratic_penalty=penalty_sum num_partial_segments++ RETURN partial_segments,num_partial_segments

[0291] The following subroutine creates a segment from the array of partial segments created by the subroutines above.

[0292] Subroutine 5: MAKE_SEGMENTS input: partial_segments num_partial_segments lower_bound upper_bound segments update_integer num_segments = num_partial_segments - 1; segment_idx = 0; IF num_segments > 0 THEN FOR i = 0, num_partial_segments - 2 DO IF partial_segments[i].lower_bound >= lower_bound THEN segments[segment_idx].lower_bound = partial_segments[i].lower_bound segments[segment_idx].upper_bound = partial_segments[i + 1].lower_bound - update_integer segments[segment_idx].quadratic_penalty = partial_segments[i].quadratic_penalty segment_idx++ ELSE num_segments -= 1 END IF END FOR ELSE # If there are no partial segments num_segments++; segments[0].lower_bound = partial_segments[0].lower_bound; segments[0].upper_bound = partial_segments[0].lower_bound; segments[0].quadratic_penalty = partial_segments[0].quadratic_penalty; END IF RETURN segments,num_segments

[0293] The final subroutine samples values ​​from the segment and differs based on whether the variable is integer (6a) or continuous (6b).

[0294] Subroutine 6a: SAMPLE_VALUE_FROM_SEGMENTS input: beta current_value segments num_segments update_integer FOR each segment DO IF update_integer THEN Calculate the segment distribution function using equation (8) ELSE Calculate the segment partition function using equations (9,10) END IF END FOR SAMPLE One segment is sampled with a probability proportional to its distribution function. SAMPLE Returns a new value using the inverse transform sampling method for this segment. RETURN new value

[0295] Subroutine 6b: SAMPLE_VALUE_FROM_SEGMENTS input: beta current_value segments num_segments update_integer Calculate log f0 from the segment containing current_value SAMPLE u=RAND(0,1) FOR each segment DO Calculate the subdomain where P(x)>=f0 Set the probability of each subregion proportional to its size END FOR SAMPLE Randomly select one subregion according to its respective probability SAMPLE Uniformly samples new values ​​from selected subregions. RETURN new value

[0296] In an alternative implementation, sampling of continuous variables can be done from a linear programming model with the values ​​of other variable types (i.e., discrete, integer, and binary) fixed. This can beneficially enable solving CQMs with continuous variables without adding multiple variables to the problem and can result in the continuous variables converging toward an optimal solution at a rate similar to the convergence of other problem variables. When solving optimization problems in which the variable types are continuous, sampling can be done based on a linear programming model. Linear programming, also known as linear optimization, is a mathematical model of variables with linear relationships. Methods for solving linear programming problems are well known in the art.

[0297] After each sweep of an optimization algorithm, such as simulated annealing, a linear programming problem can be created by fixing the values ​​of integer, binary, and discrete variables to their updated sample values. A solution to the linear programming problem can then be found to define updated values ​​for the continuous variables throughout the problem. The general structure of the linear problem does not change at each stage of the optimization algorithm, so there can be beneficially low overhead for solving a linear programming model at each stage. This can beneficially reduce the computational requirements for solving for continuous variables, since constructing a new linear programming model at each stage can incur large overhead, especially if the optimization algorithm requires many iterations to converge.

[0298] In one implementation, optimization can proceed for any integer, binary, or discrete variable in the problem, according to various implementations discussed above. Variables with quadratic interactions or interactions between continuous variables and other variable types can also have their updates sampled, such as by slice sampling as described above. Once the updates for these variables have been sampled, the remaining continuous variables can have their values ​​sampled from the linear programming problem, with the values ​​of the other variables held fixed at their updates.

[0299] It will be appreciated that linear programming models, as described below, do not support quadratic interactions between continuous variables or interactions between continuous variables and other variable types. It will be further appreciated that in some implementations, a combination of sampling from the conditional probability distributions discussed above and sampling from a linear programming model may be used. For example, continuous variables that have quadratic terms or interactions with other variable types may be sampled as discussed above using sampling from a conditional probability distribution with slice sampling. These continuous variables may then be fixed along with the other variable types, and the remaining continuous variables may have values ​​sampled from a linear programming model, as discussed in more detail below.

[0300] 14 is a flow diagram of an example of a method of operating a computing system to find a solution to a constrained or other problem that may be expressed as a constrained quadratic model, or to generate a sample solution that may be used as input to other algorithms, where the problem has one or more continuous variables. FIG. 14 shows method 1400. Method 1400 may be performed on a hybrid computing system including at least one digital or classical processor and a quantum processor, such as hybrid computing system 100 of FIG. 1, or may be performed by a classical computing system including at least one digital or classical processor.

[0301] Although method 1400 includes acts 1402-1444, one skilled in the art will understand that the number of acts shown is exemplary and that in some implementations certain acts may be omitted, additional acts may be added, and / or the order of the acts may be changed.

[0302] The method 1400 begins at 1402, for example, in response to a call or invocation from another routine or in response to input by a user.

[0303] At 1402, the processor receives a problem definition. The problem definition includes a set of variables, an objective function defined over the set of variables, and one or more constraint functions, each of the constraint functions defined by at least one variable in the set of variables. The problem definition may be received from a user (e.g., entered by an input device), transmitted from another processor, retrieved from memory, or provided as the output of another process executed by the processor. The objective function may be a quadratic function, and the problem definition may define a quadratic optimization problem. The set of variables may include one or more of continuous variables, discrete variables, binary variables, and integer variables. The constraint functions may include one or more quadratic equality or inequality constraint functions.

[0304] In one implementation, the problem definition is

number

number

number

number

number

number

number

[0305] The above problem definition does not include quadratic interactions between continuous variables and other variables. This model can be beneficially used to model mixed integer linear programming (MILP) problems as well as some cases of mixed integer quadratic programming (MIQP) problems.

[0306] At 1404, the processor initializes a sample solution to the objective function. The sample solution may be a random solution to the objective function. The random solution may be selected randomly from the entire variable space or may be selected within a range of the variable space based on known characteristics of the problem definition or other information. The sample solution may be generated by another algorithm or provided as input by a user.

[0307] At 1406, the processor initializes a progression parameter. The progression parameter can be a set of incrementally changing values ​​that define the optimization algorithm. For example, the progression parameter can be an inverse temperature, which can increment from an initial high temperature to a final low temperature. In some implementations, the inverse temperature is provided as a progression parameter for a simulated annealing algorithm. Selection of the inverse temperature can be performed as described above.

[0308] At 1408, the processor optionally initializes penalty parameters. As discussed in more detail below, in some implementations, the penalty parameters may be Lagrange multipliers. In other implementations, the penalty parameters may be selected as constant values ​​or values ​​that depend on one or more other variables.

[0309] At 1410, the processor initializes a linear programming problem of continuous variables using the initial states of the variables defined in act 1404.

[0310]

number

number

[0311] To give a linear programming model, we consider the part of the problem defined by the discrete, integer, and continuous variables as a function of constant values ​​(C0,C m ) can be set to:

number

number

[0312] As mentioned above, a penalty parameter can also be included in the model, given by:

number

number

[0313] At 1412, the processor optionally calculates the current value of each constraint function at the sample value. For example, the constraint function may be calculated as ΣA c,i x i +b c =0 or as an inequality.

[0314] At 1414, the processor increments a stage of the optimization algorithm, such as by providing a progress parameter to the optimization algorithm. The optimization algorithm may include simulated annealing, parallel tempering, Markov chain Monte Carlo techniques, branch-and-bound algorithms, and greedy algorithms, which may be executed by a classical computer. The optimization algorithm may also include algorithms executed by a quantum computer, such as quantum annealing, quantum-approximate optimization algorithms (QAOA) or other noisy intermediate-scale quantum (NISQ) algorithms, quantum-implemented fault-tolerant optimization methods, or other quantum optimization algorithms. The quantum computer may include a quantum annealing processor or a gate-model-based processor. Over successive iterations, the incremented optimization algorithm may provide samples.

[0315] At iteration k of an optimization algorithm (e.g., simulated annealing), sampling is performed for each integer, binary, and discrete variable as discussed below and in more detail above, and the problem is considered with the values ​​of the continuous variables fixed from the previous iteration (k-1):

number

number

[0316] At 1416, the processor selects an ith variable from a first set of variable sets. The first set of variables can include any integer, binary, or discrete variables, and the second set of variables can include any continuous variables. The ith variable has a current value from the sample initialized in act 1404.

[0317] At 1418, the processor determines the variable type of the variable. The variable type may be, for example, one of binary, discrete, and integer. It will be appreciated that a problem may include a mixture of these variable types, only one of these variable types, and combinations thereof.

[0318] At 1420, the processor selects a sampling distribution based on the variable type. In some implementations where the variable type is binary, the selected sampling distribution may be a Bernoulli distribution. In other implementations where the variable type is discrete, the selected sampling distribution may be a softmax distribution. Alternatively, in some implementations, discrete variables may be included as binary variables using a one-hot constraint, and the sampling distribution may be a Bernoulli distribution. In other implementations where the variable type is integer, sampling may be performed by Gibbs sampling, i.e., sampling based on conditional probability given the current state. Alternatively, Gibbs sampling can be performed for all variable types, with the conditional probability function determined by the variable type. Integer variables may also be converted to binary variables, and sampling from a Bernoulli distribution may be performed. Sampling for integer and continuous variables may also be performed by slice sampling, as described above.

[0319] At 1422, the processor determines a target energy bias to act on the variable under consideration given the current values ​​of all other variables. As noted above, the optimization problem can be structured with an objective function that defines energy, and during optimization, the processor returns a solution with the objective of reducing this energy. In some implementations, such as when the variable under consideration is a binary variable, the target energy change

number

[0320] At 1424, the processor determines a constraint energy bias acting on the variable under consideration given the current values ​​of all other variables and each of the constraint functions that include that variable. As discussed above, the constraint energy bias is Δ E =Σ c (|δ c,1 | n -|δ c,0 | n ) The penalty applied to each constraint may be determined by the magnitude of the violation. In some implementations, the processor determines a total energy change based on the objective energy change and the constraint energy change. The processor indirectly includes the constraint function in this energy by adding an energy term that takes the constraint into account.

[0321] At 1426, the processor samples updated values ​​of the variables from the sampling distributions based on the objective energy bias and constraint energy bias and the progress parameters as discussed in more detail above.

[0322] At 1428, the processor evaluates whether all of the variables in the first variable set have been considered. If all of the variables have not been considered, control returns to act 1416, where the next variable is selected. Once all of the variables have been considered, control passes to act 1430. In some implementations, the processor may incrementally consider each variable in the first variable set in turn, and evaluate all of the variables in the first variable set as having been considered when the last variable in the first variable set has been considered.

[0323] At 1430, the processor fixes the value of the first set of variables based on its updated value from above.

[0324] The processor solves a linear programming problem for the continuous variables of the second variable set and samples updated values ​​for the continuous variables at 1432. The linear programming problem can be solved using any technique known in the art, such as using an interior point method or a simplex algorithm.

[0325] When considering continuous variables, sampling is performed with the integer, binary, and discrete variables at fixed values, i.e., a given integer / binary state at iteration k, i.e., the variable

number

number

[0326] In the above equation, q m is a positive continuous variable, and u m is the Lagrange multiplier. q m is added to the left hand side of each constraint to avoid infeasibility in the linear programming problem. Then q m is penalized and added to the objective function at the current value of the Lagrange multiplier. The linear programming solver calculates q along with the continuous variables of the problem. m Define the value of independently. This mitigation technique avoids the infeasibility, but comes at a penalty.

[0327] At each stage of the optimization algorithm, the structure of the linear programming problem remains unchanged, but the right-hand side of the constraints and the variables q m The coefficients of are updated.

[0328] For large problems, a linear programming problem can be usefully divided into multiple smaller linear programming problems (subproblems), with continuous variables randomly assigned to one of the subproblems. This assignment can also be optimized to reduce the number of interactions between variables in different subproblems. When sampling for a variable in a given subproblem, the values ​​of continuous variables in all other subproblems are fixed to their most recent values. This can usefully allow the continuous variable to converge near or to a global optimum at the same rate as the other variables.

[0329] For example, given G as a set of continuous variables, n ={xi |i in G n The variables can be randomly divided into N groups of N, where N is the number of variables in a given set of variables. Variables and constraints can be defined for each group, and a linear programming problem can be constructed for each group of variables and constraints. Each linear programming problem can be solved sequentially, updating the variables involved.

[0330] At 1434, the processor updates the states of the continuous variables based on the outcome of the linear programming problem.

[0331] At 1436, the processor updates the energy biases of the objective and constraint functions based on the sampled values ​​of the continuous variables.

[0332] At 1438, the processor increments a progress parameter, for example, the temperature of the simulated annealing.

[0333] At 1440, a penalty parameter for each constraint may optionally be adjusted. In some implementations, the penalty parameter may be adjusted based on a change in energy of the constraint function defined by the i-th variable and the progress parameter. In other implementations, the penalty parameter may be adjusted as described in methods 800 and 900, discussed in more detail below.

[0334] At 1442, the processor evaluates one or more termination criteria. In some implementations, the termination criteria may be the value of a progress parameter. In other implementations, the termination criteria may include a number of iterations, an amount of time, an average energy change threshold between updates, a measure of the quality of the current value of a variable, or other metrics known in the art. If the termination criteria are not met, the method continues with act 1414. As described above, incrementing a stage of the optimization algorithm with respect to the objective function may include incrementing a simulated annealing algorithm or a parallel tempering algorithm. In other implementations, the optimization algorithm may be an MCMC algorithm or a greedy algorithm.

[0335] If one or more termination criteria are met, control passes to 1444 where the solution is output. At 1444, the method 1400 ends, e.g., until called again. The solution output in act 1444 can be passed to other algorithms, such as a quantum annealing algorithm.

[0336] In some implementations, after outputting the solution in act 1444, method 1400 can begin again with new samples that were initialized in act 1404. In some implementations, method 1400 can be run multiple times in parallel starting from different initialized samples or randomly generated sets of samples. In some implementations, the solutions can be paired and a binary problem can be constructed to evaluate the set of solutions using cross-Boltzmann updates.

[0337] An example implementation is described in the following pseudocode, which can be combined with the pseudocode above.

[0338] Pseudocode: For iteration k=0 β,L k ,x k ,y k Initialize x is initialized randomly and y is initialized to zero

number

number

number

number

number

[0339] The methods, processes, or techniques described above can be implemented by a series of processor-readable instructions stored on one or more non-transitory processor-readable media. Some example methods, processes, or techniques described above are performed in part by a dedicated device, such as an adiabatic quantum computer, or a quantum annealer, or a system for programming or controlling the operation of an adiabatic quantum computer or quantum annealer, e.g., a computer including at least one digital processor. While the methods, processes, or techniques described above can include various actions, those skilled in the art will understand that certain actions can be omitted and / or additional actions can be added in alternative embodiments. Those skilled in the art will understand that the order of actions shown is for illustrative purposes only and may be different in alternative embodiments. Some example actions or operations of the methods, processes, or techniques described above are performed iteratively. Some actions of the methods, processes, or techniques described above can be performed during each iteration, after multiple iterations, or at the end of all iterations.

[0340] The above description of exemplary implementations, including those described in the Abstract, is not intended to be exhaustive or to limit implementations to the precise form disclosed. While particular implementations and examples have been described herein for illustrative purposes, various equivalent modifications can be made by those skilled in the art without departing from the spirit and scope of the present disclosure. The teachings of the various implementations presented herein may be applied to other methods of quantum computing, not necessarily limited to the exemplary methods for quantum computing generally described above.

[0341] The various implementations described above can be combined to produce further implementations. All commonly owned U.S. published patent applications, U.S. patent applications, foreign patents, and foreign patent applications referenced herein and / or listed in Application Data Sheets are hereby incorporated by reference in their entirety, including, but not limited to, the following: U.S. Provisional Patent Application No. 62 / 951,749; U.S. Provisional Patent Application No. 63 / 174,097; U.S. Provisional Patent Application No. 63 / 250,466; U.S. Patent Application Publication Nos. 2014 / 0344322 and 2020 / 0234172, and U.S. Patent Nos. 7,533,068, 8,008,942, 8,195,596, 8,190,548, and 8,421,053.

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

Claims

1. 1. A method of operating a computing system for updating samples in an optimization algorithm to improve convergence to feasibility, said method being executed by a processor and comprising: receiving a problem definition including a set of variables, an objective function defined over the set of variables, and one or more constraint functions, each of the constraint functions defined by at least one variable of the set of variables; receiving sample values ​​of the set of variables and values ​​of a progress parameter; For each variable in the variable set, determining a variable type of said variable; selecting a sampling distribution based on said variable type; determining an objective energy bias based on the sample values ​​of the variables and one or more terms of the objective function that include the variables; determining one or more constraint energy biases based on the sample values ​​of the variables and each of the constraint functions defined by the variables; and sampling updated values ​​of the variables from the sampling distribution based on the objective energy bias, the one or more constraint energy biases, and the progress parameter; and returning an updated sample, the updated sample including the updated value for each variable in the variable set. A method comprising:

2. 2. The method of claim 1 , further comprising receiving a value of a penalty parameter, and wherein sampling the update value of the variable from the sampling distribution further comprises sampling the update value of the variable from the sampling distribution based on the value of the penalty parameter.

3. The method of claim 2 , wherein receiving a value for the penalty parameter comprises receiving a value for a Lagrangian parameter that depends on the value of the progress parameter.

4. The method of claim 1 , wherein determining the variable type of the variable comprises determining that the variable type is one of binary, discrete, integer, or continuous.

5. 5. The method of claim 4, wherein determining the variable type is one of binary, discrete, integer, or continuous comprises determining the variable type is binary, and selecting a sampling distribution based on the variable type comprises selecting a Bernoulli distribution.

6. 5. The method of claim 4, wherein determining the variable type is one of binary, discrete, integer, or continuous comprises determining the variable type is discrete, and selecting a sampling distribution based on the variable type comprises selecting a softmax distribution.

7. 5. The method of claim 4, wherein determining the variable type is one of binary, discrete, integer, or continuous comprises determining the variable type is one of integer or continuous, and selecting a sampling distribution based on the variable type comprises selecting a conditional probability distribution.

8. The method of claim 7 , wherein sampling the update values ​​of the variables from the sampling distribution comprises slice sampling from the conditional probability distribution.

9. The method of claim 1 , wherein receiving the value of the progression parameter comprises receiving an inverse temperature.

10. The method of claim 1, wherein sampling a sample set of values ​​for the set of variables from an optimization algorithm includes sampling a sample set of values ​​for the set of variables from one of a simulated annealing, parallel tempering, or quantum annealing algorithm.

11. The method of claim 1, wherein receiving a constrained quadratic optimization problem comprises a set of variables, an objective function defined over the set of variables, and one or more constraint functions, further comprising receiving a constrained quadratic optimization problem comprises a set of variables, a quadratic objective function defined over the set of variables, and one or more quadratic equality or inequality constraint functions.

12. A method of operating a hybrid computing system is provided, the hybrid computing system including a quantum processor and a classical processor, the method being performed by the classical processor; receiving a constrained quadratic optimization problem comprising a set of variables, an objective function defined over the set of variables, one or more constraint functions, each of the constraint functions defined by at least one variable of the set of variables, and progression parameters for the optimization, the progression parameters comprising a set of values ​​that increment between an initial value and a final value; iteratively until the final value of the progression parameter is reached; sampling a sample set of values ​​for said set of variables from an optimization algorithm; updating the sample set of values ​​using an update algorithm; For each variable in the variable set, determining a variable type of said variable; selecting a sampling distribution based on said variable type; determining an objective energy bias based on sample values ​​of the variable from the sample set of values ​​and one or more terms of the objective function that include the variable; determining one or more constraint energy biases based on the sample values ​​of the variables and each of the constraint functions defined by the variables; and sampling updated values ​​of the variables from the sampling distribution based on the objective energy bias, the one or more constraint energy biases, and the progress parameter; updating a sample set of said values, including: returning an updated sample, the updated sample including the updated value for each variable in the variable set; incrementing the progress parameter; transmitting the one or more final samples to a quantum processor; instructing the quantum processor to refine the sample; and outputting a solution including the refined sample. A method comprising:

13. 13. The method of claim 12, wherein transmitting one or more final samples to a quantum processor comprises transmitting pairs of samples to the quantum processor, and wherein instructing the quantum processor to refine the samples comprises instructing the quantum processor to perform quantum annealing to select between the samples.

14. The method of claim 12 , further comprising returning the output solution as a sample set of values ​​for the set of variables as input to the optimization algorithm.

15. The method of claim 12, further comprising receiving a value of a penalty parameter, and wherein sampling an updated value of the variable from the sampling distribution further comprises sampling an updated value of the variable from the sampling distribution based on the value of the penalty parameter.

16. The method of claim 15, wherein receiving the value of the penalty parameter comprises receiving a value of a Lagrangian parameter that depends on the value of the progression parameter.

17. The method of claim 12, wherein receiving the value of the progression parameter includes receiving an inverse temperature.

18. The method of claim 12, wherein determining the variable type of the variable includes determining that the variable type is one of binary, discrete, integer, or continuous.

19. A hybrid computing system, comprising: a quantum processor and a classical processor; at least one non-transitory processor-readable medium storing at least one of processor-executable instructions and data, the non-transitory processor-readable medium communicatively coupled to the classical processor, the non-transitory processor-readable medium configured to perform the method of any one of claims 12 to 18 in response to executing the at least one of the processor-executable instructions and data; A hybrid computing system comprising:

20. The hybrid computing system of claim 19, wherein the quantum processor includes a plurality of quantum bits communicatively coupled by a plurality of couplers, and the classical processor instructs the quantum processor to refine the sample using quantized annealing.

Images