Patent Publication Number: US-11645442-B1

Title: Method for reducing cost and increasing accuracy of variational quantum circuit optimization

Description:
CROSS-REFERENCE TO RELATED APPLICATION 
     This application is a continuation of U.S. patent application Ser. No. 16/840,935, filed on Apr. 6, 2020, which claims priority under 35 U.S.C. § 119(e) to U.S. Provisional Patent Application Ser. No. 62/829,745, “Method for Reducing the Cost and Increasing the Accuracy of Variation Quantum Circuit Optimization by Fixed-Point Jacobi Sweeps with Analytical Cluster Tomography and Augmented by Anderson Acceleration,” filed on Apr. 5, 2019, each of which is incorporated by reference in its entirety. 
    
    
     BACKGROUND 
     1. Technical Field 
     This disclosure relates to variational quantum circuit optimization, and in particular it relates to a Jacobi diagonalization and Anderson acceleration method for variational quantum algorithm parameter optimization. 
     2. Description of Related Art 
     Variational quantum algorithms are a general family of methods that offer a tractable route to obtaining quantum advantage over classical computing by utilizing the entanglements of a quantum circuit to produce a high-quality approximation to the ground state wavefunction of a Hamiltonian operator which encodes a physically or industrially interesting optimization problem. A challenge with variational quantum algorithms is the hybrid quantum-classical optimization of the parameters of the entangling quantum circuit to minimize the expectation value of the Hamiltonian observable. For variational quantum algorithms, it is generally found that converging the observable expectation value to the global minimum increases (e.g., maximizes) the quality of the solution. The optimizations are, however, remarkably difficult, with many local minima, large indefinite regions of the potential surface, and vast regions of “barren plateaus” with small gradients precluding straightforward use of conventional optimization techniques. 
     SUMMARY 
     Some embodiments relate to a hybrid method for improving an observable expectation value of a quantum circuit. The quantum circuit includes N nonoverlapping groups of G one-qubit rotation gates (G may be different for each group), where the G gates within each group are characterized by a parameter θ n  and the observable expectation value is a function of the parameters θ n  for all N groups. One step improves a subset M of the N groups. The selected subset of gates is characterized by parameters θ m . A non-quantum processor determines at least (2G+1) M  grid points sampling a domain of the parameters θ m . A quantum processor configures the parameters θ m  of the quantum circuit for each grid point and samples the observable expectation value of the quantum circuit configured at each grid point while holding the other parameters constant. A non-quantum processor determines a tomography function for the observable expectation value as a function of the parameters θ m . The tomography function is determined based on the sampled observable expectation values. A non-quantum processor determines values of the parameters θ m  that optimize the tomography function. 
     Some embodiments repeat this step to sweep through all of the parameters θ n . Each sweep includes a fixed number of iterations that sweep through all of the parameters θ n . Each iteration follows the step described above. In some embodiments, each sweep iterates through the same subsets in the same order. 
     Some embodiments may further improve the sweeps as follows. A non-quantum processor predicts a set of parameters θ n  that has an improved observable expectation value compared to the parameters θ n  resulting from the immediately previous sweep. The prediction is based on the parameters θ n  from previous sweeps and error values associated with the parameters θ n  from previous sweeps. 
     Other aspects include components, devices, systems, improvements, methods, processes, applications, computer readable mediums, and other technologies related to any of the above. 
    
    
     
       BRIEF DESCRIPTION OF THE DRAWINGS 
       Embodiments of the disclosure have other advantages and features which will be more readily apparent from the following detailed description and the appended claims, when taken in conjunction with the examples in the accompanying drawings, in which: 
         FIG.  1    is a flowchart of a method for optimizing a quantum circuit, according to an embodiment. 
         FIG.  2    illustrates an example quantum circuit diagram, a contour plot of a tomography function associated with the quantum circuit, and a flowchart of a method for obtaining a tomography function for the quantum circuit, according to an embodiment. 
         FIG.  3    is a detailed view of the contour plot in  FIG.  2   , according to an embodiment. 
         FIG.  4    illustrates a Jacobi sweep fixed-point iteration for the quantum circuit, according to an embodiment. 
         FIGS.  5 A- 5 B  illustrated additional Jacobi sweeps for the quantum circuit, according to some embodiments. 
         FIG.  6    is a flowchart of a method for performing an Anderson/Pulay DIIS Sequence Acceleration, according to an embodiment. 
         FIGS.  7 A- 7 D  are generalized quantum circuit diagrams, according to some embodiments. 
     
    
    
     The figures depict various embodiments for purposes of illustration only. One skilled in the art will readily recognize from the following discussion that alternative embodiments of the structures and methods illustrated herein may be employed without departing from the principles described herein. 
     DETAILED DESCRIPTION 
     The figures and the following description relate to preferred embodiments by way of illustration only. It should be noted that from the following discussion, alternative embodiments of the structures and methods disclosed herein will be readily recognized as viable alternatives that may be employed without departing from the principles of what is claimed. 
     Introduction 
     This disclosure is in the technical field of circuit-model quantum computing. It describes a set of strategies for optimizing the parameters of variational quantum algorithms (also known as “hybrid quantum-classical algorithms”), including but not limited to the quantum approximate optimization algorithm (QAOA), the variational quantum eigensolver (VQE) algorithm, and quantum neural network (QNN) algorithms. The strategies aim to improve the performance of variational quantum algorithms, where performance may be measured in terms of time-to-solution, solution accuracy, or a mixture of both. 
     This disclosure describes a new hybrid quantum-classical optimization technique that merges optimization methods with domain knowledge of the nature of the functional form or “tomography” of the quantum objective function to produce a new family of optimization methods which provide higher performance for variational quantum algorithm optimization. 
     A problem in variation quantum algorithm deployment is the optimization of the observable expectation value   (a real scalar function) as a function of the quantum circuit parameters {θ}. The quantum nature of the quantum circuit used to produce   implies that it is intractable to classically (e.g., analytically) determine the explicit function  ({θ}) for all possible values of {θ}. However, it is tractable to sample   a finite number of times at different values of {θ}. From such measurements, it is possible to reconstruct partial information about the general functional form of  ({θ}) (a process called “tomography”), and to use such information to perform the desired optimization of the quantum circuit parameters. To this point, gradient-based techniques such as L-BFGS (Limited-memory Broyden-Fletcher-Goldfarb-Shanno algorithm) and gradient-free techniques such as Powell&#39;s method have been applied to this problem. 
       FIG.  1    is a flowchart of a method for optimizing a quantum circuit, according to an embodiment. In particular,  FIG.  1    illustrates a method for reducing (e.g., finding the global minimum of) an observable expectation value   of a quantum circuit by optimizing parameters of the circuit.  FIG.  1    includes three nested stages: Stage 1 (aka, Jacobi step or local Jacobi move), Stage 2 (aka, Jacobi sweep), and Stage 3 (aka, Anderson/Pulay acceleration), which are further described below. Stage 1 is described with reference to  FIGS.  2  and  3   . Stage 2 is described with reference to  FIGS.  4 - 5   . Stage 3 is described with reference to  FIG.  6   . 
     Stage 1 Description 
       FIG.  2    illustrates a quantum circuit diagram  201 , a contour plot  203  of the tomography function associated with the quantum circuit  201 , and a flowchart  205  for obtaining the tomography function, according to an embodiment.  FIG.  3    is a detailed view of the contour plot  203 , according to an embodiment. 
     The circuit  201  is a quantum circuit (or a portion of a quantum circuit) with gates characterized by circuit parameters {θ}. It may be desirable to determine the observable expectation value  ({θ}) for all possible values of {θ}. The circuit  201  includes two qubits, six one-qubit rotational gates R y  (each with a circuit parameter θ), and two Controlled Not (CNOT) gates. In the example of  FIG.  2   , the gates R y  are generalized y-axis rotation gates with the form R y (θ)=e iθY , where Y is the Pauli-Y gate. The gates R y  are used for simplicity. Other one-qubit gates may be used, such as generalized x-axis and z-axis rotation gates. The circuit  201  is an example quantum circuit that will be referenced for example purposes throughout the rest of the disclosure. Other quantum circuits (e.g., with different gates, qubits, depth, etc.) may be solved by the methods described herein. For example, a given quantum circuit can be re-written to use entirely R y  gates (e.g., via a gate transformation). In another example, the methods described herein can be adapted to use the parametrization, tomography formulae, and sampling grids unique to the native gate set of the given quantum circuit. 
     In Stage 1, the solution of the optimization problem is solved locally for only a handful or “cluster” of active quantum circuit parameters, while the rest of the parameters are held constant (referred to as “frozen”). For example, in  FIG.  1   , the top middle gate and the bottom left gate are gates with active quantum circuit parameters while the circuit parameters of the remaining gates are frozen during the Stage 1 process. In some embodiments, the analytical form of  ({θ}) within the active cluster parameters is definitionally restricted to a discrete Fourier series of definite length, with Fourier coefficients (also referred to as “tomography coefficients”). The tomography coefficients are problem- and quantum-circuit-specific. The determination of the tomography coefficients can be done exactly (e.g., in the absence of device or statistical noise channels) by evaluating  ({θ}) on a small rectilinear Fourier quadrature grid of the active cluster parameters, followed by a discrete Fourier transformation. This fully determines the analytical tomography function for  ({θ}) within the active cluster parameters. Since the maxima of some of the analytical tomography functions cannot be determined analytically (e.g., when a cluster includes more than one parameter), classical optimization techniques, such as Newton-Raphson, L-BFGS, or iterative analytical solutions in alternating directions, are then applied (e.g., offline from the quantum processor) to optimize the circuit parameters within the cluster based on the tomography function. 
     As previously stated, the contour plot  203  illustrates an example tomography function for circuit  201 . Since all but two of the circuit parameters θ are held constant, the plot only includes two circuit parameter domains, labeled θ A  and θ B . The tomography function includes two maxima  301  and two minima  303 . The optimization steps  305 A and  305 B illustrate iterations of a classical optimization algorithm to find minima  303  (e.g., a global minima) of the tomography function. Due to the sinusoidal nature of the tomography function, the two minima  301  may be the same value. 
     The contour plot  203  also includes a Fourier quadrature grid. Since, the tomography function is restricted to a discrete Fourier series of definite length, the tomography function can be determined by sampling the observable expectation value of the quantum circuit at the grid points  307 . Since circuit  201  has two active circuit parameters, plot  203  includes nine grid points  307  (described further below). The grid points  307  are evenly spaced apart at values of θ=0, ±π/3, although this not required. The evenly spaced choice may provide the best mode for (e.g., rapid) statistical convergence of the tomography parameters with respect to number of quantum measurements. In other embodiments, three different values of θ may be sufficient to determine the tomography function. 
     Example Method Steps of Stage 1 
     In some embodiments, Stage 1 is a quadrature-based method to experimentally obtain the analytical tomography formula for the observable expectation value of a quantum circuit across all possible quantum circuit parameters by evaluating the observable expectation value at a small discrete set of parameters (a “quadrature”). For example, see flowchart  205 . The steps of this method may include: 
     A subset of active circuit parameters (referred to as a “cluster”) is selected. The subset maybe a proper subset of the circuit parameters. If the total number of circuit parameters being optimized is small enough (e.g., less than five) or the computational power of the system is high enough, then all of the circuit parameters may be selected. In these situations, Stages 2 and 3 may not be necessary to optimize the circuit parameters. 
     A Fourier quadrature grid of appropriate size and dimension is defined for the considered subset of active circuit parameters. A general cluster of M sets of circuit parameters (the sets of circuit parameters may be overlapping sets), each determining the parameter for G one-qubit (e.g., rotation) gates may use a rectilinear grid with a total of (2G+1) M  Fourier grid points evenly spaced in the periodic domain of the circuit parameters (e.g., see grid points  307 ). Note that G represents the number of one-qubit gates in each of the M sets that are pinned together (e.g., have the same circuit parameter). 
     The analytical tomography function (also referred to as formula) is defined for the observable expectation value of the quantum circuit as a function of (any possible) combinations of parameters within the active cluster. In general, the tomography function is the rectilinear M-dimensional discrete Fourier series, with 2G+1 Fourier basis functions (complex roots of unity) in each dimension. The linear coefficients of the Fourier series are the quantum-circuit-specific “tomography coefficients.” 
     The observable expectation value   of the quantum circuit is sampled at each Fourier quadrature grid point. The frozen parameters may have predefined values (e.g., θ i =0 at the beginning of an optimization, or prior determined values during optimization). 
     A discrete Fourier transformation is performed on the set of observable expectation values   sampled at the Fourier quadrature grid points to obtain the tomography coefficients and complete the determination of the analytical tomography formula for the active cluster of circuit parameters. 
     Locally optimized circuit parameters for the cluster are analytically determined based on the determined tomography function. Optimized circuit parameters may result in the tomography function for the active cluster being at a minimum value within a threshold error. 
     Example Cases of Stage 1 
     A 3-point Fourier quadrature grid may be determined for an M=1, G=1 cluster using a single one-qubit gate (as often encountered in VQE). The analytical tomography formula for this single one-qubit gate is: O(θ)=α+β cos(2θ)+γ sin(2θ), where α, β, and γ are the tomography coefficients. By determining the tomography coefficients using the 3-point quadrature grid, the tomography function may be determined. 
     A 9-point Fourier quadrature grid may be determined for an M=2, G=1 cluster using a pair of one-qubit gates (as often encountered in VQE). 
     A 3 M -point Fourier quadrature grid may be determined for an M=M, G=1 cluster using M one-qubit gates (as often encountered in VQE). Thus, this grid resolves the analytical tomography formula for a cluster of M one-qubit gates. 
     A (2G+1) M -point Fourier quadrature may be determined for an M=M, G=G cluster using M groups of G one-qubit gates (as often encountered in QAOA). Thus, this grid resolves the analytical tomography formula for a cluster of M groups of G one-qubit gates. 
     Example Tomography Functions 
     1. One-Gate Tomography 
     Consider a general R y  gate quantum circuit of the form seen in  FIG.  7 A  (e.g., as encountered in VQE). The corresponding observable expectation value is,
 
 (θ A )≡ 0| Û   †   R   y   † (θ A ) {circumflex over (V)}   †   Ô{circumflex over (V)}{circumflex over (R)}   y (θ A ) Û| 0   (1)
 
Here Ô is a Hermitian operator, and Û and {circumflex over (V)} are unitary operators, all in the 2 N -dimensional Hilbert space for N qubits. These matrices are problem-specific, and their construction may be determined before the invocation of the quantum algorithm. {circumflex over (R)} y (θ A )≡e −iθ     A     Ŷ     A    is a one-qubit rotation gate acting on an arbitrary qubit corresponding to index A. The observable expectation value  (θ A ) is periodic for
 
     
       
         
           
             
               
                 
                   θ 
                   A 
                 
                 ∈ 
                 
                   [ 
                   
                     
                       - 
                       
                         π 
                         2 
                       
                     
                     , 
                     
                       
                         + 
                         π 
                       
                       / 
                       2 
                     
                   
                 
               
               ) 
             
             . 
           
         
       
     
     By inspection, this has the definitional tomography,
 
 (θ A )≡α+β cos(2θ A )+γ sin(2θ A )  (2)
 
Here, {α, β, γ} are the “tomography parameters,” whose values depend on the specifics of the problem (e.g., the operators Û, {circumflex over (V)}, and Ô).
 
     To determine the tomography parameters, we can define the three-point Fourier quadrature,
 
   − ≡ (−π/3),   0 ≡ (0),   + ≡ (+π/3)  (3)
 
This has the matrix form,
 
                           [         1           -   1     /   2             3   /   2               1       1       0           1         1   /   2             3   /   2             ]     ︸         T   ^       (   1   )         [         α           β           γ         ]     =     [           𝒪   -               𝒪   0               𝒪   +           ]             (   4   )               
The inverse is
 
     
       
         
           
             
               
                 
                   
                     
                       T 
                       ˆ 
                     
                     
                       ( 
                       1 
                       ) 
                     
                     
                       - 
                       1 
                     
                   
                   = 
                   
                     
                       1 
                       3 
                     
                     [ 
                     
                       
                         
                           1 
                         
                         
                           1 
                         
                         
                           1 
                         
                       
                       
                         
                           
                             - 
                             1 
                           
                         
                         
                           2 
                         
                         
                           
                             - 
                             1 
                           
                         
                       
                       
                         
                           
                             - 
                             
                               3 
                             
                           
                         
                         
                           0 
                         
                         
                           
                             3 
                           
                         
                       
                     
                     ] 
                   
                 
               
               
                 
                   ( 
                   5 
                   ) 
                 
               
             
           
         
       
     
     2. Two-Gate Tomography 
     Now consider a general two R y  gate quantum circuit of the form seen in  FIG.  7 B  (e.g., as encountered in VQE). The corresponding observable expectation value is,
 
 (θ A ,θ B )≡ 0| Û   †   {circumflex over (R)}   y   † (θ A ) {circumflex over (V)}   †   {circumflex over (R)}   y   † (θ B ) Ŵ   †   ×ÔŴ{circumflex over (R)}   y (θ B ) {circumflex over (V)}{circumflex over (R)} (θ A ) Û| 0   (6)
 
     By inspection, this has the tomography
 
 (θ A θ B )≡αα+αβ cos(2θ B )+αγ sin(2θ B )+βα cos(2θ A )+ββ cos(2θ A )cos(2θ B )+βγ cos(2θ A )sin(2θ B )+γα sin(2θ A )+γβ sin(2θ A )cos(2θ B )+γγ sin(2θ A )sin(2θ B )  (7)
 
     We can define a nine-point Fourier quadrature,
 
   −− ≡ (−π/3,−π/3),   −1 ≡ (−π/3,0),
 
   −+ ≡ (−π/3,+π/3),   0− ≡ (0.−π/3),
 
   00 ≡ (0,0),   0+ ≡ (0,+π/3)
 
   +− ≡ (+π r/ 3,−π/3),   +0 ≡ (+π/3,0)
 
   ++ ≡ (+π/3,+π/3)  (8)
 
     This has the matrix form, 
     
       
         
           
             
               
                 
                   
                     
                       
                         T 
                         ˆ 
                       
                       
                         ( 
                         2 
                         ) 
                       
                     
                     [ 
                     
                       
                         
                           
                             α 
                             ⁢ 
                             α 
                           
                         
                       
                       
                         
                           
                             α 
                             ⁢ 
                             β 
                           
                         
                       
                       
                         
                           
                             α 
                             ⁢ 
                             γ 
                           
                         
                       
                       
                         
                           
                             β 
                             ⁢ 
                             α 
                           
                         
                       
                       
                         
                           
                             β 
                             ⁢ 
                             β 
                           
                         
                       
                       
                         
                           
                             β 
                             ⁢ 
                             γ 
                           
                         
                       
                       
                         
                           
                             γ 
                             ⁢ 
                             α 
                           
                         
                       
                       
                         
                           
                             γ 
                             ⁢ 
                             β 
                           
                         
                       
                       
                         
                           
                             γ 
                             ⁢ 
                             γ 
                           
                         
                       
                     
                     ] 
                   
                   = 
                   
                     [ 
                     
                       
                         
                           
                             𝒪 
                             -- 
                           
                         
                       
                       
                         
                           
                             𝒪 
                             
                               - 
                               0 
                             
                           
                         
                       
                       
                         
                           
                             𝒪 
                             
                               - 
                               + 
                             
                           
                         
                       
                       
                         
                           
                             𝒪 
                             
                               o 
                               - 
                             
                           
                         
                       
                       
                         
                           
                             𝒪 
                             
                               0 
                               ⁢ 
                               0 
                             
                           
                         
                       
                       
                         
                           
                             𝒪 
                             
                               0 
                               + 
                             
                           
                         
                       
                       
                         
                           
                             𝒪 
                             
                               + 
                               - 
                             
                           
                         
                       
                       
                         
                           
                             𝒪 
                             
                               + 
                               0 
                             
                           
                         
                       
                       
                         
                           
                             𝒪 
                             
                               + 
                               + 
                             
                           
                         
                       
                     
                     ] 
                   
                 
               
               
                 
                   ( 
                   9 
                   ) 
                 
               
             
           
         
       
     
     Here {circumflex over (T)} (2) ≡{circumflex over (T)} (1) ⊗{circumflex over (T)} (1)  and the corresponding inverse is, {circumflex over (T)} (2)   −1 ={circumflex over (T)} (1)   −1 ⊗{circumflex over (T)} (1)   −1 . 
     3. M-Gate Tomography 
     Now consider a general M Ry gate quantum circuit of the form seen in  FIG.  7 C  (e.g., as encountered in VQE). The corresponding observable expectation value is,
 
 (θ A ,θ B , . . . ,θ Z )≡ 0| Û   †   {circumflex over (R)}   y   † (θ A ) {circumflex over (V)}   †   {circumflex over (R)}   y   † (θ B ) Ŵ   †    . . . {circumflex over (R)}   y   † (θ Z ) {circumflex over (Z)}   †   ×Ô{circumflex over (Z)}{circumflex over (R)}   y (θ Z ) . . .  Ŵ{circumflex over (R)}   y (θ B ) {circumflex over (V)}{circumflex over (R)}   y (θ A ) Û| 0   (10)
 
     Note that the alphabetical labeling θ A , θ B , . . . , θ Z  is merely illustrative—there can be more or less than 26 gates depending on M. This has the 3 M -parameter tomography, 
                     𝒪   ⁡   (       θ   A     ,     θ   B     ,   …         ,     θ   Z       )     =       ∑     I   →︀           C     I   →︀       ⁢       ∏       i   D     ∈   I           ∅     i   D       (     θ   D     )                   (   11   )               
where {right arrow over (I)} ranges over the set of trinary strings of length M, e.g., 000, 001, 002, 010, 011, 012, 020, 021, 022, 100, etc., and the basis functions are Ø 0 (θ)≡1, Ø 1 (θ)≡cos(2θ), and Ø 2 (θ)≡sin(2θ).
 
     The tomography coefficients {c {right arrow over (I)} } can be computed from a 3 M -point quadrature grid consisting of a Cartesian grid of {−π,3,0+π/3} in each angle. The transfer matrix is {circumflex over (T)} (M) ≡⊗ M {circumflex over (T)} (1) , and the corresponding inverse is {circumflex over (T)} (M)   −1 =⊗ M {circumflex over (T)} (1)   −1 . 
     4. QAOA Tomography 
     In QAOA, we often encounter a generalization or simplification where the parameters of multiple R y  gates are pinned together, e.g., an M-stage quantum circuit where each stage D has G D  {circumflex over (R)} y (θ D ) gates as seen in  FIG.  7 D . The corresponding observable expectation value is,
 
 (θ A , . . . ,θ Z )≡ 0| Û   †   {circumflex over (R)}   y   † (θ A ) {circumflex over (V)}   †   {circumflex over (R)}   y   † (θ A ) Ŵ   †    . . . {circumflex over (R)}   y   † (θ Z ) Ŷ   †   {circumflex over (R)}   y   † (θ Z ) {circumflex over (Z)}   †   ×Ô{circumflex over (Z)}{circumflex over (R)}   y (θ Z ) Ŷ{circumflex over (R)}   y (θ Z ) . . .  Ŵ{circumflex over (R)}   y (θ A ) {circumflex over (V)}{circumflex over (R)}   y (θ A ) Û| 0   (12)
 
     This has the Π D (2G D +1)-parameter tomography [(2G+1) M  if G D  is independent of D], 
                     𝒪   ⁡   (       θ   A     ,   …         ,     θ   Z       )     =       ∑     I   →︀           C     I   →︀       ⁢       ∏       i   D     ∈     I   →︀             ∅     i   D       (     θ   D     )                   (   13   )               
where each digit of {right arrow over (I)} ranges from −G D  to +G D  (inclusive. The basis functions are sin(G D ·2θ), . . . , sin(2·2θ), sin(1·2θ) 1, cos(1·2θ) cos(2·2θ), . . . , cos(G D ·2θ).
 
     Two technical notes with the circuit above: (1) in QAOA, often the interstitial operators Û, {circumflex over (V)}, . . . are the identity, e.g., within commuting layers of 1-qubit driver terms and (2) above we have drawn the θ A  and θ Z  stages as disjoint, but they may be interleaved without changing the analysis. 
     For this case, a Cartesian product of (2G D +1)-point Fourier grids is sufficient to resolve the tomography. The task is to redefine the grid, transfer matrix, and transfer matrix inverse. For general G D , the Fourier grid is, 
     
       
         
           
             
               
                 
                   
                     
                       
                         θ 
                         p 
                       
                       = 
                       
                         
                           
                             
                               ( 
                               
                                 p 
                                 + 
                                 1 
                               
                               ) 
                             
                             ⁢ 
                             π 
                           
                           
                             
                               2 
                               ⁢ 
                               
                                 G 
                                 D 
                               
                             
                             + 
                             1 
                           
                         
                         - 
                         
                           π 
                           2 
                         
                       
                     
                     , 
                     
                       p 
                       ∈ 
                       
                         [ 
                         
                           0 
                           , 
                           
                             
                               2 
                               ⁢ 
                               
                                 G 
                                 D 
                               
                             
                             + 
                             1 
                           
                         
                       
                     
                   
                   ) 
                 
               
               
                 
                   ( 
                   14 
                   ) 
                 
               
             
           
         
       
     
     The transfer matrix is, 
     
       
         
           
             
               
                 
                   
                     
                       T 
                       
                         i 
                         ⁢ 
                         p 
                       
                       G 
                     
                     ≡ 
                     
                       
                         ∅ 
                         i 
                       
                       ( 
                       
                         θ 
                         p 
                       
                       ) 
                     
                   
                   = 
                   
                     
                       [ 
                       
                         
                           
                             sin 
                             ⁡ 
                             ( 
                             
                               2 
                               ⁢ 
                               i 
                               ⁢ 
                               
                                 θ 
                                 p 
                               
                             
                             ) 
                           
                           1 
                         
                         
                           cos 
                           ⁡ 
                           ( 
                           
                             2 
                             ⁢ 
                             i 
                             ⁢ 
                             
                               θ 
                               p 
                             
                           
                           ) 
                         
                       
                       ] 
                     
                     _ 
                   
                 
               
               
                 
                   ( 
                   15 
                   ) 
                 
               
             
           
         
       
     
     As before, the M-stage transfer matrix is {circumflex over (T)} (M) ≡⊗ D   M {circumflex over (T)} (1)   G     D   , and the corresponding inverse is {circumflex over (T)} (M)   −1 ≡⊗ D   M {circumflex over (T)} (1)   G     D     ,−1 . 
     5. Further Discussion of the Tomography Function 
     In closing this section, we make two observations regarding the Fourier quadrature resolution of the tomography parameters. The first observation is that the quadrature points may be widely spaced in θ, and the tomography coefficients may be resolved by the transfer matrix inverse with coefficients that are nearly unity. This means that roughly the same statistical convergence of the observable expectation value at each quadrature grid point may be required to obtained a given absolute accuracy in the observable across the full tomography formula as is required to obtain the same accuracy in the observable at a specific point. This is in marked contrast to, e.g., finite difference derivatives, where the observable must be resolved to much higher precision at the stencil grid points to obtain a given accuracy in the approximated derivative, due to subtractive cancelation. Since we are free to analytically differentiate the tomography formula, we obtain a recipe for the analytical gradient of the observable expectation which has the same number of required observable expectation value evaluations as the second-order symmetric finite difference formula, but with markedly reduced precision requirements. 
     Stage 2 Description 
     Stage 2 of the present disclosure relates to the sequential use of the cluster optimization technique described above (Stage 1) to develop a global optimization technique known as the “Jacobi sweep” fixed-point iterative procedure. For a given circuit, a Jacobi sweep may be defined as a specified sequence of active clusters. Each circuit parameter of the circuit may be specified in at least one of the clusters in the sequence. Thus, Stage 1 may be repeated for each of the clusters specified in the sweep (e.g., see Stage 2 in  FIG.  1   ).  FIGS.  4 - 5    illustrate example Jacobi sweeps for the quantum circuit  201 , according to some embodiments. Note, that  FIGS.  5 A and  5 B  do not illustrate all of the active clusters in the sequence. 
     Different types of Jacobi sweeps may be performed in Stage 2. For instance, the “Jacobi-1” method (e.g., see  FIG.  5 A ) sweeps over single parameter clusters. The “Jacobi-2” method (e.g., see  FIG.  5 A ) sweeps over pairs of parameters. Jacobi-2 may provide faster per-sweep convergence compared to Jacobi-1 but at the cost of more clusters. For example, for circuit  201 , a complete Jacobi-1 sweep includes six different iterations, while a complete Jacobi-2 sweep includes fifteen iterations (example iterations for both sweeps are illustrated in  FIG.  5 A ). Furthermore, since two parameters are optimized in each iteration of Jacobi-2, the iterations may be more computationally expensive than the single parameter iterations of the Jacobi-1 sweep. 
     The “Jacobi-A” and “Jacobi-B” methods (e.g., see  FIGS.  5 A and  5 B ) are modified versions of the Jacobi-2 sweep. The Jacobi-A and Jacobi-B sweeps lie between Jacobi-1 and Jacobi-2 (in terms of computational cost) by using pairs of parameters that are restricted to a single qubit wire (Jacobi-A) or to single or adjacent qubit wires (Jacobi-B). Depending on the circuit being optimized, the global convergence rate of the Jacobi-A or Jacobi-B method may be comparable to the global convergence rate of the Jacobi-2 method. 
     Higher order sweeps, such as Jacobi-3 (e.g., see  FIG.  5 B ) or higher, may also be used. However, as described when comparing Jacobi-2 with Jacobi-1, each iteration of a higher order sweep may be more computationally expensive. For example, as described above with reference to Stage 1, the number of grid points needed to determine the tomography function in a cluster iteration may increase according to 3 M  (e.g., if G=1). 
     Multiple Jacobi sweeps may be performed, for example until a global convergence is reached. Thus, Stage 2 may result in a global convergence progress for quantum circuits with dozens or more parameters while making only a series of local cluster moves, where each cluster may be optimization over only a handful of circuit parameters. If multiple Jacobi sweeps are performed, a sweep may be different or the order of the clusters in the sweep may be different compared to the other sweeps. For example, a first sweep may use a Jacobi-1 method and a second subsequent sweep may use a Jacobi-2B method. In another example, a Jacobi-1 sweep is performed multiple times, and the order of the clusters may be different for each sweep. However, if Stage 3 is performed, the type of sweep and the order of the clusters may be the same for each sweep. 
     Example Method Steps of Stage 2 
     In some embodiments, Stage 2 is a Jacobi-sweep fixed-point iterative method to perform global optimization of the quantum circuit by considering active clusters of circuit parameters and sweeping iteratively over the clusters. The steps of the method may include: 
     A sweep is specified. As previously stated, a sweep is a sequence of clusters of circuit parameters. The sweep covers some, all, or even repeated subsets of the full set of circuit parameters. 
     For each cluster in the sweep, the active circuit parameters may be optimized. To do this, the analytical tomography formula for a cluster of active circuit parameters may be determined (the other circuit parameters are frozen). For example, the Fourier quadrature-based method described in Stage 1 is used to determine an analytical classical tomography function for the cluster. Afterwards, as also described in Stage 1, the parameters in the cluster may be optimized based on the tomography function. Thus, the observable expectation value may be optimized for the active parameters. Note that Stage 2 is not required to perform the steps described in Stage 1. Other methods for optimizing active circuit parameters in a cluster may be used. The frozen parameters may have predefined initial values (e.g., θ i =0), optimized values determined from a previous sweep, and/or optimized values determined from a previous iteration in the current sweep (e.g., one or more circuit parameters are updated after each iteration in a sweep). 
     One or more sweeps are repeated until the observable expectation value has converged to desired tolerances. For example, changes in the observable expectation value between sweeps are less than a threshold value, changes in the optimized parameters between sweeps are bellow a threshold value, or the norm of the gradient of the tomography function with respect to the optimized parameters is below a threshold value. After a sweep is completed, initial circuit parameters for a new sweep may be the optimized circuit parameters determined from the previous sweep. Thus, the observable expectation value may converge toward a global solution (e.g., minimum value). 
     Example Cases of Stage 2 
     Specific realizations of this method as might be seen in VQE are: 
     Jacobi-1: all M=1 single one-qubit gates are selected to define the clusters in the sweep. 
     Jacobi-2: all M=2 pairs of one-qubit gates are selected to define the clusters in the sweep. 
     Jacobi-A: all M=2 pairs of one-qubit gates within the same qubit wire are selected to define the clusters in the sweep. This may yield a method with cost and accuracy between that of Jacobi-1 and Jacobi-2. 
     Jacobi-B: all M=2 pairs of one-qubit gates within the same and adjacent qubit wires are selected to define the clusters in the sweep. This may yield a method with cost and accuracy between that of Jacobi-1 and Jacobi-2. 
     Jacobi-Gen: generic clusters of user-specified one-qubit gates, including clusters of different sizes and clusters with M&gt;2 active circuit parameters. 
     Other variations include extensions of these specific realizations to G&gt;1 groups of one-qubit gates as might be seen in QAOA. Additional variations include simple variations of the selection and ordering of clusters. 
     Stage 3 Description 
     Stage 3 of the present disclosure relates to the addition of the Anderson acceleration methodology to accelerate the convergence of the Stage 2 Jacobi sweep optimization (e.g., see Stage 3 in  FIG.  1   ). Anderson acceleration examines the iterative history of circuit parameters and corresponding estimates of the error vectors. The relationships of the error between different history vectors are considered, and a new extrapolated proposal for the optimized circuit parameters is proposed as a linear combination of the circuit parameters in the iterative history (e.g., see  FIG.  6   ). The Anderson acceleration process may add insignificant cost to each Jacobi sweep iteration and may improve convergence behavior. 
     Two variants of Stage 3 include: a case where the total displacement of the Jacobi move (e.g., the difference between of the parameters of a previous iteration and parameters of the current iteration) is used as the error vector (referred to as “Anderson”) and a case where the gradient vector of the partial derivatives of the observable expectation value with respect to the circuit parameters is used as the error vector (referred to as “Pulay DIIS”). These are found to substantially improve convergence behavior, particularly when combined with the Jacobi-1 method. 
     Example Method of Stage 3 
     In some embodiments, Stage 3 is an “Anderson acceleration” extension to the Jacobi-sweep algorithm (Stage 2) that examines the iteration history of previous Jacobi sweeps and proposes an accelerated/extrapolated sequence of Jacobi sweeps. The proposed sequence of moves may be based on relationships between the “error vectors” of moves or other error proxies in the iterative history. The steps of the method may include: 
     At each completion of the Jacobi-sweep algorithm (Stage 2), store an entry in the iterative history by recording the current set of circuit parameters (the “state vector”) and some approximation of the errors of the current circuit parameters (the “error vector”). Note that Stage 2 is not required to perform the Anderson acceleration steps described in Stage 3. Other methods for iteratively optimizing circuit parameters may be used instead. 
     From the contents of the iterative history, construct an extrapolated estimate for a proposed set of circuit parameters that is composed of a linear combination of the iterative history of state vectors and that is predicted to have lower error (the proposed set of parameters is labeled {θ′ n } in  FIG.  6   ). 
     Example Cases of Stage 3 
     The error vector is determined based on the total displacement of the circuit parameters during each Jacobi sweep. This may be referred to as “Anderson.” 
     The error vector is determined based on the gradient vector of the partial derivatives of the observable expectation value with respect to the circuit parameters. This may be referred to as “Pulay DIIS.” 
     The specific constituents of the state and error vectors or the mechanism of extrapolation of the Anderson acceleration procedure are not critical, and many variations of these methods exist. 
     Additional Explanation of Stage 3 
     Stage 3 may substantially accelerate the global convergence of the Jacobi-type algorithm developed (Stage 2) above, at little or no additional cost (e.g., no additional quantum measurements, and very few classical operations). 
     The idea of Anderson is, for a given set of non-linear equations f({right arrow over (θ)})=0 (the “vector” symbol indicates the dimension of the parameter space), one is given an iterative sequence over iteration index i of state vectors {right arrow over (θ)} i  and error vectors {right arrow over (∈)} i ≡{right arrow over (θ)} i −{right arrow over (θ)} so that {right arrow over (θ)} i ≡{right arrow over (θ)} i +{right arrow over (∈)} i . Anderson acceleration posits that at iteration k, an improved state vector {right arrow over (θ)} ik  can replace {right arrow over (θ)} k , and is composed of a linear combination of the current history of iterative state vectors, 
                       θ   →       ′   ⁢   k       ≡       ∑   i         c   i   k     ⁢       θ   →     i                 (   16   )               
The corresponding improved error vector is,
 
     
       
         
           
             
               
                 
                   
                     
                       ∈ 
                       → 
                     
                     
                       ′ 
                       ⁢ 
                       k 
                     
                   
                   
                     ≡ 
                     
                       
                         ∑ 
                         i 
                       
                       
                         
                           c 
                           i 
                           k 
                         
                         ⁢ 
                         
                           
                             θ 
                             → 
                           
                           i 
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   17 
                   ) 
                 
               
             
           
         
       
     
     Anderson acceleration proposes that the coefficients c i   k  be chosen to minimize the square of the 2-norm of {right arrow over (∈)} k . 
                       0   k     ⁢     (     c   i     )       =         ∑     i   ⁢   j           c   i   k     ⁢     c   j   k             ∈   →     i       ·       ∈   →     j                 (   18   )               
subject to the normalization condition,
 
     
       
         
           
             
               
                 
                   
                     
                       ∑ 
                       i 
                     
                     
                       c 
                       i 
                       k 
                     
                   
                   = 
                   1 
                 
               
               
                 
                   ( 
                   19 
                   ) 
                 
               
             
           
         
       
     
     In the limit of a sufficiently large iterative space to span the vector space for {right arrow over (θ)}, this is guaranteed to produce the desired result of {right arrow over (θ)} ik ≡θ k , and will provide a least-squares error approximation in a less-complete limit. An analytical solution for the linearly-constrained least-squares problem is obtained, 
                       [           B   11         …         B     1   ⁢   k             -   1             ⋮       ⋱       ⋮       ⋮             B     k   ⁢   1           …         B   kk           -   1               -   1         …         -   1         0         ]     [           c   1   k             ⋮             c   k   k             λ         ]     =     [         0           ⋮           0             -   1           ]             (   20   )               
Here B ij ≡{right arrow over (∈)} i ·{right arrow over (∈)} j  and γ is a Langrange multiplier corresponding to the normalization constraint. The numerical cost of predicting θ′ k  from the iterative history is minimal (and all classical): only the inner-products inner-products in B ij , a matrix inversion of the dimension of k+1, and the vector addition to form θ′ k  are required.
 
     The cognizant reader will have noticed that the above manipulations would all be for naught if we actually had the iterative history of the error vector {right arrow over (∈)} k : at any iteration (including the first iteration), we could simply obtain the full solution as θ=θ k −{right arrow over (∈)} k . In practice, Anderson acceleration replaces the iterative history of the exact error vector {right arrow over (∈)} k  with a proxy quantity   that is pragmatically chosen to approximate a scaled error vector. For instance, the displacement of the parameters during a fixed-point iterative move δ{right arrow over (θ)} k ≡{right arrow over (θ)} k ={right arrow over (θ)}′ k-1  (the original Anderson and Pulay DIIS proposals) may be effective slowly converging fixed-point series which are making large numbers of small moves in roughly the same direction. Alternatively, an approximately preconditioned gradient may prove more useful for accelerating nonlinear optimization procedures posed as gradient root finding or linear solve procedures posed as residual zero finding. 
     Note that the invocation of Anderson or Pulay DIIS sequence acceleration techniques may invalidate the monotonic convergence property of the Jacobi diagonalization-type procedure. However, we note that Anderson/Pulay DIIS are robust methods such that if a “bad move” is proposed, it will usually have a large associated error vector, and therefore its DIIS coefficient will be small. 
     1. Anderson-Style Acceleration 
     There are some implementation subtleties regarding the order of Jacobi and Anderson steps, the initiation of the iterative sequence, and the timing and history of the Anderson subspace. For Anderson-style acceleration, the first few iterations of the sequence may be: 
     
       
         
           
             
               
                 
                   
                     
                       θ 
                       → 
                     
                     0 
                   
                   ⁢ 
                   
                     Jacobi 
                     → 
                   
                   ⁢ 
                       
                   
                     
                       θ 
                       → 
                     
                     1 
                   
                   ⁢ 
                   D 
                   ⁢ 
                   I 
                   ⁢ 
                   I 
                   ⁢ 
                   
                     S 
                     [ 
                     
                       
                         
                           θ 
                           → 
                         
                         1 
                       
                       , 
                         
                       
                         
                           
                             
                               
                                 δ 
                                 ⁢ 
                                 
                                   
                                     θ 
                                     → 
                                   
                                   1 
                                 
                               
                               ≡ 
                             
                             → 
                           
                           ⁢ 
                           
                             
                               θ 
                               → 
                             
                             1 
                           
                         
                         - 
                         
                           
                             θ 
                             → 
                           
                           0 
                         
                       
                     
                     ] 
                   
                 
               
               
                 
                   ( 
                   21 
                   ) 
                 
               
             
           
         
       
       
         
           
             
               
                 θ 
                 → 
               
               ′1 
             
             ⁢ 
             
               Jacobi 
               → 
             
             ⁢ 
                 
             
               
                 θ 
                 → 
               
               2 
             
             ⁢ 
             D 
             ⁢ 
             I 
             ⁢ 
             I 
             ⁢ 
             
               S 
               [ 
               
                 
                   
                     θ 
                     → 
                   
                   2 
                 
                 , 
                   
                 
                   
                     
                       
                         
                           δ 
                           ⁢ 
                           
                             
                               θ 
                               → 
                             
                             2 
                           
                         
                         ≡ 
                       
                       → 
                     
                     ⁢ 
                     
                       
                         θ 
                         → 
                       
                       2 
                     
                   
                   - 
                   
                     
                       θ 
                       → 
                     
                     ′1 
                   
                 
               
               ] 
             
           
         
       
       
         
           
             
               
                 θ 
                 → 
               
               ′2 
             
             ⁢ 
             
               Jacobi 
               → 
             
             ⁢ 
                 
             
               
                 θ 
                 → 
               
               3 
             
             ⁢ 
             D 
             ⁢ 
             I 
             ⁢ 
             I 
             ⁢ 
             
               
                 S 
                 [ 
                 
                   
                     
                       θ 
                       → 
                     
                     3 
                   
                   , 
                     
                   
                     
                       
                         
                           
                             δ 
                             ⁢ 
                             
                               
                                 θ 
                                 → 
                               
                               3 
                             
                           
                           ≡ 
                         
                         → 
                       
                       ⁢ 
                       
                         
                           θ 
                           → 
                         
                         3 
                       
                     
                     - 
                     
                       
                         θ 
                         → 
                       
                       ′2 
                     
                   
                 
                 ] 
               
               . 
             
           
         
       
     
     Here DIIS[θ,∈] is a subroutine call that adds state vector θ and error vector {right arrow over (∈)} to the iterative history, and then returns an extrapolated result for {right arrow over (θ)} from the current contends of the iterative history. 
     2. Pulay-Style DIIS Acceleration 
     Similarly, for Pulay-style DIIS the first few iterations of the sequence may be: 
                       θ   →     0     ⁢   D   ⁢   I   ⁢   I   ⁢     S   [             θ   →     0     ,     →     ⁢       G   →       (     θ   0     )         ]     ⁢       θ   →     ′0     ⁢     Jacobi   →             (   22   )                         θ   →     1     ⁢   D   ⁢   I   ⁢   I   ⁢     S   [             θ   →     1     ,     →     ⁢       G   →       (     θ   1     )         ]     ⁢       θ   →     ′1     ⁢     Jacobi   →                     θ   →     2     ⁢   D   ⁢   I   ⁢   I   ⁢     S   [             θ   →     2     ,     →     ⁢       G   →       (     θ   2     )         ]     ⁢       θ   →     ′2     ⁢     Jacobi   →           
Here {right arrow over (G)}(θ) is the gradient vector
 
               ∂   𝒪       ∂     θ   g             
at the current parameter set.
 
     While the above written description enables one of ordinary skill to make and use this method, those of ordinary skill will recognize that variations of the specific embodiment are possible. The patent should therefore not be limited by the particular description above. 
     Although the detailed description contains many specifics, these should not be construed as limiting the scope of the invention but merely as illustrating different examples. It should be appreciated that the scope of the disclosure includes other embodiments not discussed in detail above. Various other modifications, changes and variations which will be apparent to those skilled in the art may be made in the arrangement, operation and details of the method and apparatus disclosed herein without departing from the spirit and scope as defined in the appended claims. Therefore, the scope of the invention should be determined by the appended claims and their legal equivalents.