Method and system for black box optimization using quantum computing

Abstract

The invention relates to a method and system for black box optimization using quantum computing. A computer-implemented method for solving an optimization problem for a given function, said method comprising the steps of: obtaining an approximation of the function in the form of an approximation matrix product operator (MPO) representation of the function; selecting an approximate rank based on the approximate MPO representation; determining a rank-based orthogonal approximation of the approximate MPO representation in the form of an orthogonal MPO having an equidistant sub-tensor of the approximate rank; encoding the orthogonal approximation into a quantum circuit based on encoding the equidistant sub-tensor as a quantum gate; and encoding any power of the orthogonal approximation by applying quantum circuits to the input quantum state multiple times in the form of quantum circuits in a series sequence.

Description

Technical Field

[0001] The present invention belongs to the field of encoding computational operations into quantum circuits. More precisely, the present invention relates to computer-implemented optimization of desired computational operations for encoding into quantum circuits. Background Art

[0002] Quantum computers provide a platform for controllable quantum mechanical systems whose states and interactions can be controlled to perform computations. Computations are performed through the deterministic evolution of controllable quantum mechanical systems (e.g., qubits, which are quantum analogs of classical bits), and the states of quantum mechanical systems can be measured to determine the results of the computations.

[0003] The control operations on these qubits are called quantum gates. Quantum gates can act coherently on qubits to induce changes in the state of a single qubit (so-called single-qubit gates) and on multiple qubits (so-called multi-qubit gates) (for example, in order to entangle the states of multiple qubits) and any combination thereof. For example, a single-qubit gate can induce the spin state of an electron to rotate to a selectable value (for example, π / 2). Multi-qubit gates can act coherently on two or more qubits, such as performing a coherent CNOT operation on the states of two qubits. Multiple quantum gates can be applied to the qubits of a quantum computer in parallel or sequentially to perform calculations. Ultimately, the state of the qubit can be repeatedly measured after applying a series of quantum gates to determine the probability of each possible calculation result.

[0004] Essentially, the operation of a quantum computer can be thought of as relying on encoding a starting state into the internal quantum state of a qubit, then multiplying the internal quantum state of the qubit by a matrix (which can be implemented through a combination of different quantum gates), and measuring the result.

[0005] To compute solutions to problems that are considered intractable on classical computers, quantum computers can exploit the special properties of quantum mechanical states (particularly superposition and entanglement of different quantum states) to find solutions with a relatively small number of computational steps. In addition, since quantum operations can create complex superposition states of qubits, quantum computers can in principle access large internal memories to process computational states.

[0006] However, as Shende et al. discovered ("Minimal universal two-qubit controlled-not-based circuits"; Physical Review A, 69:062321, June 2004), converting an arbitrary operation into a combination of single-qubit gates and two-qubit gates typically requires a number of operations that is exponential in the number of qubits.

[0007] Lubasch et al. ("Variational quantum algorithms for nonlinear problems") teach solving nonlinear differential equations by variational quantum computation. The outputs of quantum variational circuits are combined by quantum nonlinear processing units and compared to an MPS anatomy, where a Hamiltonian representable as MPS is encoded as a quantum circuit according to the MPS quantum circuit encoding scheme. Summary of the invention

[0008] However, current approaches remain limited due to the significant complexity involved in implementing the corresponding steps of encoding the state into qubits, encoding the matrix through quantum gates, and measuring. Specifically, the exponential complexity of factoring arbitrarily multi-qubit gates with respect to the number of qubits, combined with the need for full connectivity between qubits, presents a significant obstacle. Although algorithms exist to convert large matrices into a combination of single-qubit gate operations and two-qubit operations (e.g., CNOT gates), the complexity of the resulting circuits typically scales exponentially with the number of qubits. An example of an applicable problem using the Hamiltonian-limited algorithm that can be expressed as an MPS.

[0009] In particular, higher powers of the desired operation cannot be directly generated using arbitrary tensor network representations. However, higher powers of computational operations can be exploited in certain computational schemes, such as the power method for determining the highest eigenvalue of a matrix. If the matrix encodes an approximation of a black-box function, the eigenvalues ​​can indicate the maximum value of the black-box function, allowing a large class of optimization procedures to be solved using quantum processing resources.

[0010] In view of this prior art, the object of the present invention is to provide a computer implemented method and a corresponding system for determining a quantum coding for an approximation to a function in the form of an encoding matrix and its arbitrary powers, wherein the quantum coding can closely approximate a black box function and its arbitrary powers but can also be implemented on quantum hardware without the depth of the quantum circuit growing exponentially with the number of qubits.

[0011] This object is solved by a method and a system according to the independent claims. The dependent claims relate to preferred embodiments.

[0012] According to a first aspect, the present invention relates to a computer-implemented method for solving a function optimization problem. The method comprises the following steps: obtaining an approximation of a function in the form of an approximate matrix product operator MPO representation of the function, and selecting an approximation rank based on the approximate MPO representation. The method further comprises determining a rank-based orthogonal approximation of the approximate MPO representation in the form of an orthogonal MPO of an isometric subtensor with the approximate rank. The method further comprises: encoding the orthogonal approximation into a quantum circuit based on encoding the isometric subtensor as a quantum gate acting on a computational qubit and a set of initialization auxiliary qubits; and encoding an arbitrary power of the orthogonal approximation based on applying the quantum circuit multiple times to an input quantum state in the form of a series of quantum circuits.

[0013] The function may be a black box function and may be provided with a set of inputs and associated outputs, or may be based on a discretization of a smooth function. The function may typically be mapped to a matrix, such as based on a diagonal matrix, where each entry corresponds to a function value at an input associated with a corresponding computational basis state. In some examples, the function is provided as a matrix or matrix representation.

[0014] Furthermore, given a set of function values ​​for corresponding input values, there are a variety of algorithms for generating an approximate representation that implements the function (such as the TT cross approximation), which can generate a tensor network approximation, typically in MPS or MPO format.

[0015] In a preferred embodiment, the approximate MPO representation is an approximate tensor network representation of the function, wherein the method comprises in particular obtaining a sequence of sub-tensors via cross approximation of the function.

[0016] A tensor network approximation can interpolate and / or approximate a function based on a set of input and corresponding output pairs of the function, such as by sampling the function at a randomly selected set of inputs. For example, a tensor train (TT) approximation of a function can be obtained using a TT cross approximation. Although in principle the dimensionality of the underlying matrix corresponding to the discretized function is large, such a tensor network representation can still be efficiently stored in memory. Preferably, the tensor network approximation has a relatively small rank.

[0017] In a preferred embodiment, the method further comprises obtaining an approximate MPO representation based on a plurality of sample values ​​of the function for a corresponding plurality of different input values.

[0018] For example, a discretized sample values ​​of a function (eg, a multivariate black-box function) for associated inputs may be used to obtain a tensor network approximation and / or an MPS approximation of the function using skeleton decomposition or crossover approximation.

[0019] In a preferred embodiment, obtaining the approximate MPO representation includes obtaining a matrix product-state MPS representation of the function by a cross approximation technique, and obtaining the approximate MPO representation based on the MPS representation by enforcing diagonal restrictions.

[0020] In general, any function can be approximated by first recovering the MPS, which can then be expanded in a straightforward way to obtain an approximate MPO representation. The MPS representation can be converted to MPO by enforcing diagonal restrictions and using the same set of subtensors obtained for the MPS representation of the function.

[0021] The approximate MPO representation of a function can have a relatively small bond dimension. However, the approximate MPO representation / approximation has no direct correspondence with a quantum circuit that can be determined without exponential computational difficulty. Instead, an approximation rank is chosen based on the approximate matrix product operator representation, which should be equal to or greater than the rank of the approximate MPO representation in order to approximate the function using an orthogonal approximation for which more efficient encoding strategies exist.

[0022] In some embodiments, the approximation rank is greater than the rank of a matrix product operator representation of the function, wherein the rank of the matrix product operator representation is specifically the key dimension of the matrix product operator representation.

[0023] Isometric constraints on isometric subtensors may benefit from an increase in rank due to the advantages associated with an increase in degrees of freedom, such as maintaining sufficient approximation accuracy for common problems.

[0024] It is noted that the exact matrix product operator representation may be determined for certain matrices / functions, but the matrix product operator representation may generally be obtained by an approximation algorithm (such as a cross approximation or an algorithm based on singular value decomposition) to approximate the function using MPO without significantly affecting the results of the method. The skilled person will further recognize that the matrix product operator representation / approximation may already be known, and it may be sufficient to receive the matrix product operator representation from a suitable source.

[0025] The approximate MPO representation may include non-isometric subtensors and may be expressed with a certain error (which may be less than 1% or less than 10 -6 , e.g. depending on the application) to approximate the function, and may also be referred to as a non-orthogonal MPO approximation of the function. For example, for solving an optimization problem, an approximation of the function with an error of about 1% may be sufficient. The skilled person will recognize that not every matrix can be represented by an MPO representation with an arbitrarily small accuracy, but the skilled person may still choose to apply this method in view of the increased inaccuracy of the MPO approximation, e.g. if there is no other way to perform the calculation, or if a relatively low accuracy is sufficient for the respective computational problem.

[0026] Based on the bond dimension (rank) of the approximate MPO representation, the approximate rank can be chosen to be the next power of two (i.e., 2 N , where N is an integer), or any power of two that is greater than the bond dimension of the approximate MPO representation. As an example, if the rank of the approximate MPO representation of the function is 3, then the rank chosen for the isometric subtensors of the orthogonal approximation may be 4 or 8. However, if the approximate rank is equal to or less than the rank of the approximate MPO representation of the function, then in principle a suitable approximation can be made, but in this case the orthogonal approximation is more likely to be an inaccurate representation of the function.

[0027] An implementation based on the orthogonal approximation of the approximate MPO representation of the function can obtain the maximum eigenvalue of the function / approximate MPO representation by sequentially applying the quantum circuit encoding of the orthogonal approximation to the input state multiple times in a manner similar to the classical "power method" and measuring the output.

[0028] To determine an orthogonal MPO approximation, the method includes determining an initial guess of an orthogonal MPO corresponding to an MPO having isometric subtensors, the isometric subtensors having a suitable rank. The initial guess may be composed of random subtensors that satisfy an isometry condition.

[0029] In a preferred embodiment, obtaining an orthogonal approximation in the form of an orthogonal MPO includes: determining an initial guess of the orthogonal approximation using isometric subtensors of the approximate rank; and starting from the initial guess, iteratively optimizing the orthogonal approximation based on an optimization algorithm, which minimizes a cost function under the isometry constraint of the isometric subtensors, wherein the cost function assigns a cost to the orthogonal approximation based on the quality of the orthogonal approximation with respect to the approximate MPO representation.

[0030] In some embodiments, the orthogonal approximation is a matrix product operator, and the orthogonal approximation is mathematically equivalent to a tensor network in which each intermediate isometric sub-tensor has two outer uncontracted indices and two inner indices that are contracted with adjacent isometric sub-tensors in a chained manner.

[0031] The unmerged indices may correspond to the dimensions of the input and output qubit states of each isometric subtensor, while the internal indices may be the approximate rank.

[0032] In some embodiments, the approximate rank is the common key dimension of all sub-tensors of the tensor network.

[0033] Subsequently, the initial guess is optimized based on an iterative optimization algorithm so that the difference between the orthogonal approximation and the function approximated by the approximate MPO representation is minimized. The difference can be based on the norm of the difference, such as the Frobenius norm or the cosine norm of the difference between the orthogonal approximation and the approximate MPO representation of the function (and in particular the subtensors of the approximate MPO representation).

[0034] In some embodiments, the cost function is based on a difference function of the function and a renormalized orthogonal approximation evaluated based on tensor network calculus, wherein the renormalized orthogonal approximation is based on the orthogonal approximation multiplied by a renormalization constant, wherein the cost function is specifically based on the Frobenius norm of the difference function evaluated based on tensor network calculus, and / or wherein the renormalization constant is selected such that the approximate MPO representation of the function divided by the normalization constant does not increase the trace of the state on which the expected matrix acts, and / or such that the trace of the density matrix on which the renormalized orthogonal approximation acts is equal to or less than 1, regardless of the initial state.

[0035] The renormalization constant can preserve the constraint that quantum operations cannot increase the trace of the density matrix of the quantum state, that is, given any input density matrix ρ, the orthogonal approximation A should satisfy:

[0036] Given an arbitrary matrix M representing the function, the orthogonal approximation A can be renormalized by the renormalization constant c, so that the cost function of the difference is expressed as C = ||cA-M|| 2 , where ||·|| represents the chosen difference norm, such as the Frobenius norm. In other words, the operation to be implemented in the quantum computer can be the renormalized matrix c representing the function -1 M.

[0037] In some embodiments, for Alternatively, the orthogonal approximation A of the approximate MPO representation is optimized according to a gradient-based optimization algorithm so as to obtain a renormalization constant c at each iterative step of minimizing the cost function.

[0038] For example, the renormalization constant may be updated using steepest descent on c to minimize the cost function.

[0039] Based on a cost function that quantifies the difference between the (renormalized) orthogonal approximation and the approximate MPO representation of the function, elements of the orthogonal approximation may be updated to minimize the cost function.

[0040] In some embodiments, the optimization algorithm comprises gradient descent, and in particular an optimization algorithm based on Riemannian gradient descent or derived therefrom.

[0041] In other words, the method may comprise (stochastically) determining the gradient of the difference between the orthogonal approximation and the approximate MPO representation of the function, and subsequently updating the orthogonal approximation according to the determined gradient.

[0042] The skilled person will recognize that in some embodiments, minimizing the difference between the orthogonal approximation and the approximate matrix product operator representation of the function may be preferred because this can reduce the computational complexity of the optimization. For example, the difference between the subtensor of the orthogonal approximation and the approximate matrix product operator representation of the function can be determined, and the corresponding isometric subtensor can be updated based on the difference. Any update to the isometric subtensor should preserve the isometry constraint, that is, the updated isometric subtensor should be isometric.

[0043] In some embodiments, the orthogonal approximation is a matrix product operator, which can be expressed as

[0044]

[0045] Among them, V k is isometric, that is, for any index k of the isometric subtensor, And where each isometric subtensor V k Specifically, an iterative optimization is performed by alternating descent to optimize the orthogonal approximation.

[0046] The alternating descent may be implemented as a stochastic gradient descent of each sub-tensor, for example in a predetermined or random order, such that each isometric sub-tensor of the orthogonal approximation approaches a sub-tensor of an approximate matrix product operator representation of the function.

[0047] In some embodiments, the gradient of the cost function is determined based on multiple partial derivatives of the sub-tensors of the orthogonal approximation (e.g., by determining the partial derivative of the cost function with respect to the value of each sub-tensor, and by updating the elements of the orthogonal approximation based on the partial derivatives of multiple sub-tensors (e.g., all sub-tensors). The skilled person will understand that the optimization may include stochastic optimization, for example, partial derivatives need not be determined for all parameters of the orthogonal approximation at each iteration, but rather a randomly determined subset of parameters may be updated, such as in stochastic gradient descent and related optimization algorithms.

[0048] However, the updated subtensors should preserve the isometries, for example by updating the isometries according to Riemannian gradient descent. Riemannian gradient descent can preserve the isometries, for example by projecting the update vectors and / or the updated subtensors onto the Stiefel manifold of isometries.

[0049] For example, optimizing the orthogonal approximation under the isometry constraint may include determining the partial derivative of the cost function or the optimization vector obtained thereby at the point V k Stiefel manifolds in m×p isometric matrices Projection on the space tangent

[0050] For example, the subtensor V k An optimization algorithm similar to steepest descent optimization can be used according to The update is performed by , where α is the step size. However, the partial derivatives may also be used initially to determine the search direction, for example based on adaptive moment estimates, and the search direction may be projected onto the space tangent of the Stiefel manifold.

[0051] In some embodiments, optimizing the orthogonal approximation under the isometry constraint includes determining the subtensor V based on partial derivatives of the cost function and optionally based on information from a previous iteration of the optimization algorithm. k The search direction p k ; and determine the search direction p k At point V k Stiefel manifolds with m×p isometric matrices Projection on tangent space Among them, when reshaped into a two-dimensional matrix, the isometric subtensor V k The dimension is m×p.

[0052] The partial derivatives and / or the search direction may then be used to determine an update direction, and the isometric subtensor and / or the orthogonal approximation may be modified based on the update direction. For example, the modified subtensor and / or the modified orthogonal approximation may be retracted onto the Stiefel manifold to obtain an updated subtensor and / or an updated orthogonal approximation.

[0053] In some embodiments, optimizing the orthogonal approximation under the isometry constraint includes: projecting the modified subtensor and / or the modified orthogonal approximation onto a Stiefel manifold of the subtensor and / or the orthogonal approximation, respectively, wherein the modified subtensor and / or the modified orthogonal approximation is modified based on the search direction and / or the update direction; and / or includes determining the subtensor V k The update to the retraction on the manifold To update the sub-tensor V k , where the retraction is based on the update direction The update direction is based on the partial derivative of the cost function with respect to the subtensor V k Encodes the direction of the update.

[0054] Retraction can update the vector Projected onto the product of the Stiefel manifold and the isometric submatrices, for example, according to Thus, isometry constraints are enforced in all iteration steps. For example, the updated intermediate subtensor can be determined based on the projection of the updated gradient on the space tangent of the Stiefel manifold, and the updated intermediate subtensor can be retracted onto the Stiefel manifold.

[0055] The skilled person will appreciate that this technique can be extended to an optimization technique (which may include momentum), such that the update vector can be a function of the projection of the partial derivatives onto the tangent of the Stiefel manifold and the calculation results from the previous iteration step. For example, the iterative optimization can be based on an adaptive moment based optimization algorithm. Therefore, more generally, the update of the isometric subtensors can be based on an update function that is based on the projection of the partial derivatives of the cost function onto the tangent of the Stiefel manifold, wherein the update function determines the intermediate updates; and / or the update of the isometric subtensors can be based on the retraction of the intermediate updates on the Stiefel manifold, so that the resulting subtensors / orthogonal approximations are isometric.

[0056] In some embodiments, optimizing the orthogonal approximation under the isometry constraint includes updating the orthogonal approximation by updating all isometric subtensors, where each isometric subtensor V k All pass through point V k Stiefel manifold of m×p isometric matrices along the update direction The retraction is performed on the update, where, when reshaped into a two-dimensional matrix, the isometric subtensor V k The dimension of is m×p, and the update direction is the direction that minimizes the cost function based on the partial derivative of the cost function, and in particular uses the projection p k,tangent As used to update the isometric subtensor V k Update direction.

[0057] The skilled person will appreciate that different deflations may be used in embodiments, such as the Cayley deflation or the singular value decomposition, to name just two examples, and that the skilled person may select a suitable deflation from a variety of known deflations.

[0058] During iterative optimization of the orthogonal approximation, the isometric subtensors may be reshaped into matrices based on isometric constraints on the subtensors, the matrices being square for intermediate subtensors of the orthogonal approximation and the reshaped matrices corresponding to the intermediate subtensors being unitary.

[0059] In some embodiments, before the orthogonal approximation is expanded, when the isometric subtensors and intermediate tensors of the orthogonal approximation at the ends of the orthogonal approximation are reshaped into a two-dimensional matrix, the isometric subtensors are isometric and the intermediate tensors are unitary, so that after the isometric subtensors are reshaped into a square matrix, all the isometric subtensors are unitary.

[0060] To obtain a set of unitary matrices that can be implemented in a quantum circuit, auxiliary qubits can be used to extend the subtensors at the ends of the tensor chain (sometimes also called boundary cores), for example, by preparing an auxiliary qubit in the zero state |0> to multiply with the first isometric subtensor and by implementing a post-selection at the opposite end of the tensor chain.

[0061] Thus, the method may include expanding the bounding kernel of the orthogonal approximation so that after reshaping the isometric subtensors into square matrices, all isometric subtensors are unitary. Expanding the orthogonal approximation may be part of determining an initial guess for the orthogonal approximation, or may be performed at the beginning of an iterative optimization of the orthogonal approximation or after it has been iteratively optimized. For example, the initial guess may consist of subtensors that may be reshaped into unitary matrices, and the determination of the cost function may take into account an implicit expansion of the orthogonal approximation.

[0062] In some embodiments, the extended orthogonal approximation is based on providing one or more additional auxiliary qubits as input to the first end of the tensor network, wherein the method specifically provides at least log R additional auxiliary qubits, wherein R is the approximation rank.

[0063] In some embodiments, extending the orthogonal approximation includes measuring a qubit at a second end of the tensor network, wherein selecting a result having a predetermined measurement result of the qubit measured at the second end of the tensor network is part of implementing the approximate MPO representation of the function.

[0064] Based on the addition and post-selection of auxiliary qubits, all sub-tensors can be square after appropriate reshaping, and can thus form a set of unitary matrices that can be implemented in quantum circuits by decomposing the isometric sub-tensors into quantum gates of the corresponding quantum architecture (e.g., a combination of single-qubit rotations and CNOT gates).

[0065] The skilled person will recognize that the expansion can be performed before or after the iterative optimization of the orthogonal approximation. In other words, the initial guess can be defined based on the expanded matrix product operator, or the orthogonal approximation can be expanded after it has been iteratively optimized, for example by introducing and filling additional elements of the boundary core so that the boundary core obtained after appropriate reshaping is unitary.

[0066] In some embodiments, the method further comprises implementing the orthogonal approximation of the function as a quantum circuit on quantum hardware based on encoding the isometries as quantum gates.

[0067] The decompositions of the corresponding isometric subtensors can be combined to implement an orthogonal approximation based on sequential application of the corresponding isometric subtensors, where the output state of log R qubits produced by the unitary operation corresponding to one of the subtensors can be used as input to the unitary operation corresponding to the next isometric subtensor in the tensor network.

[0068] In some embodiments, implementing an orthogonal approximation of a function as a quantum circuit is based on sequentially applying isometries to a set of qubits and rearranging unitary matrices into a quantum network.

[0069] As a result, the depth of a quantum circuit can grow only polynomially with the number of qubits, rather than exponentially as is typically done by factoring matrices into single-qubit and two-qubit gates.

[0070] In a preferred embodiment, applying the quantum circuit a plurality of times in the form of a serial sequence of quantum circuits comprises preparing a number of auxiliary qubit groups to an initial state, wherein the number is at least equal to the number of times the quantum circuit is applied.

[0071] Specifically, for each extended orthogonal approximation, corresponding inputs of auxiliary qubits are prepared, and each auxiliary qubit is initialized to a corresponding initial state of the auxiliary qubit, such as a "zero", "one" or "ground" state.

[0072] In a preferred embodiment, each group of auxiliary qubits includes a number of qubits equal to the logarithm of the approximate rank.

[0073] In a preferred embodiment, a subsequent quantum circuit in a serial sequence of quantum circuits operates on a corresponding set of ancillary qubits and a set of computational qubits that were operated on by a previous quantum circuit.

[0074] In other words, the method may begin with initializing a set of initial computational qubits and a first set of auxiliary qubits acted upon by a first quantum circuit in a series of quantum circuits, and a second quantum circuit in the series of quantum circuits may act upon the first set of auxiliary qubits, the second set of auxiliary qubits, and a subset of the initial computational qubits, which subset may be smaller than the initial computational qubits by the number of auxiliary qubits in the second set of auxiliary qubits.

[0075] In a preferred embodiment, each time a quantum circuit is applied, the method comprises measuring a set of post-selected qubits on which the corresponding quantum circuit has acted and on which a subsequent application of the quantum circuit in the serial sequence of quantum circuits will not act.

[0076] The post-selection can be based on the measurement results of the post-selected qubits of each quantum circuit in the series of quantum circuits. For example, when all post-selected qubits are measured to be in the |0> state, the results can be recorded to determine the maximum eigenvalue of the approximate MPO representation of the function.

[0077] In a preferred embodiment, each group of post-selected qubits includes a number of qubits equal to the logarithm of the approximate rank.

[0078] In a preferred embodiment, preparing the initial quantum state comprises applying a Hadamard gate to each qubit in a set of initial computational qubits.

[0079] Applying a Hadamard gate to each qubit in the initial computational qubits prepares a vector of computational ground states "one" to support a quantum power method for determining the largest eigenvalue of a function.

[0080] A quantum computing system implementing the method may include a number of initial computational qubits equal to the number of isometries, wherein, in particular, each quantum operation of one of the isometries implementing a first quantum circuit in a series of quantum operations may act on a corresponding one of the initial computational qubits. A first quantum operation in a quantum operation cascade implementing the quantum circuit may act on the set of auxiliary qubits. Subsequent quantum operations in the cascade may act on a corresponding one of the computational qubits and on a number of processed qubits equal to the logarithm of the approximate rank, wherein a previous quantum operation in the quantum operation cascade has acted on the processed qubit.

[0081] Each quantum operation may have as output a quantum state of a qubit that will not be acted upon by subsequent quantum operations in the quantum operation cascade. The output of the quantum operation may form the computational qubit of the next quantum circuit in the series of quantum circuits. The last quantum operation in the quantum operation cascade may include as output a quantum state for processing by subsequent quantum circuits and a number of processed qubits equal to the logarithm of the approximate rank, which may be measured as post-selected qubits.

[0082] Those skilled in the art will appreciate that the output of each quantum operation that is not acted upon by subsequent quantum operations in the cascade may include more than one qubit, and a corresponding number of initial computational qubits may be prepared for each quantum operation of the first circuit. However, in the following, reference will be made primarily to an example in which each quantum operation has a single initial computational qubit.

[0083] In a preferred embodiment, applying the quantum circuit p times to the input quantum state in the form of a serial sequence of quantum circuits includes preparing p groups of auxiliary qubits, wherein each group of auxiliary qubits is acted upon by a respective one of the quantum circuits and includes a number of auxiliary qubits equal to the logarithm of the approximate rank.

[0084] Each quantum circuit after the first in the series of quantum circuits may act on a corresponding set of those sets of auxiliary qubits and on processed computational qubits acted upon by a previous quantum circuit in the series, wherein the processed computational qubits did not form part of the post-selected qubits of the previous quantum circuit.

[0085] According to a second aspect, the present invention relates to a processing system for solving an optimization problem. The system is configured to obtain an approximation of a function in the form of an approximate matrix product operator MPO representation of the function, and to select an approximation rank based on the approximate MPO representation. The system is further configured to determine a rank-based orthogonal approximation of the approximate MPO representation in the form of an orthogonal MPO of an isometric subtensor with an approximate rank. The system is further configured to: determine an implementation of an orthogonal approximation to an arbitrary power in a quantum circuit based on encoding the isometric subtensor as a quantum gate acting on a computational qubit and a set of initialization auxiliary qubits, wherein in the implementation, the quantum circuit is applied multiple times to an input quantum state in the form of a series of quantum circuits.

[0086] The system may implement the method according to the first aspect or any combination of the embodiments of the first aspect. In particular, the system according to the second aspect may also benefit from any features of the preferred embodiments of the first aspect.

[0087] The processing system may include a single processing unit, or may include multiple processing units that may be functionally connected. The processing unit may include a microcontroller, an ASIC, a PLA (CPLA), an FPGA, or other processing device (including a processing device that operates based on software, hardware, firmware, or a combination thereof). The processing device may include an integrated memory, or communicate with an external memory, or both, and may further include an interface for connecting to a sensor, device, apparatus, integrated logic circuit, other controller, etc., wherein these interfaces may be configured to receive or send signals such as electrical signals, optical signals, wireless signals, acoustic signals, etc.

[0088] The system may be configured to obtain the approximate MPO representation, for example by executing a cross-approximation or decomposition algorithm on the function, or by receiving the approximate MPO representation from a suitable source, such as a database.

[0089] In a preferred embodiment, in order to obtain an orthogonal approximation in the form of an orthogonal MPO, the system is configured to: determine an initial guess for an orthogonal approximation of a function in the form of a tensor network of isometric subtensors having approximate rank; and starting from the initial guess, iteratively optimize the orthogonal approximation based on an optimization algorithm that minimizes a cost function under the isometry constraint of the isometric subtensors, wherein the cost function assigns a cost to the orthogonal approximation based on the quality of the orthogonal approximation with respect to the approximate MPO representation.

[0090] In some embodiments, the system is further configured to extend the orthogonal approximation so that after reshaping the isometric subtensors into square matrices, all isometric subtensors are unitary.

[0091] In a preferred embodiment, the system is further configured to transmit the implementation to a quantum computing system to execute the implementation on quantum hardware.

[0092] This embodiment may allow determining the largest element of an orthogonal approximation of a function, in particular the largest eigenvalue.

[0093] According to a third aspect, the present invention relates to a hybrid quantum-classical computing system, the hybrid quantum-classical computing system comprising the system according to the second aspect and quantum computing hardware. The hybrid quantum-classical computing system is configured to receive the embodiment from the processing system, implement the embodiment in the quantum computing hardware, and receive the calculation result from the quantum computing hardware.

[0094] According to a fourth aspect, the invention relates to a non-transitory medium comprising machine-readable instructions which, when executed by a processing system, implement the method according to the first aspect and / or the system according to the second and / or third aspect. BRIEF DESCRIPTION OF THE DRAWINGS

[0095] The features and numerous advantages of the method and system according to the present invention will be best understood from the detailed description of the preferred embodiments with reference to the accompanying drawings, in which:

[0096] Figure 1 An example of a quantum computing system is schematically illustrated;

[0097] Figure 2 An example of a computer-implemented method for encoding higher powers of a function approximation into a quantum circuit is illustrated;

[0098] Figure 3 schematically illustrates a flow chart for determining a matrix product operator representation based on a given expected matrix;

[0099] Figure 4 schematically illustrates a strategy for determining a quantum circuit encoding an approximation of an expected matrix according to an example;

[0100] Figure 5 and Figure 6 The diagram shows the combination Figure 4 The results obtained when the discussed encoding strategy is applied to two different expectation matrices; and

[0101] Fig. 7A and Figure 7B illustrates experimental runtime results from implementing a method for determining a quantum circuit when approximately diagonal matrices according to an example;

[0102] FIG. 8A to FIG. 8F illustrates a schematic flow chart of a quantum computer implementation of an algorithm for determining a maximum eigenvalue of a matrix that approximates a function according to an example;

[0103] Fig.9A and Fig. 9B illustrates a detailed implementation of a quantum power method for finding the maximum eigenvalue of a matrix of an approximate function based on quantum circuit encoding of computational operations according to corresponding examples; and

[0104] FIG. 10A to FIG. 10D Examples of measurement probabilities obtained for quantum circuit encoding and the simulated quantum power method for two-dimensional Ackley functions are illustrated in the form of contour plots. DETAILED DESCRIPTION

[0105] Figure 1 An example of a quantum computing system 10 is schematically illustrated. The system 10 includes a qubit register 12 including computational qubits, a state preparation module 14 for preparing the quantum state of the computational qubits to a predetermined quantum state, a quantum matrix multiplication module 16 for transforming the predetermined quantum state according to a multi-qubit operation, and a measurement module 18 for measuring the result of the quantum computation, for example, by performing a projection measurement on the state of some or all of the qubits. The processing system 20 may initialize the computational qubits in the qubit register 12, define the operations of the state preparation module 14 and / or the quantum matrix multiplication module 16, and may obtain the measurement results from the measurement module 18, for example, by determining a control operation and sending the control operation to the hardware of the quantum computing system 10.

[0106] The state preparation module 14 and the quantum matrix multiplication module 16 may be similar in construction, and may implement a combination of single-qubit quantum gates and multi-qubit quantum gates to affect the state of the qubits of the qubit register 12. It should be noted that in some embodiments, the state preparation module 14 may not be required, and the combination of quantum gates determined by the quantum matrix multiplication module 16 may directly act on the quantum state of the qubits initialized to a certain initialization quantum state in the qubit register 12. The type of quantum gate implemented in the quantum circuit by the state preparation module 14 and / or the quantum matrix multiplication module 16 may depend on the architecture of the specific quantum hardware used, such as a native gate based on quantum hardware. The state preparation module 14 and the quantum matrix multiplication module 16 may calculate the quantum gate arrangement and / or send control operations to control the quantum hardware so that the quantum hardware implements the desired computational operation on a certain quantum state.

[0107] However, decomposing arbitrary computational operations into quantum circuits is generally a difficult computational problem, whose complexity usually grows exponentially with the number of qubits involved. In addition, for general operations, the depth of the circuit (i.e., the number of consecutive quantum gates used to implement the quantum operation) may also grow exponentially with the number of qubits. Although it is known that for certain types of problems, it is possible to efficiently encode problems as combinations of quantum gates, encoding many desired problem types into quantum circuits may still be a significant obstacle to achieving quantum computing.

[0108] Figure 2 An example of a computer-implemented method for encoding arbitrary powers of a matrix product operator approximation of a function into a quantum circuit is illustrated. The method includes obtaining an approximate matrix product operator representation of the function (S10), and selecting an approximate rank based on the MPO representation (S12). The method further includes determining a rank-based orthogonal approximation of the approximate MPO representation in the form of an orthogonal MPO of isometric subtensors having an approximate rank. The method further includes: encoding the orthogonal approximation into the quantum circuit based on encoding the isometric subtensors as quantum gates (S16); and encoding arbitrary powers of the orthogonal approximation based on applying the quantum circuit multiple times to an input quantum state in the form of a serial sequence of quantum circuits (S18).

[0109] Determining the orthogonal approximation may include determining an initial guess for an orthogonal approximation of the function in the form of a tensor network of isometric subtensors having approximate rank, and starting from the initial guess, iteratively optimizing the orthogonal approximation based on an optimization algorithm that minimizes a cost function under isometry constraints of the isometric subtensors.

[0110] The cost function may assign a cost to an orthogonal approximation of a function based on the quality of the orthogonal approximation with respect to the approximate MPO representation. Thus, instead of directly implementing the approximate MPO representation in a quantum circuit, the method comprises approximating the approximate MPO representation using an orthogonal approximation constructed in the form of a tensor network with isometric subtensors (hereinafter also referred to as isometric cores).

[0111] The isometry constraints on the subtensors of the orthogonal approximation can be enforced at each step of the iterative optimization (e.g., based on a modified gradient descent method and / or scaling back any updated subtensors to isometry subtensors) to cause the optimization algorithm to converge toward an approximation of the function with isometry subtensors. In general, the approximation rank should be at least equal to the rank of the approximation matrix product operator representation of the function, so that the constraint of isometry subtensors can be at least partially compensated by greater freedom in the merge indexing of the orthogonal approximation.

[0112] For any subtensor between the ends of the tensor network (the so-called boundary cores), the isometric subtensors can be reshaped into square matrices, which can be unitary after appropriate reshaping. These subtensors that can be reshaped into unitary matrices can be directly implemented as quantum circuits using known decomposition algorithms, which can be efficiently computed due to the small size of the reshaped matrices / isometric subtensors. The skilled person will recognize that the subtensors can be implemented without explicitly reshaping the isometric subtensors into matrices.

[0113] To implement boundary cores (i.e., isometric subtensors arranged at the ends of a tensor network in the form of tensor columns), the ends of the tensor network representing the orthogonal approximation can be extended so that all isometric subtensors can be reshaped into unitary matrices, for example by providing auxiliary qubits at one end of the tensor network and introducing measurements with post-selection at the other end so that the boundary cores can likewise be reshaped into square matrices. This extension can be pre-configured as part of an initial guess for the orthogonal approximation, or can be configured when the orthogonal approximation has been iteratively optimized, for example by completing the boundary cores with additional indices so that they are square matrices when the appropriate reshaping is applied, and filling in any missing matrix elements so that the resulting matrix is ​​unitary.

[0114] Therefore, orthogonal approximations can be implemented in quantum circuits based on the decomposition of isometric subtensors that require less memory than the matrix representation of the full function. As a result, the circuit depth of a quantum circuit for the full function can grow only polynomially with the number of qubits, allowing for more complex computations on near-term quantum hardware, where the number of qubits determines the size of the memory space used to hold the computational state.

[0115] Figure 3Schematically illustrates a method for generating a matrix based on a matrix 24 (shown as a square matrix M) according to an example. NxN ) and a flowchart for determining a matrix product operator representation 22 for a given function provided by a Penrose graphical representation 24a. The matrix 24 may correspond to a discretized function value of a function to be implemented on quantum hardware with a computational basis state of a qubit register 12 having n qubits, where N=2 n .

[0116] The matrix 24 associated with the function may be represented in a Penrose diagram representation 24a by a circle with two unmerged N-dimensional indices. The function 24 may be approximated or represented by a matrix product operator representation 22 consisting of a plurality of sub-tensors 26, where summing over the merged indices recovers the individual elements of the matrix 24 representing the original function, for example according to:

[0117]

[0118] Among them, matrix 24A l m The MPO representation 22 of can be composed of log(N) sub-tensors 26, where the unmerged index is l k and m k , and the merge index is j k 、j k-1 .exist Figure 3 In FIG. 2 , each subtensor 26 indicated by a circle has four indices, except for the boundary core at the opposite end of the MPO representation 22 (which has three indices), where the merged indices (extending horizontally) are the MPO bond indices with a bond dimension of r. In the example, the unmerged indices have the dimension of the quantum state of one qubit, i.e., 2.

[0119] Such an MPO representation 22 can be found for a matrix 24 or an approximation thereof using a variety of known algorithms. In principle, any matrix 24 can be represented in a matrix product operator (MPO) representation 22 by decomposing the matrix 24 into a set of matrix products whose characteristic is that their matrix product is equal to a particular element of the matrix 24.

[0120] If only a few function values ​​of a black-box function are known, a matrix product operator / state approximation of the function can be obtained by interpolation, for example using a cross approximation technique (such as skeleton decomposition or TT cross approximation) to obtain an approximate tensor network of the function based on the known function values. The approximate tensor network can correspond to an MPO or MPS approximation of the function.

[0121] However, while MPS states can be efficiently encoded into quantum circuits (e.g. as shown by Lubasch et al.), no corresponding strategy exists for general MPO.

[0122] despite this, Figure 2 The method shown in the example of can be used to provide a quantum circuit for the approximation of a function.

[0123] Figure 4 Schematically illustrates a strategy for determining a quantum circuit 28 encoding an approximation of a function according to an example. The function may initially be represented by an approximate MPO representation 22 having log(N) (i.e., the logarithm of the order of the square matrix to be represented) (or a fraction thereof) of rank r subtensors 26, which may be obtained based on a suitable decomposition and / or approximation algorithm.

[0124] The approximate rank R is selected based on the rank r of the approximate MPO representation 22, which should be equal to or greater than the rank r of the approximate MPO representation 22 and can be a power of two, i.e., 2 k , where k is an integer. Based on the selected approximation rank R, the method can generate an initial guess for the orthogonal approximation 30 using the isometric subtensors 32, 32a, 32b so that all isometric subtensors 32, 32a, 32b are isometric when reshaped into a two-dimensional matrix, where the bond dimension of the orthogonal approximation 30 is equal to the approximation rank R. The isometric subtensors 32, 32a, 32b used for the initial guess of the orthogonal approximation 30 can be randomly generated, for example, selected from random isometric subtensors 32, 32a, 32b.

[0125] Starting from an initial guess of the orthogonal approximation 30, the method can proceed by iteratively optimizing the orthogonal approximation 30 based on an optimization algorithm (whose goal is to minimize the distance metric with the approximate MPO representation 22), wherein the intermediate step size of the optimization algorithm can maintain the isometry of the sub-tensors 32, 32a, 32b.

[0126] Since the approximate MPO representation 22 may not be unitary when reshaped into a matrix, it cannot be directly implemented on a quantum computer. In addition, the quantum operation A will have to preserve the trace

[0127]

[0128] Therefore, it would be preferable to minimize the difference metric with respect to the renormalized orthogonal approximation 30, for example according to

[0129] C=||cA-M|| 2 , (3)

[0130] So that the inequality of equation (2) is satisfied. Therefore, the renormalized approximate MPO representation (e.g., c -1M), rather than the approximate MPO representation 22 (i.e., M), and A indicates the orthogonal approximation 30. From a practical point of view, some calculations may require remembering this normalization factor c and taking it into account at the end of the calculation, for example to interpret the results when they depend on complex output quantum states.

[0131] For the optimization algorithm, the task of the processing system (such as a computer) may be to perform the multiplication of all sub-tensors 32, 32a, 32b (abbreviated as V k ) finds the minimum value of the cost function under the isometry constraint, that is,

[0132]

[0133] Wherein, ||*|| represents a suitable norm, such as the Frobenius norm, and wherein both c and A take optimal values. The optimization can be implemented as an alternating descent by determining the partial derivatives of the cost function with respect to any isometric subtensors 32, 32a, 32b under an isometry constraint, and by optimizing the orthogonal approximation 30 based on the partial derivatives, such as by a gradient descent-based algorithm, such as stochastic gradient descent or a momentum-based algorithm (such as adaptive moment estimation).

[0134] The isometry constraint may be obeyed by ensuring the isometry of the sub-tensors 32, 32a, 32b in each update step of the optimization algorithm and / or by considering the isometry constraint as part of the gradient direction.

[0135] In general, any m-by-p isometric matrix belongs to the Stiefel manifold

[0136]

[0137] This is a Riemannian manifold. Therefore, the cost function C to be minimized can be defined on the Cartesian product of n Stiefel manifolds. In order to find a set of isometric matrices V that minimize this cost function k , we can use gradient-based optimization methods, such as the Riemann gradient descent method, which can solve the optimization problem under constraints in this case.

[0138] In conventional gradient-based methods (such as gradient descent), the variable being optimized is updated by moving in the opposite direction of the gradient of the cost function:

[0139]

[0140] where α is the step size. In the Riemann generalization of these methods, to determine the direction of this movement, one can use its position at point V k The projection δ k RC, without using δ k C. Further, in order to update the matrix V based on the partial derivatives and / or their projections k The value of can be used to shrink on the manifold:

[0141]

[0142] In optimization methods that use information from previous optimization steps (such as gradient descent with momentum), the concept of vector transfer is also used, which is a generalization of the concept of parallel transfer and can therefore be implemented directly based on the above-mentioned scheme.

[0143] The Riemann formula for gradient descent can be implemented as an iterative procedure that computes at each step the Riemann gradient of the optimization function C at the current point on the manifold, and can then find the next point using a retraction chosen in the direction of the negative gradient.

[0144] The algorithm can then be implemented as follows:

[0145] Step 0: Get an initial guess of the orthogonal approximation 30;

[0146] Step 1: Find a renormalization constant c that minimizes the cost function of, for example, equation (3), where the orthogonal approximation 30 is fixed (ie, A is constant in this step);

[0147] Step 2: For the problem of equation (4), find a new approximation A via a step of Riemann gradient descent, where the renormalization constant c can be constant during this step, and then return to step 1.

[0148] The iterative algorithm may terminate at any point, such as after a predetermined number of iterations, or when the cost function has reached a predetermined cost threshold.

[0149] Once the optimized orthogonal approximation 30 is obtained (which may be an orthogonal approximation 30 of a renormalized function, c -1 M), all internal sub-tensors 32 of the orthogonal approximation 30 should already be unitary matrices after corresponding reshaping into two-dimensional matrices.

[0150] However, the sub-tensors 32a, 32b at the boundary (i.e., V 1 and V n, also called boundary cores) can only be isometric. However, the boundary cores 32a, 32b can be expanded to unitary matrices by expanding the subtensors under an isometry condition with an addition of a |0> state 36 at one end of the tensor network and a projection onto a |0> state 38 at the other end of the tensor network to maintain the correspondence with the approximate MPO representation 22. By expanding the orthogonal approximation 30, an expanded orthogonal approximation 34 can be obtained, which is characterized only by having subtensors 32 that can be reshaped as unitary matrices.

[0151] Since all sub-tensors 32 can be reshaped into unitary matrices, they can be obtained by performing a quantum operation 40 (implemented as a unitary operation U 1 toU n ) in a cascade of isometric subtensors 32 to implement the conversion to quantum circuits, such as Figure 4 As shown in the quantum circuit 28 in the example. Here, the unmerged indices of the isometric subtensor 32 may correspond to the input states and output states of the corresponding quantum operation 40, while the merged indices of the isometric subtensor 32 correspond to the quantum state that is the input of the next unitary operation 40 in the serial chain of unitary operations 40. Therefore, each unitary operation 40 has 1+log(R) qubit inputs and qubit outputs, and log(R) qubits are processed by the subsequent unitary operation 40 in the serial chain of unitary operations 40. The decomposition of these subtensors 32 into quantum gates can be performed based on known decomposition techniques for the corresponding quantum hardware architecture.

[0152] The upper |0> state 36 as the input of the upper boundary core 32a of the extended orthogonal approximation 34 may correspond to preparing the additional auxiliary qubit 42 to the "zero" state. At the same time, the lower |0> state 38 as the input of the upper boundary core 32a of the extended orthogonal approximation 34 may correspond to a projection onto the zero state, i.e., a measurement with post-selection.

[0153] The skilled person will recognize that the actual quantum states prepared or measured for these operations can be changed, as long as the correspondence with the approximate MPO representation 22 is maintained, such as by preparing and post-selecting based on the |1> state or a different state, based on an appropriate basis change. However, for the sake of simplicity, the preparation of the "|0>" state and the corresponding post-selection will be considered below.

[0154] Therefore, if Figure 4 As shown in the example of FIG. 2 , the quantum circuit 28 can be implemented by the following steps of the quantum processing system:

[0155] Step 1: Prepare auxiliary qubit 42 to zero state The remaining qubits remain "branch open" (solid horizontal line on the left side of quantum circuit 28), ie the qubits will be encoded according to some initial state based on the input to the computation.

[0156] Step 2: Apply gate U sequentially 1 ,...,U n , these gates should be decomposed into quantum gate implementations based on the quantum hardware employed, such as a combination of single-qubit operations and two-qubit operations.

[0157] Step 3: Measure the auxiliary quantum bit 44 in the calculation ground state, and if all the measurement results in step 3 are equal to "0", record the results, otherwise return to step 1.

[0158] Note that for post-selection, the qubits originally introduced as auxiliary qubits 42 into the initial boundary core 32a should not be measured, but rather the qubits from the other end of the orthogonal approximation 30 (i.e., corresponding to the opposite boundary core 32b (see Figure 4 ))'s quantum bit.

[0159] Only when all auxiliary qubit measurements 44 are equal to |0> (e.g., a string of zero bits for the corresponding qubits is obtained), the algorithm implements the quantum circuit 28 corresponding to the approximate MPO representation 22. The state of the quantum system before the measurement determines the probability of this "successful" outcome. This means that the initial state of the system to which the quantum scheme is applied affects this probability.

[0160]

[0161] The initial state of the system is in It can be determined by the specific computing task. The probability values ​​of all possible initial states of the system can be averaged. <Pr suc > to ignore this dependence. To this end, ρ in can be viewed as a maximally mixed state, i.e., a uniform mixture of states from an orthogonal basis. In this case, ρ in =1 / 2 n I, using <0| a2 U|0> a1 =A, and the average probability can be obtained as

[0162]

[0163] That is, the success probability can be mainly determined by the norm of the orthogonal approximation 30. At the same time, since A is approximately equal to c -1 M, we can see how the success probability is related to the normalization coefficient c, i.e. the higher the normalization coefficient, the lower the success probability. In other words, there may be a trade-off between the accuracy of the approximation and the success probability. It is found that the optimal theoretical value of the normalization constant c is c = λ max(M), which is equal to the maximum eigenvalue of the approximate MPO representation 22. As long as the approximation rank R of the orthogonal approximation 30 is large enough, the function can be approximated accurately enough.

[0164] Note that the depth of quantum circuit 28 can grow polynomially with the approximation rank R, and the size of the approximation rank R can further increase the complexity of decomposing the isometric subtensor 32 into a combination of quantum gates, since each unitary operation U k Both work on log 2 R+1 qubits, where R is the rank of the orthogonal approximation A.

[0165] Note also that the numerical strategy of approximating the decomposition of general multi-qubit unitary operations using a combination of single-qubit and CNOT gates can achieve a number of CNOT gates close to the theoretical lower limit proposed by Shende et al. (Minimal universal two-qubit controlled-not-based circuits, Physical Review A, 2004), i.e., 4 for any decomposition of m qubits. m-1 CNOT gates, so in this case, a total of n*4 log R+1 CNOT gates, which is about n*R 2 Therefore, in general, the lower the approximation rank R selected for the orthogonal approximation 30, the simpler the quantum circuit 28 is, but the less accurate the orthogonal approximation 30 is.

[0166] To demonstrate the capability of this approach, we consider commonly used matrices and encode them into multi-qubit quantum circuits28. Specifically, we consider D and the Laplace operator M L The matrix with discrete linear functions (corresponding to the discretized second-order derivative matrix) on the y-axis, neither of which is unitary.

[0167]

[0168] Wherein, in this case, n is the number of columns N minus one, i.e., N-1.

[0169] Figure 5 and Figure 6 The diagram shows the combination Figure 2 and Figure 4 The algorithm discussed is applied to the diagonal matrix MD and the Laplacian matrix ML of equation (1) respectively. For any n>1, the MPO form of the MD and ML matrices has a rank r equal to 2 and 3 respectively.

[0170] The upper diagram illustrates the relationship between the diagonal matrix ( Figure 5 ) and the Laplacian matrix ( Figure 6 ) as a function of the approximation rank R (in the figure, n is the number of qubits) for the number of auxiliary qubits 42 used equal to log2 R. The lower graph illustrates the corresponding average probability of success based on the auxiliary qubit measurement 44. The dashed line represents the value of ||M||2λmax-2(M)2-n (see equation (9)), which for both matrices is independent of n when the number of qubits n is large.

[0171] For these two matrices, the method can achieve reasonable accuracy (M D and M L The errors of about 0.01% and 0.5% for qubits and 50 qubits, respectively, and a sufficient success probability (over 5%) up to 50 qubits are not a limitation of the method. The important property found in the experiments is that the errors and success probability are almost independent of the number of qubits, making the method very likely to be applicable to more complex tasks.

[0172] To assess computational efficiency, the algorithm can be considered to consist of three main steps, namely, finding an orthogonal approximation 30, expanding the boundary core 32a, 32b to allow reshaping into a unitary matrix, and decomposing the resulting 2R*2R unitary matrix into a series of single-qubit and two-qubit gates.

[0173] The time required for the last step depends on the decomposition method used and will not be discussed here. In addition, in the example, the expansion step is skipped and the value corresponding to U is directly optimized in an iterative optimization algorithm based on an initial guess of the orthogonal approximation 30. 1 and U n The expanded sub-tensor 32 instead of the optimized kernel V 1 and V n itself.

[0174] Therefore, the complexity of the algorithm used to determine the quantum circuit 28 can be mainly determined by the complexity of the iterative optimization process. The execution time of the optimization program is generally proportional to the number of iterations and the execution time of each iteration. The number of iterations required for the optimization algorithm to converge may depend on multiple factors, including the required accuracy, the number of optimization parameters, the selected learning rate and the optimization starting point. Therefore, it may be difficult to determine a sufficient number of iterations in advance, and the following will focus on the complexity of one iteration of the optimization program consisting of the following elements:

[0175] Element 1: The calculation of the normalization constant c can have O(nR 3 ), which is based on scaling the rank r1 and r 2 The computational complexity of multiplying two MPOs and taking their traces at the same time is O(nr 1 r 2 max(r 1 ,r 2 ))$).

[0176] Element 2: The computation of the cost function can have O(n(r+R) 3 ). After summing two MPOs, their ranks are added and then the resulting r 1 +r 2 The rank MPO is combined with its complex conjugate, and its computational complexity is O(n(r 1 +r 2 ) 3 ).

[0177] Element 3: Use automatic differentiation to evaluate the derivative of the cost function. When automatic differentiation is used, its algorithmic complexity is theoretically guaranteed to be no greater than the algorithmic complexity of the original program, which is equivalent to the previous element (2) and is equal to O(n R 3 ).

[0178] Element 4: The optimization process includes a Riemann optimization step, which may involve for each V k Compute the Riemann gradient, retraction, and vector transfer. For a 2R*2R matrix, the first and third steps of the Riemann optimization procedure have O(R 3 ) is calculated with the same complexity. The complexity of the second step depends on the selected shrinking method. In the example implementation, SVD shrinking can be used, and for a 2R*2R matrix, its complexity is also O(R 3 ). Therefore, the Riemann optimization step per iteration can have O(n R 3 ) of the image.

[0179] Fig. 7A and Figure 7B The diagonal matrix M according to the example approximation equation (10) is shown. D When implementing the combination Figure 2 and Figure 4 Experimental runtime results of the discussed method for determining quantum circuit 28.

[0180] Fig. 7A The execution time of one iteration of the iterative optimization algorithm for different approximation ranks R of the orthogonal approximation 30 is illustrated as a function of the number of qubits n (which depends on the size of the matrix 24 ).

[0181] Figure 7BThe experimentally determined execution time of each iteration and its elements (1)-(4) (as discussed in the previous section) is plotted as a function of the number of qubits n.

[0182] In full agreement with theoretical estimates of the complexity of the algorithm, these running times are polynomial in the number of qubits n and therefore in the target matrix size N=2 for the matrix 24. n A multi-logarithmic relationship.

[0183] Therefore, the method for determining quantum circuit 28 is expected to provide an efficient algorithm for determining the quantum gate decomposition of matrices in quantum computing applications. It is noted that the depth of quantum circuit 28 grows only linearly with the number of qubits n, and therefore logarithmically with the matrix size of matrix 24.

[0184] The algorithm can exploit the relative strengths of quantum and classical hardware to generate an optimal algorithm for a given computation. Specifically, an orthogonal approximation 30 can be efficiently optimized on classical hardware with respect to a tensor network representation / approximation 22 of the matrix 24, even if the matrix 24 is large. At the same time, the resulting orthogonal approximation, which consists of compact sub-operations, can be efficiently implemented on quantum hardware where the circuit depth can grow only linearly with the number of qubits.

[0185] In addition, the orthogonal approximation 30 encoded as the quantum circuit 28 implements quantum operations that can be applied to the initial quantum state multiple times to implement higher orders of the approximate MPO representation 22, such as finding the largest eigenvalue of the matrix of the approximate function by measuring the output state of the state after the higher orders of the approximate MPO representation 22 have acted on the input state.

[0186] The task of finding the eigenvalues ​​of a diagonalizable matrix can be viewed as finding the largest element in a tensor (a tensor of values ​​on the diagonal of a matrix), and thus the power method can be applied to almost any optimization problem after corresponding preprocessing, in which the function representation is processed to obtain an MPS / MPO approximation.

[0187] Let f(x) denote some real-valued function that we wish to maximize / optimize over the input domain (x∈[a,b]). Given access to n qubits, we can use 2 n discrete points to represent f(x) as a quantum state,

[0188]

[0189] where |k> is the computational ground state of n qubits, and normalization forces |<f|f> | 2 =1. It can be seen that the maximization problem becomes

[0190]

[0191] Among them, p k =|<k|f> | 2 is the measurement probability of the state |f> with respect to the projection onto the computational basis. Therefore, if such a state can be prepared, measuring it will allow to recover an approximately optimal value of f(x) - up to the error introduced in the discretization. Power iteration can improve this idea. Consider the diagonal matrix with the function values ​​as diagonal elements

[0192]

[0193] The power iteration is as follows:

[0194]

[0195] where ||.|| is the 2-norm. Ideally, when j-→∞, the resulting state becomes |k * >, k * = arg max k p k , thereby recovering the optimal value. In order to actually implement this simple method on a quantum computer, a quantum method can be used to implement the power iteration as in equation (13) (hereinafter also referred to as the "quantum power method").

[0196] In classical computers, multiplication of matrices with large dimensions can be computationally very expensive. However, with quantum computers, the size of the matrices can be compressed exponentially due to the larger computational state space, making this approach feasible for optimization methods.

[0197] FIG. 8A to FIG. 8F A schematic flow chart of a quantum computer implementation of an algorithm for determining a maximum eigenvalue of a matrix that approximates a function according to an example is illustrated.

[0198] exist Fig. 8A In the first step shown, the corresponding input value x of the function f(x) is determined l 、x r The function value f 0 ,…,f 2 n -1 , which is usually a black-box function. In the example shown, the function values ​​are sampled from a smooth curve at equally spaced x values.

[0199] Based on the function value, the MPS state is determined according to a suitable algorithm (e.g. based on a repeated truncated SVD procedure or a crossover approximation). Figure 8B As shown, the function values ​​of the equidistant sampling function can form a function value vector, based on which the value of the MPS state vector can be determined. The MPS state can represent a vector with at most 2nr2 2 of the function value vector of elements n value.

[0200] However, the MPS approximation cannot be directly used in the classical algorithm implementing the power method, because the rank of the MPS will grow exponentially with the power of the power method. Therefore, since the MPS rank grows exponentially with iterations, a simple classical implementation of the power method quickly goes out of scope.

[0201] Nevertheless, based on the MPS representation of a function (which can approximate a function based on discretized function values) (or the function itself), an approximate MPO representation of the function can be obtained22, as Figure 8C Shown schematically.

[0202] Use the combination above Figures 2 to 4 The discussed optimization then allows to obtain an (extended) orthogonal approximation 30 of the approximate MPO representation 22, as Figure 8C , Fig.8D Indicated in.

[0203] for Fig.8D The orthogonal approximation 30 shown in , determines the corresponding quantum circuit 28, which can be composed of the unitary operator 40 for each isometric subtensor 32 obtained by the optimization algorithm.

[0204] To implement the quantum power algorithm, quantum circuit 28 acts multiple times on the input states of the qubits in qubit register 12, such as Fig. 8E A schematic example is shown.

[0205] The input state of the qubit can be obtained based on the state initialization quantum circuit 46, which can include applying a Hadamard gate to each qubit starting from the |0> ground state, thereby initializing a "vector of one" as the computational basis spanned by the qubit. Then, the quantum circuit 28 acts on the qubit multiple times to obtain the final result state, which can be measured using the detector of the measurement module 18. By repeating the quantum circuit 28, quantum power iteration can be achieved.

[0206] like Figure 8F As schematically shown in , the obtained result 48 is measured (eg, measured a plurality of times) to obtain the output of the quantum power algorithm in the form of the maximum eigenvalue of the effective matrix implemented by the quantum circuit 28 based on the orthogonal approximation 30 .

[0207] like Fig.9AAs shown in more detail in the example of , for each power p of the power method, the algorithm implementation includes one quantum circuit 28a-28p applied to the quantum bits in sequence. In addition, for each of the quantum circuits 28a-28p, a corresponding auxiliary quantum bit 42 is initialized and introduced. In the illustration, a single |0> state is shown, but the skilled person will understand that for each of the quantum circuits 28a-28p, log(R) auxiliary quantum bits 42 can be prepared, as discussed above.

[0208] In some examples, where sufficiently fast measurement and initialization techniques may be available, the measured qubit 44 of each of the quantum circuits 28a-28p may be reinitialized and used as an auxiliary qubit 42 for a subsequent quantum circuit 28a-28p. In other examples, a new auxiliary qubit 42 is added to each of the quantum circuits 28a-28p.

[0209] The algorithm program for implementing the quantum power method can be implemented as follows (e.g. Fig.9A As shown in the picture):

[0210] 1. Prepare the initial state: where n can be the number of computational qubits, and na is the number of auxiliary qubits (e.g., log(R)). In this example, no |f> MPS As the initial state, because we prepare the unit vector Might be easier.

[0211] 2. Application of quantum circuits 28, U MPO , which is decomposed into quantum operations 40U according to the isometric subtensors 32 of the extended orthogonal approximation 30 k .

[0212] 3. Regarding adding a new auxiliary space, quantum circuits 28 are added in sequence (the new auxiliary space is initialized as: until the desired power p is reached.

[0213] 4. Measure all auxiliary qubits 44 that are in the computational ground state.

[0214] 5. If the measurement result of the auxiliary qubit 44 corresponds to the bit string 0 ×(p·na) (Successful result), the detection module 18 is used to record the measurement state of the remaining master qubits and restore the bit string corresponding to the candidate optimal value. In the case of an unsuccessful result, return to step 1.

[0215] 6. Repeat steps 1 to 4 to collect measurement statistics and infer a candidate optimal value as the most likely outcome.

[0216] Due to post-selection on the auxiliary space, the success of the algorithm may depend on the post-selection probability, which can be shown to be independent of the number of computational qubits (n) and may scale with the power of the power method according to 1 / p. Taking all steps into account, the running time of the quantum method can be given by t∝n×R 2 ×p. The complexity is only linear in p, not exponential as in the classical (tensor network) implementation of the power method. Therefore, in this case, the quantum method can provide an exponential runtime advantage.

[0217] like Fig. 9B As shown, by executing some quantum operations of quantum circuits 28a-28p in parallel, the quantum circuits used to implement the quantum power method can be efficiently compressed.

[0218] In this example, four quantum circuits 28 are connected in series, each of which includes four quantum operations 40 based on the corresponding isometric subtensors 32 of the orthogonal approximation, U 1 toU 4 However, the second quantum operation 40U of the first quantum circuit 28a 2 Can be used with the first quantum operation 40U of the second quantum circuit 28b 1 Therefore, the running time of the entire quantum operation 28a-28d can be less than the sum of the running times of the individual constituent quantum circuits 28a-28d.

[0219] In some examples, for each power of the quantum power method, the running time of the algorithm on the quantum computer may increase by the running time of one quantum operation 40 corresponding to the last isometric subtensor 32b of the orthogonal approximation 30.

[0220] In some examples, a first output of a first quantum operation 40 of a first quantum circuit 28a in the series sequence of quantum circuits 28a-28p is a second input to a first quantum operation 40 of a second subsequent quantum circuit 28b in the series sequence of quantum circuits 28a-28p, and a second output of the first quantum operation 40 of the first quantum circuit 28 is a first input to a second quantum operation 40 of the second circuit 28b. In some examples, a first output of a second quantum operation 40 of a first quantum circuit 28a in the series sequence of quantum circuits 28a-28p is a second input to a second quantum operation 40 of a second subsequent quantum circuit 28b, and a second output of the first quantum operation 40 of the second quantum circuit 28 is a first input to a second quantum operation 40 of the second circuit 28b, such as Fig. 9B shown.

[0221] FIG. 10A to FIG. 10DExamples of measurement probabilities obtained for a quantum circuit code and a simulated quantum power method for a two-dimensional Ackley function are illustrated in the form of contour plots, with the highest probability located at the center of the corresponding plot.

[0222] Fig. 10A and Fig. 10B The measurement probability of measuring a certain state for the Ackley function implemented with the orthogonal approximation 30 of rank R=8 and R=4 are illustrated respectively. Fig. 10C and Fig. 10D A higher power of an orthogonal approximation 30 of rank R=4 implemented as a quantum circuit 28 on a simulated quantum computer is illustrated, where Fig. 10C and Fig. 10D Examples of power iterations are shown for p=10 and p=25, respectively.

[0223] The 2D Ackley function is defined as:

[0224]

[0225] Consider the domain x,y∈[-5,5] where the function reaches a global minimum f(0,0)=0. Finding this global minimum is a well-known non-convex optimization problem, which poses difficulties for gradient-based optimizers. Therefore, this function is often used to benchmark optimization algorithms. Quantum power iteration is considered here as a method to find the global minimum. To do this, the minimization problem is first transformed into a maximization problem by standard inversion techniques. Specifically, g(x,y)=1–tanh(κf(x,y)).

[0226] Note that for an appropriate choice of κ, minimization of the Ackley function as shown in equation (15) is equivalent to maximization of g(x,y). For simplicity, this modified function g(x,y) is referred to as the Ackley function below. In the example, a total of n = 10 computational qubits (5 qubits per dimension) are used for domain discretization, and the MPS representation is obtained based on the discretized values ​​of the maximum rank r = 6, thereby obtaining the MPS kernel using the cross approximation method in TTPy. From the MPS kernel, the diagonal MPO H is obtained in a straightforward manner by a known algorithm as an approximate MPO representation of the function 22.

[0227] Next, an orthogonal approximation 30U of the approximate MPO representation 22 is obtained by minimizing the cost function as in equation (4) * The algorithm is used to find isometric MPO cores 32 (i.e. isometric subtensors 32) of approximate ranks R=8 and R=4, and it is observed that equation (4) can indeed be minimized by Riemannian gradient descent, where the learning rate is 0.02 and the number of iterations is 10000 in total.

[0228] Although the orthogonal approximation rank R = 8 is 30U * can be projected into the Qiskit quantum circuit, but it requires a compilation step to decompose the 4-qubit operation into a series of two-qubit gates and single-qubit gates. Therefore, the resulting circuit depth becomes relatively large. In order to reduce the simulation time on a classical computer, equation (11) is minimized with an approximate rank R = 4. In this case, convergence is not guaranteed, but the motivation for this approach is still the fact that the power method finds the maximum as specified by equation (14). Although the minimization of equation (4) may produce large errors due to the lack of approximate rank, if the maximum is solved exactly, the orthogonal approximation 30U * Still can be used for power iteration.

[0229] exist Fig. 10A and Fig. 10B In , the orthogonal approximation 30 of the Ackley function with different ranks R = 4 and R = 8 are compared respectively, and the following are observed: 1) when the approximation rank R = 8, the orthogonal approximation 30 qualitatively recovers the exact function; and 2) when the approximation rank R = 4, the orthogonal approximation solves the optimal value of the exact function.

[0230] When R = 4, the orthogonal MPO is projected into the Qiskit quantum circuit and achieves power up to 25 (see Fig. 10C , Fig. 10D ). Fig. 10C , Fig. 10D The results in show that the quantum power method does converge to the optimal value of the Ackley function at (0,0).

[0231] Therefore, a low-rank tensor network approximation 22 of a given function can be prepared, and the low-rank tensor network approximation can be used to obtain an orthogonal approximation 30 through a variational optimization algorithm executed on a classical computer, which can be projected into a quantum circuit 28 consisting of multiple shallow quantum operations with a small number of qubits and using auxiliary qubits 42, 44. Power iteration can be achieved by repeatedly connecting such quantum circuit blocks 28a-28p in series. The linear growth of the quantum power method with the power p shows that quantum algorithms can be efficiently used on recent quantum hardware to solve optimization problems that are difficult to solve on classical hardware.

[0232] In the foregoing discussion, reference has been made primarily to functions for finding extreme values. However, the skilled person will appreciate that the function may be provided in matrix form, and that the method may therefore efficiently find the highest eigenvalue of the matrix as a function. On the other hand, although the examples shown involve equidistant sampling of functions, based on known interpolation and approximation algorithms (such as cross approximation), an MPO / MPS representation of the function may be obtained for any distribution of function values, and the MPO / MPS representation may be used as Figure 2 As part of the method, the approximate MPO representation as a function is presented in 22.

[0233] Further, the skilled person will appreciate that although the illustrated example shows that the qubit measured to determine the result 48 is different from the computational qubit initially prepared to the |+> state, in some examples, the same qubit prepared to the |+> state is also measured to determine the result 48, for example by using additional gates (e.g., a SWAP gate) to swap the quantum states of the qubits.

[0234] The description of the preferred embodiments and the accompanying drawings is only for illustrating the present invention and the beneficial effects associated therewith, and should not be understood as implying any limitation. The scope of the present invention will be determined only by the appended claims.

[0235] List of Reference Numerals

[0236] 10 System

[0237] 12 qubit register

[0238] 14 State Preparation Module

[0239] 16 Quantum Matrix Multiplication Module

[0240] 18 Measurement Module

[0241] 20 Processing System

[0242] 22 MPO approximation

[0243] 24 Matrix Representation of Functions

[0244] 24a Penrose diagram representation of a matrix

[0245] 26 Subtensors

[0246] 28 Quantum Circuits

[0247] 30 Orthogonal approximation

[0248] 32, 32a, 32b isometric subtensors

[0249] 34 Extended Orthogonal Approximation

[0250] 36, 38 |0> state of the boundary core

[0251] 40 Quantum Operations

[0252] 42 auxiliary qubits

[0253] 44 Measurement of auxiliary qubits

[0254] 46 State Initialization Quantum Circuit

[0255] 48 Results

Claims

1. A computer-implemented method for solving an optimization problem for a given function (24, 24a), the method comprising the following steps: obtaining an approximation of the function (24, 24a) in the form of an approximate matrix product operator MPO representation (22) of the function (24, 24a); selecting an approximate rank based on the approximate MPO representation (22); determining a rank-based orthogonal approximation (30, 34) of the approximate MPO representation (22) in the form of an orthogonal MPO of isometric subtensors (32, 32a, 32b) having the approximate rank; encoding the orthogonal approximation (30, 34) into a quantum circuit (28, 28a-p) based on encoding the isometries (32, 32a, 32b) as quantum gates acting on computational qubits and a set of initialization auxiliary qubits (42); and Arbitrary powers of the orthogonal approximation (30, 34) are encoded based on applying the quantum circuit (28, 28a-p) multiple times to an input quantum state in a serial sequence of quantum circuits (28, 28a-p).

2. The method of claim 1, wherein: The approximate MPO representation (22) is an approximate tensor network representation of the function (24, 24a), wherein the method particularly comprises obtaining a series of sub-tensors (26) via a cross approximation of the function (24, 24a).

3. A method as claimed in any one of the preceding claims, wherein: The method further comprises obtaining the approximate MPO representation (22) based on a plurality of sample values ​​of the function (24, 24a) for a corresponding plurality of different input values.

4. A method as claimed in any one of the preceding claims, wherein: Obtaining the approximate MPO representation (22) includes obtaining a matrix product state MPS representation of the function (24, 24a) by a cross approximation technique, and obtaining the approximate MPO representation (22) based on the MPS representation by enforcing diagonal restrictions.

5. A method as claimed in any one of the preceding claims, wherein: Obtaining the orthogonal approximation (30, 34) in the form of the orthogonal MPO comprises: determining an initial guess of the orthogonal approximation (30, 34) using the isometric subtensors (32, 32a, 32b) of the approximate rank; Starting from the initial guess, the orthogonal approximation (30, 34) is iteratively optimized based on an optimization algorithm, which minimizes a cost function under the isometry constraint of the isometric sub-tensors (32, 32a, 32b), wherein the cost function assigns a cost to the orthogonal approximation (30, 34) based on the quality of the orthogonal approximation (30, 34) with respect to the approximate MPO representation (22).

6. A method as claimed in any one of the preceding claims, wherein: Applying the quantum circuit (28, 28a-p) a plurality of times in the form of a series of quantum circuits (28, 28a-p) comprises: preparing a certain number of groups of auxiliary quantum bits (42) to an initial state, wherein the number is at least equal to the number of times the quantum circuit (28, 28a-p) is applied, Each group of auxiliary quantum bits (42) specifically comprises a number of quantum bits equal to the logarithm of the approximate rank.

7. The method of claim 6, wherein: A subsequent quantum circuit (28a-p) in the serial sequence of quantum circuits (28, 28a-p) acts on a corresponding set of auxiliary quantum bits (42) and a set of computational quantum bits that were acted on by the previous quantum circuit.

8. A method as claimed in any one of the preceding claims, wherein Each time the quantum circuit (28, 28a-p) is applied, the method includes measuring a set of post-selected qubits (44) on which the corresponding quantum circuit (28, 28a-p) has acted and on which a subsequent application of the quantum circuit (28, 28a-p) in the series sequence of quantum circuits (28, 28a-p) will not act, Each group of post-selected quantum bits (44) specifically includes a number of quantum bits equal to the logarithm of the approximate rank.

9. A method as claimed in any one of the preceding claims, wherein: Preparing the initial quantum state involves applying a Hadamard gate to each qubit in a set of initial computational qubits.

10. A method as claimed in any one of the preceding claims, wherein: Applying the quantum circuit (28, 28a-p) to an input quantum state p times in the form of a serial sequence of quantum circuits (28, 28a-p) includes preparing p groups of auxiliary quantum bits (42), wherein each group of auxiliary quantum bits (42) is acted upon by a respective one of the quantum circuits (28, 28a-p) and includes a number of auxiliary quantum bits (42) equal to the logarithm of the approximate rank.

11. A processing system (20) for solving an optimization problem for a given function (24, 24a), wherein: The system (10) is configured as follows: obtaining an approximation of the function (24, 24a) in the form of an approximate matrix product operator MPO representation (22) of the function (24, 24a); selecting an approximate rank based on the approximate MPO representation (22); determining a rank-based orthogonal approximation (30, 34) of the approximate MPO representation (22) in the form of an orthogonal MPO of isometric subtensors (32, 32a, 32b) having the approximate rank; An implementation of the orthogonal approximation (30, 34) raised to any power in a quantum circuit (28, 28a-p) is determined based on encoding the isometries (32, 32a, 32b) as quantum gates acting on computational qubits and a set of initialization auxiliary qubits (42), wherein in the implementation, the quantum circuit (28, 28a-p) is applied multiple times to an input quantum state in the form of a series of quantum circuits (28, 28a-p).

12. The system (10) of claim 11, wherein: In order to obtain the orthogonal approximation (30, 34) in the form of the orthogonal MPO, the system (10) is configured to: determining an initial guess of an orthogonal approximation (30, 34) of said function (24, 24a) in the form of a tensor network of isometric sub-tensors (32, 32a, 32b) of said approximation rank; Starting from the initial guess, the orthogonal approximation (30, 34) is iteratively optimized based on an optimization algorithm, which minimizes a cost function under the isometry constraint of the isometric sub-tensors (32, 32a, 32b), wherein the cost function assigns a cost to the orthogonal approximation (30, 34) based on the quality of the orthogonal approximation (30, 34) with respect to the approximate MPO representation (22).

13. The system (10) of claim 11 or 12, wherein: The system (10) is further configured to transmit the implementation to a quantum computing system to execute the implementation on quantum hardware.

14. A hybrid quantum-classical computing system, comprising a system (10) as claimed in any one of claims 11 to 13 and quantum computing hardware, wherein: The hybrid quantum-classical computing system is configured as: receiving the implementation from the processing system (20), implementing said embodiments in said quantum computing hardware, and A computation result is received from the quantum computing hardware.

15. A non-transitory medium comprising machine-readable instructions which, when executed by a processing system (20), implement the method according to any one of claims 1 to 10 and / or the system according to any one of claims 11 to 14.

Images