Method and system for black-box optimization using quantum computation

Abstract

To provide a method and system for solving optimization problems using quantum computation.SOLUTION: A computer-implemented method for solving an optimizing problem given a function comprises the steps of: obtaining an approximation of the function in the form of an approximated matrix product operator (MPO) representation of the function, selecting an approximation rank based on the approximated MPO representation; determining a rank-based orthogonal approximation of the approximated MPO representation in the form of an orthogonal MPO with isometric sub-tensors of the approximation rank; encoding the orthogonal approximation into a quantum circuit based on encoding of the isometric sub-tensors into quantum gates; and encoding an arbitrary power of the orthogonal approximation based on applying the quantum circuit in a concatenated sequence of quantum circuits multiple times to an input quantum state.SELECTED DRAWING: Figure 2

Description

[Technical field]

[0001] The present invention relates to the field of encoding computational operations into quantum circuits. More precisely, the invention relates to the optimization of the computer implementation of computational operations aimed at encoding in quantum circuits. [Background technology]

[0002] Quantum computers provide a platform for controllable quantum mechanical systems whose states and interactions can be controlled to perform computations. Computation is realized through the deterministic evolution of controllable quantum mechanical systems, e.g., qubits, which are the quantum analogues of classical bits, whose states can be measured to determine the outcome of the computation.

[0003] Control operations on these qubits are called quantum gates. Quantum gates can act coherently on qubits, for example to entangle the states of multiple qubits and any combination of them, both to cause a change in the state of a single qubit (so-called single-qubit gates) and to act on multiple qubits (so-called multi-qubit gates). For example, a single-qubit gate may cause a rotation of the spin state of an electron by a selectable value, for example π / 2. A multi-qubit gate may act coherently on two or more qubits, such as a coherent CNOT operation on the states of two qubits. Multiple quantum gates can be applied to qubits in a quantum computer in parallel or in sequence to perform a computation. Finally, the state of the qubits may be repeatedly measured after applying a sequence of quantum gates to determine the probability for each possible outcome of the computation.

[0004] In essence, the operation of a quantum computer can be thought of as relying on encoding a starting state into an internal quantum state of a qubit, followed by multiplying the internal quantum state of the qubit by a matrix that can be implemented by a combination of different quantum gates, and measuring the resulting result.

[0005] To compute solutions to problems that are considered intractable for classical computers, quantum computers can exploit the special properties of quantum mechanical states, in particular the superposition and entanglement of different quantum states, to find the solution in a relatively small number of computational steps. Moreover, because quantum operations can create complex superposition states of qubits, quantum computers may, in principle, have access to large internal memories for processing computational states.

[0006] However, as found by [1], converting any operation into a combination of single and two-qubit gates generally requires a number of operations that is exponential with respect to the number of qubits.

[0007] Non-Patent Document 2 teaches solving nonlinear differential equations by variational quantum computing. The outputs of quantum variation circuits are combined by a quantum nonlinear processing unit, and a Hamiltonian expressible as an MPS is compared with an MPS ansatz encoded as a quantum circuit according to the MPS quantum circuit encoding method. [Prior art documents] [Non-patent literature]

[0008] [Non-Patent Document 1] Shende et al., Minimal universal two-qubit controlled-not-based Circuits, Phys. Rev. A, 69:062321, June 2004 [Non-Patent Document 2] Lubasch et al., Variational quantum algorithms for nonlinear problems Summary of the Invention

[0009] However, current approaches are still limited by the significant complexity involved in performing the steps of encoding the state into qubits, encoding the matrix through quantum gates, and measuring. In particular, the exponential complexity of decomposing any multi-qubit gate with respect to the number of qubits, combined with the need for full connectivity between the qubits, presents a significant obstacle. Although algorithms exist for converting large matrices into combinations of single quantum gates and two-qubit operations (e.g., CNOT gates), the complexity of the resulting circuits generally scales exponentially with the number of qubits. The use of a Hamiltonian expressible as an MPS constrains the applicable problem instances for the algorithm.

[0010] In particular, with any tensor network representation, it is not easy to generate high powers of an intended operation. However, some computational schemes can take advantage of access to high power computational operations, such as power methods to determine 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, thereby enabling the solution of a large class of optimization processes using quantum processing resources.

[0011] In view of this state of the art, it is an object of the present invention to provide a computer-implemented method, and corresponding system, for determining a quantum encoding of an approximation of a function, in the form of an encoded matrix and any powers thereof, that can closely approximate the black-box function and any powers thereof, but that can also be implemented on quantum hardware without the depth of the quantum circuit scaling exponentially with the number of qubits.

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

[0013] According to a first aspect, the present invention relates to a computer-implemented method for solving optimization problems of a function. The method includes obtaining an approximation of the function in the form of an approximated matrix multiplication operator (MPO) representation of the function, and selecting an approximation rank based on the approximated MPO representation. The method further includes determining a rank-based orthogonal approximation of the approximated MPO representation in the form of an orthogonal MPO having an isometric subtensor of the approximate rank. The method further includes encoding the orthogonal approximation into a quantum circuit based on encoding the isometric subtensor into quantum gates operating on an initialized set of computation qubits and ancilla qubits, and encoding any power of the orthogonal approximation based on applying the quantum circuit multiple times to an input quantum state in a concatenated sequence of the quantum circuit.

[0014] The function may be a black-box function, may be provided in terms of a set of inputs and correlated outputs, or may be based on a discretization of a smooth function. The function may generally be mapped to a matrix, for example, 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 a representation thereof.

[0015] Furthermore, given a set of function values ​​for corresponding input values, there are several algorithms for generating approximated representations that implement the function, such as the TT cross approximation, which can generate tensor network approximations, typically in the form of MPS or MPO.

[0016] In a preferred embodiment, the approximated MPO representation is an approximate tensor network representation of the function, and the method includes, inter alia, the step of obtaining a set of sub-tensors via cross approximation of the function.

[0017] Tensor network approximations may interpolate and / or approximate a function based on a set of pairs of inputs and corresponding outputs of the function, such as by sampling the function with a set of randomly selected inputs. For example, a tensor train (TT) approximation of a function may be obtained using TT cross approximation. Such tensor network representations can, in principle, be stored efficiently in memory despite the large dimensions of the underlying matrices corresponding to the discretized functions. Preferably, the tensor network approximation has a relatively small rank.

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

[0019] For example, discretized sample values ​​of a function for the relevant inputs, e.g., a multivariate black-box function, may be used to obtain a tensor network and / or MPS approximation of the function using skeleton decomposition or cross approximation.

[0020] In a preferred embodiment, obtaining an approximated MPO representation includes obtaining a matrix product state (MPS) representation of the function by a cross approximation technique, and obtaining the approximated MPO representation based on the MPS representation by enforcing diagonal restrictions.

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

[0022] The approximated MPO representation of a function can have a relatively small bond dimension. However, the approximated MPO representation / approximation does not directly correspond to a quantum circuit that can be determined without the difficulty of exponential computation. Rather, an approximation rank is selected that should be greater than or equal to the rank of the approximated MPO representation to approximate the function with an orthogonal approximation for which efficient encoding strategies exist, based on the approximated matrix multiplication operator representation.

[0023] In a preferred embodiment, the approximation rank is greater than the rank of a matrix multiply operator representation of the function, which rank is in particular the bond dimension of the matrix multiply operator representation.

[0024] Isometry constraints on isometric subtensors may benefit from increased rank due to advantages associated with the increased degrees of freedom, such as maintaining sufficient approximation accuracy for general problems.

[0025] It should be noted that for some matrices / functions, an exact matrix multiply operator representation may be determined, but in general, the matrix multiply operator representation may be obtained by an approximation algorithm, such as cross approximation, or an algorithm based on singular value decomposition, to approximate functions in MPO without significantly affecting the results of the method. Those skilled in the art will further appreciate that matrix multiply operator representations / approximations are already known and it may be sufficient to receive the matrix multiply operator representation from a suitable source.

[0026] The approximated MPO representation may include non-isometric subtensors, e.g., less than 1% or 10, depending on the application. -6 The method approximates a function with an error that may be less than 1%, and may also be referred to as a non-orthogonal MPO approximation of a function. For example, an approximation of a function with an error of about 1% may be sufficient to solve an optimization problem. Although not every matrix can be represented by an MPO representation with an arbitrarily small accuracy, those skilled in the art will understand that they may further choose to apply the method taking into account the increased inaccuracy of the MPO approximation, for example, when other ways of performing the calculation are not possible or when a relatively low accuracy is sufficient for the respective calculation problem.

[0027] Based on the bond dimension (rank) of the approximated MPO representation, the approximate rank is the next power of two, i.e., 2 N (N is an integer), or may be chosen as any power of two greater than the bond dimension of the approximated MPO representation. As an example, if the rank of the approximated MPO representation of a function is 3, then the rank chosen for the isometric subtensors of the orthogonal approximation might be 4 or 8. However, if the approximation rank is less than or equal to the rank of the approximated MPO representation of a function, then although in principle a good approximation may be made, the orthogonal approximation is likely to be an inaccurate representation of the function.

[0028] Based on an implementation of an orthogonal approximation of an approximated MPO representation of a function, the maximum eigenvalue of the function / approximated MPO representation may be obtained by successively applying the quantum circuit encoding of the orthogonal approximation to the input states multiple times and measuring the output, in a manner similar to the classical "power method".

[0029] To determine the orthogonal MPO approximation, the method includes determining an initial guess for the orthogonal MPO, which corresponds to an MPO having an isometric subtensor with appropriate rank. The initial guess may be constructed from random subtensors that satisfy the isometry condition.

[0030] In a preferred embodiment, obtaining an orthogonal approximation in the form of an orthogonal MPO comprises determining an initial guess for the orthogonal approximation having an isometric subtensor of approximation rank, and, starting from the initial guess, iteratively optimizing the orthogonal approximation based on an optimization algorithm that minimizes a cost function subject to an isometry constraint of the isometric subtensor, where the cost function attributes a cost to the orthogonal approximation based on the quality of the orthogonal approximation with respect to the approximated MPO representation.

[0031] In some embodiments, the orthogonal approximation is a matrix multiplication operator, and the orthogonal approximation is mathematically equivalent to a tensor network, where each intermediate isometric subtensor has two exterior uncontracted indices, as well as an interior index that is chain-contracted with adjacent isometric subtensors.

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

[0033] In some embodiments, the approximate rank is the common bond dimension of all the subtensors of the tensor network.

[0034] The initial estimate is then optimized based on an iterative optimization algorithm such that the difference between the orthogonal approximation and the function approximated by the approximated MPO representation is minimized, which may be based on a norm of the difference, for example the Frobenius norm or the cosine norm of the difference between the orthogonal approximation and the approximated MPO representation of the function, in particular its subtensors.

[0035] In a preferred embodiment, the cost function is based on a difference function of the function and a renormalized orthogonal approximation evaluated based on a tensor network calculation, where the renormalized orthogonal approximation is based on the orthogonal approximation multiplied by a renormalization constant, the cost function is in particular based on the Frobenius norm of the difference function evaluated based on the tensor network calculation, and / or the renormalization constant is selected such that the approximated MPO representation of the function divided by the normalization constant does not increase the trace of the states on which the function operates and / or such that the trace of the density matrix on which the renormalized orthogonal approximation operates is less than or equal to 1 regardless of the initial state.

[0036] The renormalization constant allows preserving the constraint that quantum operations cannot increase the trace of the density matrix of the quantum state, i.e., given an arbitrary input density matrix ρ, the orthogonal approximation A can be expressed as

number

[0037] Given an arbitrary matrix M representing a function, we can find an orthogonal approximation A with a differential cost function C=||cA-M|| 2 where || || represents a chosen norm of the difference, e.g., the Frobenius norm. In other words, the operation to be implemented in a quantum computer may be renormalized by a renormalization matrix c -1 It may be M.

[0038] In a preferred embodiment, a renormalization constant c is set to minimize the cost function:

number

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

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

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

[0042] In other words, the method may include a step of (probabilistically) determining the gradient of the difference between the orthogonal approximation and an approximated MPO representation of the function, and subsequently updating the orthogonal approximation according to the determined gradient.

[0043] Those skilled in the art will appreciate that in some embodiments, minimizing the difference between the orthogonal approximation and the approximated matrix multiplication operator representation of the function may be preferred as this may reduce the computational complexity of the optimization. For example, the difference between the subtensors of the orthogonal approximation and the approximated matrix multiplication operator representation of the function may be determined and the respective isometric subtensors updated based on the difference. The update of the isometric subtensors should preserve the isometry constraint, i.e., the updated isometric subtensors must be isometric.

[0044] In a preferred embodiment, the orthogonal approximation is a matrix multiplication operator, which can be expressed as:

[0045]

number

[0046] Here, V k is isometric, i.e., for any index k of the isometric subtensor

number

[0047] The alternating descent method may be performed as a stochastic gradient descent method on each of the subtensors, e.g., in a predetermined order or random order, such that each of the isometric subtensors of the orthogonal approximation approaches a subtensor of the approximated matrix product operator representation of the function.

[0048] In a preferred embodiment, the gradient of the cost function is determined based on multiple partial derivatives of the subtensors of the orthogonal approximation, e.g., by determining partial derivatives of the cost function with respect to the values ​​of each subtensor, and by updating the elements of the orthogonal approximation based on multiple subtensors, e.g., the partial derivatives of all the subtensors. Those skilled in the art will appreciate that the optimization includes stochastic optimization, e.g., it is not necessary to determine partial derivatives for all parameters of the orthogonal approximation at each iteration, but may update a stochastically determined subset of the parameters, such as in stochastic gradient descent and related optimization algorithms.

[0049] However, the updated subtensors should preserve the isometry constraints, for example, by updating the isometric subtensors according to a Riemann gradient descent method, which may preserve the isometry constraints, for example, by projecting the update vectors and / or the updated subtensors onto the Stiefel manifold of the isometric tensors.

[0050] For example, optimizing an orthogonal approximation under an isometry constraint is a projection of the partial derivatives of the cost function.

number

number

[0051] For example, subtensor V k is an optimization algorithm similar to steepest descent optimization,

number

[0052] In some embodiments, optimizing the orthogonal approximation under the isometry constraint involves approximating the sub-tensor V based on partial derivatives of the cost function and, optionally, information from previous iteration steps of the optimization algorithm. k Determine the search direction pk of point V k Stiefel manifold of m×p isometric matrices in

number

number

[0053] 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.

[0054] In some embodiments, optimizing the orthogonal approximation under the isometry constraint includes projecting the modified subtensor and / or the modified orthogonal approximation, where the modified subtensor and / or the modified orthogonal approximation are modified on a Stiefel manifold for the subtensor and / or the orthogonal approximation, respectively, based on a search direction and / or an update direction, and / or projecting the subtensor Vk Retraction of updates to the manifold

number

number

[0055] Retraction, for example,

number

number

[0056] Those skilled in the art will appreciate that this technique can be extended to optimization techniques that may include, for example, momentum based on gradients determined in previous iterations, thereby providing an update vector

number

[0057] In a preferred embodiment, optimizing the orthogonal approximation under the isometry constraint comprises updating the orthogonal approximation by updating all the isometric subtensors, where for each isometric subtensor V k is the update direction point V k Stiefel manifold of m×p isometric matrices in

number

[0058] Those skilled in the art will appreciate that different retractions are used in embodiments such as Cayley retraction or singular value decomposition, to name just two examples, and that those skilled in the art may select an appropriate retraction from multiple known retractions.

[0059] The isometric subtensor can be reshaped into a matrix that is square with respect to the intermediate subtensor of the orthogonal approximation, and the reshaped matrix corresponding to the intermediate subtensor is unitary based on the isometry constraint of the subtensor during the iterative optimization of the orthogonal approximation.

[0060] In some embodiments, the isometric sub-tensors of the orthogonal approximation at the ends of the orthogonal approximation are isometric and the intermediate tensors are unitary when reshaped into two-dimensional matrices before extending the orthogonal approximation, such that all isometric sub-tensors are unitary after reshaping the isometric sub-tensors into square matrices.

[0061] To obtain a set of unitary matrices that can be implemented in a quantum circuit, a sub-tensor at the end of the chain of tensors (sometimes called the boundary core) can be expanded using an ancilla qubit, for example by preparing a zero-state |0> ancilla qubit that is multiplied with the first isometric sub-tensor, and by implementing a post-selection at the opposite end of the chain of tensors.

[0062] Thus, the method may include expanding the boundary core of the orthogonal approximation such that all isometric subtensors become unitary after reshaping the isometric subtensors into square matrices. Expanding the orthogonal approximation may be part of determining an initial guess for the orthogonal approximation or may be performed at the beginning or later when the orthogonal approximation is iteratively optimized. For example, the initial guess may be constructed from subtensors that can be reshaped into unitary matrices and the determination of the cost function may take into account the implicit expansion of the orthogonal approximation.

[0063] In some embodiments, the step of extending the orthogonal approximation is based on providing one or more additional ancilla qubits as input to the first end of the tensor network, and the method specifically provides at least logR additional ancilla qubits, where R is the approximation rank.

[0064] In some embodiments, the step of extending the orthogonal approximation includes measuring a qubit at a second end of the tensor network, and selecting a result having a predetermined measurement result of the qubit measured at the second end of the tensor network is part of the implementation of the approximated MPO representation of the function.

[0065] Based on the addition of ancilla qubits and post-selection, all subtensors may be square after appropriate reshaping, and thus form a set of unitary matrices that can be implemented in quantum circuits by decomposing the isometric subtensors into combinations of quantum gates of the respective quantum architecture, e.g., single-qubit rotations and CNOT gates.

[0066] Those skilled in the art will understand that the expansion can be performed before or after iteratively optimizing the orthogonal approximation. In other words, it is possible to define an initial guess based on an expanded matrix multiplication operator, or to expand the orthogonal approximation after it has been iteratively optimized, for example by introducing and filling additional elements of the boundary core, such that the resulting boundary core becomes unitary after appropriate reshaping.

[0067] In some embodiments, the method further comprises implementing an orthogonal approximation of the function as a quantum circuit on quantum hardware based on the encoding of the isometric subtensors into quantum gates.

[0068] The decompositions of each isometric subtensor are combined to perform an orthogonal approximation based on sequential application of each isometric subtensor, and the output state of the logR qubits resulting from a unitary operation corresponding to one of the subtensors may be the input for the unitary operation corresponding to the next isometric subtensor in the tensor network.

[0069] In some embodiments, implementing an orthogonal approximation of the function as a quantum circuit is based on sequentially applying isometric sub-operations to a set of qubits and rearranging unitary matrices in a quantum network.

[0070] As a result, the depth of a quantum circuit can only scale polynomially with the number of qubits, in contrast to the exponential scaling of the typical decomposition of matrices into single- and two-qubit gates.

[0071] In a preferred embodiment, the step of applying a quantum circuit multiple times in a concatenated sequence of quantum circuits includes preparing a number of sets of ancilla qubits in an initial state, the number being at least equal to the number of times the quantum circuit is to be applied.

[0072] Specifically, for each extended orthogonal approximation, a corresponding input of an ancilla qubit is prepared, which is initialized to a respective initial state of the ancilla qubit, e.g., "0", "1", or "ground" state, respectively.

[0073] In a preferred embodiment, each set of ancilla qubits contains a number of qubits equal to the logarithm of the approximate rank.

[0074] In a preferred embodiment, a subsequent quantum circuit in a concatenated sequence of quantum circuits operates on each of the sets of ancilla qubits and computation qubits operated on by the previous quantum circuit.

[0075] In other words, the method may begin with initializing a set of initial computation qubits and a first set of ancilla qubits to be acted on by a first quantum circuit in a concatenated sequence of quantum circuits, and a second quantum circuit in the concatenated sequence of quantum circuits may act on the first set of ancilla qubits, a second set of ancilla qubits, and a subset of the initial computation qubits, which may be smaller than the set of initial computation qubits by the number of ancilla qubits in the second set of ancilla qubits.

[0076] In a preferred embodiment, each time a quantum circuit is applied, the method includes measuring a set of post-selected qubits that the respective quantum circuit acts on and that are not acted on by subsequent applications of the quantum circuit in a concatenated sequence of quantum circuits.

[0077] The post-selection may be based on the results of measurements of the post-selected qubits of each of the quantum circuits in the concatenated sequence of quantum circuits, for example, the results may be recorded for purposes of determining the largest eigenvalue of the approximated MPO representation of the function when all the post-selected qubits are measured in the |0> state.

[0078] In a preferred embodiment, each set of post-selected qubits contains a number of qubits equal to the logarithm of the approximate rank.

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

[0080] Applying a Hadamard gate to each qubit of the initial computation qubits may prepare a vector of computational basis states of "1" to support a quantum exponentiation method for determining the maximum eigenvalue of the function.

[0081] A quantum computing system implementing the method may include a number of initial computation qubits equal to the number of isometric subtensors, and in particular, each quantum operation implementing one of the isometric subtensors of a first quantum circuit in a concatenated sequence may operate on a respective one of the initial computation qubits. A first quantum operation in a cascade of quantum operations implementing a quantum circuit may operate on a set of ancillary qubits. Subsequent quantum operations in the cascade may operate on a respective one of the computation qubits and a number of processing qubits equal to the logarithm of the approximate rank, the processing qubits having been operated on by a previous quantum operation in the cascade of quantum operations.

[0082] Each quantum operation may have as an output the quantum state of one of the qubits that is not acted upon by a subsequent quantum operation in the cascade of quantum operations. The output of the quantum operation may form the computational qubit of the next quantum circuit in the concatenated sequence of quantum circuits. The final quantum operation in the cascade of quantum operations may include as an output a quantum state for processing by the subsequent quantum circuit, as well as a number of processing qubits equal to the logarithm of the approximate rank that may be measured as a post-selection qubit.

[0083] 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 be more than a single quantum bit, and a corresponding number of initial computational quantum bits may be provided for each quantum operation of the first circuit, however, in the following we will primarily refer to examples having a single initial computational quantum bit per quantum operation.

[0084] In a preferred embodiment, the step of applying the quantum circuits in the concatenated sequence of quantum circuits p times to the input quantum state includes preparing p sets of ancilla qubits, each set of ancilla qubits acted upon by a respective one of the quantum circuits and including a number of ancilla qubits equal to the logarithm of the approximate rank.

[0085] Each quantum circuit after the first quantum circuit in the concatenated sequence of quantum circuits may operate on a respective one of the set of ancilla qubits, and on a processed computation qubit operated on by the previous quantum circuit in the concatenated sequence, where the processed computation qubit does not form part of the post-selected qubit of the previous quantum circuit.

[0086] According to a second aspect, the 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 approximated matrix multiplication operator MPO representation of the function and to select an approximation rank based on the approximated MPO representation. The system is further configured to determine a rank-based orthogonal approximation of the approximated MPO representation in the form of an orthogonal MPO with an isometric subtensor of the approximate rank. The system is further configured to determine an implementation of an arbitrary power of the orthogonal approximation in a quantum circuit based on an encoding of the isometric subtensor into quantum gates acting on the computation qubits and an initialized set of ancilla qubits, in an implementation where the quantum circuit is applied multiple times to an input quantum state in a concatenation sequence of the quantum circuit.

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

[0088] The processing system may include a single processing unit or 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 processing devices that operate 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 interfaces for connecting to sensors, devices, instruments, integrated logic circuits, other controllers, etc., and the interfaces may be configured to receive or transmit signals, such as electrical signals, optical signals, radio signals, acoustic signals, etc.

[0089] The system may be configured to obtain the approximated MPO representation, for example by performing a cross approximation or decomposition algorithm of a function, or by receiving the approximated MPO representation from a suitable source, such as a database.

[0090] In a preferred embodiment, to obtain an orthogonal approximation in the form of an orthogonal MPO, the system is configured to determine an initial estimate of the orthogonal approximation of the function, in the form of a tensor network with an isometric subtensor of approximation rank, and, starting from the initial estimate, iteratively optimize the orthogonal approximation based on an optimization algorithm that minimizes a cost function subject to the isometry constraint of the isometric subtensor, where the cost function attributes a cost to the orthogonal approximation based on the quality of the orthogonal approximation with respect to the approximated MPO representation.

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

[0092] In a preferred embodiment, the system is further configured to communicate the implementation to a quantum computing system for executing the implementation on quantum hardware.

[0093] An embodiment may allow for determining the maximum element of the orthogonal approximation of the function, in particular the maximum eigenvalue.

[0094] According to a third aspect, the present invention relates to a hybrid quantum-classical computing system comprising a system according to the second aspect and quantum computing hardware, the hybrid quantum-classical computing system being configured to receive an implementation from a processing system, to implement the implementation in the quantum computing hardware and to receive a computation result from the quantum computing hardware.

[0095] According to a fourth aspect, the present invention relates to a non-transitory medium comprising machine-readable instructions which, when executed by a processing system, implements the method according to the first aspect and / or the system according to the second aspect and / or the third aspect.

[0096] The features and many advantages of the method and system according to the present invention will be best understood from the following detailed description of the preferred embodiment, when read in conjunction with the accompanying drawings. [Brief description of the drawings]

[0097] [Figure 1] FIG. 1 is a schematic diagram of an example quantum computing system. [Diagram 2] FIG. 1 illustrates an example of a computer-implemented method for encoding higher powers of function approximations in quantum circuits. [Diagram 3] FIG. 13 illustrates generally a flow chart for determining a matrix multiplication operator representation based on a given target matrix. [Figure 4] FIG. 13 illustrates a schematic diagram of a strategy for determining a quantum circuit encoding for an approximation of an objective matrix, according to an example. [Diagram 5] FIG. 5 illustrates results obtained from applying the encoding strategy described in relation to FIG. 4 to two different target matrices. [Figure 6] FIG. 5 illustrates results obtained from applying the encoding strategy described in relation to FIG. 4 to two different target matrices. [Figure 7A] FIG. 13 illustrates experimental run-time results by implementing a method for determining a quantum circuit when approximating a diagonal matrix, according to one example. [Figure 7B] FIG. 13 illustrates experimental run-time results by implementing a method for determining a quantum circuit when approximating a diagonal matrix, according to one example. [Figure 8] FIG. 1 is a schematic flow chart diagram for determining a quantum computer implementation of an algorithm for finding the largest eigenvalue of a matrix that approximates a function, according to an example. [Figure 9A] 1A-1C show detailed implementations of a quantum exponentiation method for finding the largest eigenvalue of a matrix that approximates a function based on quantum circuit encoding of computational operations, according to respective examples. [Figure 9B]1A-1C show detailed implementations of a quantum exponentiation method for finding the largest eigenvalue of a matrix that approximates a function based on quantum circuit encoding of computational operations, according to respective examples. [Figure 10A] FIG. 13 shows, as contour plots, examples of measurement probabilities obtained for quantum circuit encoding and simulated quantum exponentiation of a two-dimensional Ackley function. [Figure 10B] FIG. 13 shows, as contour plots, examples of measurement probabilities obtained for quantum circuit encoding and simulated quantum exponentiation of a two-dimensional Ackley function. [Figure 10C] FIG. 13 shows, as contour plots, examples of measurement probabilities obtained for quantum circuit encoding and simulated quantum exponentiation of a two-dimensional Ackley function. [Figure 10D] FIG. 13 shows, as contour plots, examples of measurement probabilities obtained for quantum circuit encoding and simulated quantum exponentiation of a two-dimensional Ackley function. DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS

[0098] 1 illustrates a schematic diagram of an example quantum computing system 10. System 10 includes a qubit register 12 with computation qubits, a state preparation module 14 for preparing the quantum states of the computation qubits in 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 results of the quantum computation, e.g., by projective measurement of the states of some or all of the qubits. A processing system 20 can initialize the computation qubits in qubit register 12 and define the operation of state preparation module 14 and / or quantum matrix multiplication module 16, and may retrieve the measurement results from measurement module 18, e.g., by determining and sending control actions to hardware of quantum computing system 10.

[0099] State preparation module 14 and quantum matrix multiplication module 16 may be similar in configuration and may implement a combination of single and multi-qubit quantum gates that affect the state of the qubits in qubit register 12. Note that in some embodiments, state preparation module 14 may not be required and the combination of quantum gates determined by quantum matrix multiplication module 16 may directly act on the quantum states of the qubits initialized in qubit register 12 at some initialized quantum states. The types of quantum gates implemented in the quantum circuit by state preparation module 14 and / or quantum matrix multiplication module 16 may depend on the architecture of the particular quantum hardware utilized, for example, based on the native gates of the quantum hardware. State preparation module 14 and quantum matrix multiplication module 16 may calculate quantum gate configurations and / or send control operations to control the quantum hardware so that the quantum hardware realizes the desired computational operations for the particular quantum states.

[0100] However, decomposing any computational operation into quantum circuits is generally a hard computational problem whose complexity usually scales exponentially with the number of qubits involved. Furthermore, for common operations, the depth of the circuit, i.e., the number of successive quantum gates to realize a quantum operation, may scale exponentially with the number of qubits. While efficient encoding of the problem into combinations of quantum gates is known for certain types of problems, the encoding of many desirable problem types into quantum circuits may still be a significant obstacle to realizing quantum computing.

[0101] 2 illustrates an example of a computer-implemented method for encoding an arbitrary power of a matrix multiplication operator approximation of a function in a quantum circuit. The method includes obtaining an approximated matrix multiplication operator representation of the function (S10) and selecting an approximation rank based on the MPO representation (S12). The method further includes determining a rank-based orthogonal approximation of the approximated MPO representation in the form of an orthogonal MPO having an isometric subtensor of the approximate rank. The method further includes encoding the orthogonal approximation into a quantum circuit based on encoding the isometric subtensor into a quantum gate (S16), and encoding the arbitrary power of the orthogonal approximation based on applying the quantum circuit multiple times in a concatenation sequence of the quantum circuit to an input quantum state (S18).

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

[0103] The cost function may attribute a cost to the orthogonal approximation of the function based on the quality of the orthogonal approximation with respect to the approximated MPO representation. Thus, instead of directly implementing the approximated MPO representation in a quantum circuit, the method includes approximating the approximated MPO representation with an orthogonal approximation, constructed in the form of a tensor network with an isometric subtensor (also referred to below as an isometric core).

[0104] The isometry constraint of the subtensors of the orthogonal approximation can be enforced at each step of the iterative optimization, for example based on modified gradient descent and / or retraction of any updated subtensors into isometric subtensors, so that the optimization algorithm converges towards an approximation of the function with the isometric subtensor. In general, the approximation rank should be at least equal to the rank of the approximated matrix multiplication operator representation of the function, so that the isometric subtensor constraint can be at least partially compensated for by the large degrees of freedom of the contracted index of the orthogonal approximation.

[0105] The isometric subtensors are reshaped into square matrices for any subtensors between the ends of the tensor network (so-called boundary cores), which may be unitary after appropriate reshaping. These subtensors that may be reshaped into unitary matrices can be directly implemented as quantum circuits with known decomposition algorithms, which may be computed efficiently since the size of the reshaped matrices / isometric subtensors is small. Those skilled in the art will understand that the subtensors may be implemented without explicitly reshaping the isometric subtensors into matrices.

[0106] To implement the boundary core, i.e., the isometric subtensors placed at the ends of the tensor network in the form of a tensor string, the ends of the tensor network representing the orthogonal approximation may be expanded such that all isometric subtensors can be reshaped into unitary matrices, for example, by providing an ancilla qubit at one end of the tensor network and introducing a post-selection measurement at the opposite end, thereby reshaping the boundary core equally to a square matrix. Such an expansion may be pre-configured as part of the initial guess of the orthogonal approximation, or may be configured as the orthogonal approximation is iteratively optimized, for example, by completing the boundary core with additional indices such that it is a square matrix, and filling in missing matrix elements such that when the appropriate reshaping is applied, the resulting matrix is ​​unitary.

[0107] As a result, orthogonal approximations can be implemented in quantum circuits based on decompositions of isometric subtensors that are smaller in memory than a matrix representation of the full function. As a result, the circuit depth of a quantum circuit for an entire function can only scale polynomially with the number of qubits, which determines the size of the memory space to hold the computational state, making complex computations possible on near-term quantum hardware.

[0108] FIG. 3 shows an example of a square matrix M NxN FIG. 1 shows a flow chart for determining a matrix multiplication operator representation 22 based on a given function provided in the form of a matrix 24 (denoted as x, y, z) and a Penrose graphical representation 24a. The matrix 24 may correspond to discretized function values ​​of a function implemented on quantum hardware and has computational basis states for a qubit register 12 having n qubits, where N=2 n It is.

[0109] A matrix 24 associated with a function may be represented in a Penrose diagram representation 24a by a circle with two uncontracted indices of dimension N. This function 24 can be approximated or represented by a matrix multiplication operator representation 22 composed of multiple sub-tensors 26, where a summation by the contracted indices recovers each element of the matrix 24 representing the original function, for example according to:

[0110]

number

[0111] matrix 24, A l m The MPO representation of 22 is the uncontracted index l k , m k , and the contracted index j k , j k-1 26, V kEach subtensor 26, shown as a circle in Figure 3, has four indices, except for the boundary cores at either end of the MPO representation 22, which have three, and the contracted indices (projecting horizontally) are the MPO bond indices of bond dimension r. The uncontracted indices have the dimension of the quantum state by one qubit, in this example two.

[0112] Such an MPO representation 22 may use a number of known algorithms to find the matrix 24 or an approximation thereof. In principle, any matrix 24 may be represented in a matrix multiplication operator (MPO) representation 22 by decomposing the matrix 24 into a set of matrix products that have the property that their matrix products are equal to certain elements of the matrix 24.

[0113] When only a small number of function values ​​are known for a black-box function, a matrix multiplication operator / state approximation of the function may be obtained by interpolation, for example using cross-approximation techniques such as skeleton decomposition or TT cross-approximation to obtain an approximated tensor network of the function based on the known function values. The approximated tensor network may correspond to an MPO or MPS approximation of the function.

[0114] However, although it is possible to efficiently encode MPS states in quantum circuits, as shown, for example, by Lubasch et al., no general MPO-compatible strategy exists.

[0115] Nevertheless, the method illustrated in the example of FIG. 2 may be used to provide a quantum circuit for the approximation of a function.

[0116] 4 shows a schematic of a strategy for determining which quantum circuits 28 to encode for the approximation of a function according to one example. The function may first be represented by an approximated MPO representation 22 having log(N) (i.e. the logarithm of the order of the square matrix to be represented), or by a part of it, a sub-tensor 26, having rank r, which can be obtained based on suitable decomposition and / or approximation algorithms.

[0117] Based on the rank r of the approximated MPO representation 22, an approximation rank R is selected, which should be greater than or equal to the rank r of the approximated MPO representation 22 and should be a power of two, i.e., 2 k where k is an integer. Based on the selected approximation rank R, the method may generate an initial estimate of the orthogonal approximation 30 having isometric subtensors 32, 32a, 32b such that all isometric subtensors 32, 32a, 32b are isometric when reshaped into a two-dimensional matrix and the bond dimension of the orthogonal approximation 30 is equal to the approximation rank R. The isometric subtensors 32, 32a, 32b for the initial estimate of the orthogonal approximation 30 may be generated randomly, e.g., selected from random isometric subtensors 32, 32a, 32b.

[0118] Starting from an initial estimate of the orthogonal approximation 30, the method may proceed by iteratively optimizing the orthogonal approximation 30 based on an optimization algorithm with the objective of minimizing a distance metric to the approximated MPO representation 22, wherein intermediate steps of the optimization algorithm may maintain isometry of the sub-tensors 32, 32a, 32b.

[0119] The approximated MPO representation 22, when reshaped into a matrix, may not be unitary and therefore cannot be directly implemented on a quantum computer. Furthermore, the quantum operation A must preserve the trace, which is given by

[0120]

number

[0121] For example, it is preferable to minimize a difference metric for the renormalized orthogonal approximation 30 according to:

[0122]

number

[0123] This satisfies the inequality in (2). As a result, we can use the approximated MPO representation 22, i.e., instead of M, in the quantum circuit, e.g., c -1 While a renormalized approximated MPO representation such as M may be realized, A denotes the orthogonal approximation 30. From a practical point of view, in some calculations it may be necessary to store this normalization factor c and take it into account at the end of the calculation, for example to interpret the result when it depends on the complex output quantum state.

[0124] In the case of optimization algorithms, a processing system such as a computer k One may seek to find the minimum of the cost function under an isometry constraint on all sub-tensors 32, 32a, 32b, abbreviated as

[0125]

number

[0126] where ||*|| represents an appropriate norm, such as the Frobenius norm, and the optimum is taken over both c and A. This optimization may be performed as an alternating descent method by determining partial derivatives of the cost function with respect to any isometric subtensors 32, 32a, 32b under the isometry constraint, and optimizing an orthogonal approximation 30 based on the partial derivatives, such as by a gradient descent-based algorithm, e.g., stochastic gradient descent, or a momentum-based algorithm such as adaptive moment estimation.

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

[0128] In general, any m×p isometric matrix belongs to a Stiefel manifold.

[0129]

number

[0130] This is a Riemannian manifold. Therefore, the cost function C to be minimized can be defined on a Cartesian product of n Stiefel manifolds. The isometric matrix V that minimizes this cost function is k To find the set of θ, one can use a gradient-based optimization method such as Riemann gradient descent, which works with optimization problems under constraints, as in this case.

[0131] In normal gradient-based methods, such as gradient descent, the variables to be optimized are updated by moving in the opposite direction of the cost function gradient.

[0132]

number

[0133] where α is the step size. In the Riemann generalization of these methods, δ is used to determine the direction of this movement. k Instead of C, point V k Its projection onto the space tangent to the Stiefel manifold in δ k R C can be used. Furthermore, the matrix V can be calculated based on the partial derivatives and / or their projections. k Retraction onto the manifold can be used to update the value of

[0134]

number

[0135] Optimization methods that use information from previous optimization steps, such as gradient descent with momentum, also use the concept of vector translation, which is a generalization of the concept of translation translation and may therefore be implemented directly based on the scheme described above.

[0136] The Riemann formulation of gradient descent may be implemented as an iterative procedure, where at every step one computes the Riemann gradient of the optimization function C at the current point on the manifold, and then finds the next point using a selected retraction in the direction of the negative of the gradient.

[0137] The algorithm can then be implemented as follows: Step 0: Obtain an initial estimate of the orthogonal approximation 30; Step 1: Find a renormalization constant c that minimizes the cost function in equation (3), for example using a fixed orthogonal approximation 30 (i.e., A is constant during this step); Step 2: For the problem in equation (4), find a new approximation A via one step of Riemann gradient descent, where the renormalization constant c may be constant during this step, then return to step 1.

[0138] The iterative algorithm may be terminated at any point, for example, after a predetermined number of iterations or when a predetermined cost threshold is achieved by the cost function.

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

[0140] However, the subtensors 32a, 32b at the boundaries (i.e., V1 and V2, also called boundary cores) n) may be isometric only. However, the boundary cores 32a, 32b may be extended to unitary matrices by extending the sub-tensors under the isometric condition with the addition of |0> states 36 at one end of the tensor network and a projection to |0> states 38 at the opposite end of the tensor network in order to maintain the correspondence with the approximated MPO representation 22. By extending the orthogonal approximation 30, an extended orthogonal approximation 34 can be obtained, which features only the sub-tensors 32 that can be reshaped into unitary matrices.

[0141] Since all subtensors 32 can be reshaped into unitary matrices, the transition to the quantum circuit can be made by unitary operations U1 to U n The isometric subtensors 32 may be implemented by concatenating the isometric subtensors 32 in a cascade of quantum operations 40 implemented as: where the uncontracted indices of the isometric subtensors 32 may correspond to the input and output states of the respective quantum operations 40, while the contracted indices of the isometric subtensors 32 correspond to the input quantum states for the next unitary operation 40 in the concatenated chain of unitary operations 40. Thus, each unitary operation 40 has a number of 1+log(R) qubits input and qubits output, with log(R) qubits being processed by the subsequent unitary operation 40 in the concatenated chain of unitary operations 40. The decomposition of these subtensors 32 into quantum gates may be implemented based on known decomposition techniques of the respective quantum hardware architecture.

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

[0143] Those skilled in the art will appreciate that the actual quantum states prepared or measured for these operations may be modified, for example by preparing based on the |1> state and post-selecting, or by a different state based on an appropriate basis change, as long as the correspondence with the approximated MPO representation 22 is preserved. However, for simplicity, in the following the preparation of the "|0>" state and the corresponding post-selection will be considered.

[0144] Thus, as shown in the example of FIG. 4, quantum circuit 28 may be implemented by the following steps of a quantum processing system.

[0145] Step 1: Set ancilla qubit 42 to the zero state

number

[0146] Step 2: Gates U1, ..., U2 to be decomposed into quantum gate implementations based on the quantum hardware used, such as a combination of one-qubit and two-qubit operations. n are applied in sequence.

[0147] Step 3: Measure the ancilla qubit 44 on a computational basis, and if all measurement results of step 3 are equal to “0”, record the result, otherwise return to step 1.

[0148] Note that for post-selection, one should measure the quantum bit from the other end of the orthogonal approximation 30, i.e., the quantum bit corresponding to the opposite boundary core 32b (see FIG. 4), rather than the quantum bit initially introduced into the initial boundary core 32a as the ancilla quantum bit 42.

[0149] The algorithm implements the quantum circuit 28 corresponding to the approximated MPO representation 22 only if all ancilla qubit measurements 44 are equal to |0>, e.g., a string of zeros is obtained for the corresponding qubits. The state of the quantum system immediately before the measurement determines the likelihood of this "successful" outcome. This means that the initial state of the system to which the quantum method is applied influences this probability.

[0150]

number

[0151] The initial state of the system, ρ in may be determined by the specific computational task. To ignore this dependency, <Pr suc To do this, we can average the probability values ​​over all possible initial states of ρ in can be thought of as a maximally mixed state that is a uniform mixture of states from an orthonormal basis. In this case, <0| a2 U|0> a1 Using =A, ρ in =1 / 2 n I, and the average probability is:

[0152]

number

[0153] That is, the success probability can be determined mainly by the norm of the orthogonal approximation 30. On the other hand, if A is c -1 Since M is approximately equal to λ, we can see how the success probability is related to the normalization factor c, i.e., the higher the normalization factor, the lower the success probability. In other words, there may be a trade-off between the accuracy of the approximation and the success probability. The best theoretical value for the normalization constant c is c = λ. max (M), i.e., the largest eigenvalue of the approximated MPO representation 22. If the approximation rank R of the orthogonal approximation 30 is large enough, the function can be approximated with sufficient accuracy.

[0154] The depth of the quantum circuit 28 scales polynomially with the approximation rank R, and the size of the approximation rank R is proportional to the unitary operation U k Note that this may further increase the complexity of decomposing the isometric subtensor 32 into a combination of quantum gates, since each of these operates on log2R+1 qubits, where R is the rank of the orthogonal approximation A.

[0155] Moreover, a numerical strategy that approximately decomposes general multi-qubit unitary operations40 using a combination of single qubits and CNOT gates can achieve CNOT gate counts close to the theoretical lower bound presented by Shende et al. (Minimal universal two-qubit controlled-not-based circuits, Phys, Rev. A, 2004), i.e., any decomposed m-qubit gate can be calculated as follows: m-1 CNOT gates are needed, so in this case, a total of n*4 logR+1 , i.e., approximately n*R 2 Note that this requires CNOT gates. Thus, in general, the lower the approximation rank R selected for the orthogonal approximation 30, the simpler the quantum circuit 28 will be, but the less accurate the orthogonal approximation 30 will be.

[0156] To demonstrate the power of the method, a frequently used matrix is ​​considered and encoded into the multi-qubit quantum circuit 28. Specifically, the diagonal matrix M D The matrix with the discretized linear function above and the Laplace operator M corresponding to the second derivative matrix after discretization L , both of which are not unitary.

[0157]

number

[0158] n is the number of columns N-1, or N-1 in this example.

[0159] 5 and 6 show the algorithm described in relation to FIG. 2 and FIG. 4 for the diagonal matrix M D and the Laplace matrix M L For any n>1, M D and M L The MPO forms of the matrix have ranks r equal to 2 and 3, respectively.

[0160] The top graph shows the minimum achievable error in approximating the diagonal (FIG. 5) and Laplace (FIG. 6) matrices by the orthogonal approximation 30 as a function of the approximation rank R (in terms of the number of used ancilla qubits 42, which is equal to log2R) for four different sizes of matrix 24 (n is the number of qubits in the graph legend). The bottom graph shows the corresponding average probability of success based on the ancilla qubit measurements 44. The dotted lines indicate the mean probability of success for ||M|| 2 λ max -2 (M)2 -n (see equation (9)) and for both matrices is independent of n for the number n of qubits.

[0161] For both matrices, the method achieves reasonable accuracy (M D and M L It is possible to achieve error rates of about 0.01% and 0.5%, respectively, for 10 qubits and a satisfactory probability of success (above 5%) for up to 50 qubits, which is not a limitation of the method. An important property found in the experiments is that the error and probability of success are nearly independent of the number of qubits, making the method likely to scale well for complex tasks.

[0162] To evaluate the efficiency of the computation, the algorithm can be considered to consist essentially of three steps: finding the orthogonal approximation 30; expanding the boundary cores 32a, 32b so that they can be reshaped into a unitary matrix; and decomposing the resulting 2R*2R unitary matrix into a series of 1- and 2-qubit gates.

[0163] The time required for the last step is determined by the decomposition method used and will not be described here. Furthermore, in this example, the expansion step is skipped and the kernels V1 and V2 are calculated based on the initial estimate of the orthogonal approximation 30. n Instead of optimizing itself, U1 and U n The expanded sub-tensor 32 corresponding to is directly optimized in an iterative optimization algorithm.

[0164] Therefore, the complexity of the algorithm for determining the quantum circuit 28 may be determined primarily by the complexity of the iterative optimization process. The execution time of the optimization procedure 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 several 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 in advance the number of sufficient iterations, and below we focus on the complexity of one iteration of the optimization procedure, which consists of the following elements:

[0165] Element 1: The calculation of the normalization constant c is based on the computational complexity of multiplication while simultaneously taking the traces of two MPOs with ranks r1 and r2, which has computational complexity of O(nr1r2max(r1,r2))$. 3 ) computational complexity scaling.

[0166] Element 2: The calculation of the cost function is O(n(r+R) 3 After the sum of the two MPOs, their ranks are added and the resulting MPO of rank r1+r2 is deflated by its complex conjugate, which can have a computational complexity of O(n(r1+r2) 3 ) computational complexity.

[0167] Element 3: Evaluation of partial derivatives of the cost function using automatic differentiation. When using automatic differentiation, its algorithmic complexity is theoretically guaranteed to be less than or equal to that of the original program, which is equivalent to the preceding element (2) and takes O(nR 3 ).

[0168] Element 4: The optimization process includes a Riemannian optimization step, which k For a 2R*2R matrix, the first and third steps of the Riemann optimization step are O(R 3 ) The complexity of the second step depends on the retraction method chosen. In the exemplary implementation, it is O(R 3 ) complexity. Thus, the Riemannian optimization step for each iteration is O(nR 3 ) scaling.

[0169] 7A and 7B show an example of the diagonal matrix M D 2 and 4 to determine a quantum circuit 28 when approximating

[0170] FIG. 7A shows the running time of one iteration of the iterative optimization algorithm for different approximation ranks R of the orthogonal approximation 30 as a function of the number of qubits n, which depends on the size of the matrix 24 .

[0171] FIG. 7B shows the experimentally determined execution time of each iteration and its components (1)-(4) (as explained in the sections above) as a function of the number of qubits n.

[0172] In full agreement with theoretical estimates of the complexity of the algorithm, these execution times depend polynomially on the number of qubits, n, and therefore on the target matrix size of matrix 24, N=2 n It depends polynomially on

[0173] Thus, the present method for determining quantum circuit 28 promises to provide an efficient algorithm for determining quantum gate decompositions of matrices in quantum computing applications. Note that the depth of quantum circuit 28 scales only linearly with the number of qubits n, and therefore scales logarithmically with the matrix size of matrix 24.

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

[0175] Additionally, the orthogonal approximation 30 encoded as quantum circuit 28 performs quantum operations that can be applied multiple times to the initial quantum state to perform higher orders of the approximated MPO representation 22, such as to find the maximum eigenvalue of a matrix that approximates a function by measuring the output state of the state after the higher orders of the approximated MPO representation 22 act on the input state.

[0176] The task of finding the eigenvalues ​​of a diagonalizable matrix can be thought of as a search for the maximum element in a tensor (a tensor of values ​​on the matrix diagonal), and therefore the power method can be applied to almost any optimization problem after the corresponding preprocessing, where a functional representation is processed to obtain an MPS / MPO approximation.

[0177] Let f(x) represent some real-valued function that we wish to maximize / optimize over the input domain, x ∈ [a, b]. Given access to n qubits, we can calculate the number of iterations required to represent f(x) as a quantum state: n discretization points are available.

[0178]

number

[0179] where |k> is the computational basis state of n qubits, and by normalization, |<f|f> | 2 = 1. It turns out that the maximization problem is as follows:

[0180]

number

[0181] Here, p k =|<k|f> | 2 is the measurement probability of state |f with respect to the projection onto the computational base. Therefore, if such a state can be prepared, by measuring it it is possible to recover an approximate optimum value of f(x) up to the error introduced by the discretization. Power iteration can improve on this idea. A power iteration considering a diagonal matrix with function values ​​as diagonal elements looks like this:

[0182]

number

[0183]

number

[0184] where ∥.∥ is the 2-norm. Ideally, as j-→∞, the resulting state is |k * >,k * =arg max k p k and thus recovering the optimal value. To practically realize this simple technique on a quantum computer, a quantum method can be used to realize the power iteration as in (13), hereafter also called the "quantum power method".

[0185] In classical computers, matrix multiplication with large dimensions can be computationally expensive, but in quantum computers, the size of the matrices can be exponentially compressed due to the large computational state space, making this method viable for optimization techniques.

[0186] FIG. 8 illustrates a schematic flow chart for determining a quantum computer implementation of an algorithm for finding the maximum eigenvalue of a matrix that approximates a function, according to one example.

[0187] In the first step shown in FIG. 8(A), the corresponding input value x of the function f(x) is l , x r Function value of

number

[0188] Based on the function values, the MPS state is determined according to a suitable algorithm, for example based on an iterative truncated SVD procedure or a cross approximation. As shown in FIG. 8(B), the function values ​​of the evenly spaced sampled functions may form a function value vector, based on which the values ​​of the MPS state vector may be determined. The MPS states are up to 2nr 2 2 of a function-valued vector with elements n It can represent values.

[0189] However, the MPS approximation cannot be directly used in classical algorithms to implement power methods because the rank of MPS grows exponentially with the power of the power method. Thus, simple classical realizations of power methods quickly go out of range because MPS rank scaling exponentially with iterations.

[0190] Nevertheless, based on the MPS representation of the function (or the function itself), which may approximate the function based on discretized function values, an approximated MPO representation 22 of the function may be obtained as shown diagrammatically in FIG. 8(C).

[0191] Using the optimization described in connection with Figures 2 to 4 above, a (now expanded) orthogonal approximation 30 of the approximated MPO representation 22 may then be obtained, as shown in Figures 8(C) and 8(D).

[0192] In the orthogonal approximation 30 shown in FIG. 8(D), a corresponding quantum circuit 28 is determined, which can be constructed from unitary operators 40 for each isometric subtensor 32 obtained by the optimization algorithm.

[0193] For an implementation of the quantum exponentiation algorithm, quantum circuit 28 operates multiple times on the input states of the qubits in qubit register 12, as shown in a schematic example in FIG. 8(E).

[0194] The input states of the qubits may be obtained based on state initialization quantum circuit 46, which may include application of a Hadamard gate to each of the qubits starting in a |0> basis state, thereby initializing a "vector of ones" in the computational basis spanned by the qubits. Quantum circuit 28 then operates on the qubits multiple times to obtain a final result state that may be measured using the detectors of measurement module 18. Quantum exponentiation iterations can be achieved by repeating quantum circuit 28.

[0195] As shown diagrammatically in FIG. 8(F), the resulting result 48 is measured (e.g., multiple times) to obtain the output of the quantum exponentiation algorithm in the form of the maximum eigenvalue of the effective matrix implemented by the quantum circuit 28 based on the orthogonal approximation 30.

[0196] As shown in detail in the example of Figure 9A, for each power p of the exponentiation scheme, the algorithm implementation includes one quantum circuit 28a-28p that is applied sequentially to the qubits. Additionally, for each quantum circuit 28a-28p, a corresponding ancilla qubit 42 is initialized and installed. Although the figure shows a single |0> state, one skilled in the art will appreciate that a log(R) ancilla qubit 42 may be prepared for each quantum circuit 28a-28p as described above.

[0197] In some examples, for each quantum circuit 28a-28p, the measured qubit 44 may be re-initialized and used as the ancilla qubit 42 for the subsequent quantum circuit 28a-28p, if a sufficiently fast measurement and initialization technique may be available. In other examples, a new ancilla qubit 42 is added for each quantum circuit 28a-28p.

[0198] The algorithmic procedure in implementing the quantum exponentiation method can be implemented as follows (represented diagrammatically in FIG. 9A).

[0199] 1. Prepare the initial state:

number

number

[0200] 2. According to the isometric subtensor 32 of the extended orthogonal approximation 30, the quantum operation 40, U k Quantum circuits decomposed into 28,U MPO applies.

[0201] 3. Until you get the power p you need,

number

[0202] 4. Measure all ancilla qubits44 on a computational basis.

[0203] 5. Measure Ancila qubit 44. The measurement result is a bit string

number

[0204] 6. Repeat steps 1-4 to collect measurement statistics and infer the candidate optimum as the most likely result.

[0205] With post-selection on an ancillary space, algorithmic success can be shown to depend on the post-selection probability, which can be shown to be independent of the number of computational qubits (n), and can be proportional to the power of the power law by 1 / p. Taking all steps into account, the running time of the quantum method is t ∝ n × R 2 ×p. The complexity is only linear in p, not exponential with respect to classical (tensor network) implementations of the power method. Thus, in this case, quantum methods can offer an exponential run-time advantage.

[0206] As shown in FIG. 9B, a quantum circuit for implementing the quantum exponentiation method can be efficiently compressed by executing some of the quantum operations of quantum circuits 28a to 28p in parallel.

[0207] In this example, four quantum circuits 28, each including four quantum operations 40, U1-U4, are concatenated based on corresponding isometric subtensors 32 of the orthogonal approximation. However, the second quantum operation 40, U2, of the first quantum circuit 28a may be executed in parallel with the first quantum operation 40, U1, of the second quantum circuit 28b. Thus, the execution time of the entire quantum operations 28a-28d may be less than the sum of the execution times of the individual constituent quantum circuits 28a-28d.

[0208] In some examples, for each power of the quantum exponentiation method, the execution time of the algorithm on the quantum computer may increase by the execution time of one quantum operation 40 corresponding to the last isometric sub-tensor 32b of the orthogonal approximation 30.

[0209] In some examples, the first output of the first quantum operation 40 of the first quantum circuit 28a in the concatenated sequence of quantum circuits 28a-28p is the second input for the first quantum operation 40 of the second subsequent quantum circuit 28b in the concatenated sequence of quantum circuits 28a-28p, and the second output of the first quantum operation 40 of the first quantum circuit 28 is the first input for the second quantum operation 40 of the second circuit 28b. In some examples, as shown in FIG. 9B, the first output of the second quantum operation 40 of the first quantum circuit 28a in the concatenated sequence of quantum circuits 28a-28p is the second input for the second quantum operation 40 of the second subsequent quantum circuit 28b, and the second output of the first quantum operation 40 of the second quantum circuit 28 is the first input for the second quantum operation 40 of the second circuit 28b.

[0210] 10A-10D show example measurement probabilities obtained for quantum circuit encoding and simulated quantum power law of a two-dimensional Ackley function as contour plots, with the highest probabilities in the center of each graph.

[0211] Figures 10A and 10B show measurement probabilities for measuring particular states for an Ackley function implemented with an orthogonal approximation 30 of rank R=8 and R=4, respectively. Figures 10C and 10D show higher powers of an orthogonal approximation 30 of rank R=4 implemented as a quantum circuit 28 on a simulated quantum computer, with Figures 10C and 10D showing examples of power iterations with p=10 and p=25, respectively.

[0212] The two-dimensional Ackley function is defined as follows:

[0213]

number

[0214] Consider the domain x,y ∈ [-5,5], the function reaches a global minimum f(0,0) = 0. Finding this global minimum is a well-known non-convex optimization problem that poses difficulties for gradient-based optimizers. Hence, the function is commonly used to benchmark optimization algorithms. Here, we consider quantum power iteration as a method to find the global minimum. To do this, we first convert the minimization problem into a maximization problem by standard inversion techniques. Specifically, g(x,y) = 1-tanh(κf(x,y)).

[0215] It should be noted that for an appropriate choice of κ, minimizing the Ackley function as in equation (15) is equivalent to maximizing g(x,y). In the following, for simplicity, this modified function g(x,y) is called the Ackley function. In the example, a total of n=10 computational qubits (5 qubits per dimension) are used for the domain discretization, and the MPS representation is obtained based on the discretized values ​​with a maximum rank r=6 using the cross approximation method in TTPy to obtain the MPS core. From the MPS core, the diagonal MPO, H, as the approximated MPO representation 22 of the function is obtained in a straightforward manner by known algorithms.

[0216] Next, an orthogonal approximation 30 of the approximated MPO representation 22, U * is obtained by minimizing the cost function as in equation (4). The algorithm is used to search for an isometric MPO core 32, i.e., an isometric subtensor 32, with approximate rank R=8 and R=4, and equation (4) can indeed be minimized by Riemann gradient descent with a learning rate of 0.02 and a total of 10,000 iterations.

[0217] Orthogonal approximation with approximation rank R = 8, U * can be cast to a Qiskit quantum circuit, but requires a compilation step that decomposes the 4-qubit operations into a sequence of 2-qubit and single-qubit gates. The resulting circuit depth is therefore relatively large. To reduce simulation time on a classical computer, we perform a minimization of equation (11) with approximation rank R = 4. Although convergence cannot be guaranteed in this case, such an approach is still motivated by the fact that the power method finds the maximum prescribed by equation (14). Although the minimization of equation (4) may tolerate large errors due to approximation rank deficiency, if the maximum is solved exactly, the orthogonal approximation 30, U * can be further used for power iteration.

[0218] In Figures 10A and 10B, orthogonal approximations 30 of the Ackley function with different ranks R=4 and R=8, respectively, are compared, and it is observed that: 1) the orthogonal approximation 30 of approximation rank R=8 qualitatively recovers the exact function, and 2) the orthogonal approximation of approximation rank R=4 solves the optimal value of the exact function.

[0219] The R=4 orthogonal MPO is cast to a Qiskit quantum circuit and exponentiation up to 25 is realized (see Fig. 10C, D). The results are shown in Fig. 10C, D, where the quantum exponentiation method indeed converges to the optimum value of the Ackley function at (0,0).

[0220] Thus, a low-rank tensor network approximation 22 of some given function can be prepared and used to obtain an orthogonal approximation 30 through a variational optimization algorithm running on a classical computer, which can be cast to a quantum circuit 28 having a small qubit count and consisting of multiple shallow quantum operations using ancilla qubits 42, 44. Exponentiation iteration can be achieved by repeated concatenation of such quantum circuit blocks 28a-28p. The linear scaling of quantum exponentiation with power p indicates that quantum algorithms can be efficiently utilized on near-term quantum hardware to solve optimization problems that are difficult on classical hardware.

[0221] The above description primarily refers to functions for which extrema are sought. However, those skilled in the art will appreciate that the function may also be provided in matrix form, and thus the method may effectively search for the highest eigenvalue of the matrix as a function. While the illustrated example relates to an equidistant sampling of the function, 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 used as part of the method of FIG. 2 as an approximated MPO representation 22 of the function.

[0222] Additionally, although the illustrated example shows that the quantum bit measured to determine outcome 48 is different from the computation quantum bit initially prepared in the |+> state, one skilled in the art will understand that in some examples, the same quantum bit prepared in the |+> state is also measured to determine outcome 48, for example, through the use of an additional gate, e.g., a SWAP gate, to swap the quantum states of the quantum bits.

[0223] The description of the preferred embodiment and the drawings merely serve to illustrate the invention and its associated beneficial effects, and should not be understood as implying any limitation, the scope of the invention being to be determined solely by the appended claims. [Explanation of symbols]

[0224] 10. System 12 qubit register 14 State Preparation Module 16 Quantum Matrix Multiplication Module 18 Measurement Module 20 Processing System 22 MPO approximation 24 Matrix Representation of Functions 24a Penrose diagram representation of a matrix 26 Subtensor 28 Quantum circuit 30 Orthogonal Approximation 32, 32a, 32b Isometric Subtensors 34 Extended Orthogonal Approximation 36,38 Boundary core |0> state 40 Quantum Computation 42 Ancila Qubit 44 Measured Ancillar Qubits 46 State Initialization Quantum Circuits 48 results

Claims

1. 1. A computer-implemented method for solving an optimization problem given a function (24, 24a), the computer-implemented method comprising: obtaining an approximation of said function (24, 24a) in the form of an approximated matrix product operator (MPO) representation (22) of said function (24, 24a); selecting an approximate rank based on said approximated MPO representation (22); determining a rank-based orthogonal approximation (30, 34) of said approximated MPO representation (22) in the form of an orthogonal MPO having isometric subtensors (32, 32a, 32b) of said approximate rank; encoding the orthogonal approximation (30, 34) into a quantum circuit (28, 28a-p) based on encoding the isometric subtensors (32, 32a, 32b) into quantum gates that operate on an initialized set of computation qubits and ancilla qubits (42); encoding an arbitrary power of the orthogonal approximation (30, 34) based on applying the quantum circuits (28, 28a-p) to an input quantum state multiple times in a concatenated sequence of the quantum circuits (28, 28a-p); Including, method.

2. the approximated MPO representation (22) is an approximate tensor network representation of the function (24, 24a), the method comprising, in particular, the step of obtaining a set of sub-tensors (26) via cross approximation of the function (24, 24a), The method of claim 1.

3. The method further comprises obtaining the approximated MPO representation (22) based on a plurality of sample values ​​of the function (24, 24a) for a corresponding plurality of different input values. The method according to claim 1 or 2.

4. obtaining the approximated MPO representation (22) comprises obtaining a matrix product state (MPS) representation of the function (24, 24a) by a cross approximation technique, and obtaining the approximated MPO representation (22) based on the MPS representation by enforcing diagonal constraints. The method according to claim 1 or 2.

5. Obtaining the orthogonal approximation (30, 34) in the form of the orthogonal MPO, determining an initial estimate of the orthogonal approximation (30, 34) having an isometric subtensor (32, 32a, 32b) of the approximation rank; starting from the initial estimate, iteratively optimizing the orthogonal approximation (30, 34) based on an optimization algorithm that minimizes a cost function subject to an isometry constraint of the isometric subtensors (32, 32a, 32b), the cost function attributing a cost to the orthogonal approximation (30, 34) based on its quality with respect to the approximated MPO representation (22); Including, The method according to claim 1 or 2.

6. the step of applying the quantum circuits (28, 28a-p) a plurality of times in a concatenated sequence of quantum circuits (28, 28a-p) comprises the step of preparing a set of a number of ancilla qubits (42) in an initial state, said number being at least equal to the number of times the quantum circuits (28, 28a-p) are applied; Each set of ancilla qubits (42) includes, in particular, a number of qubits equal to the logarithm of the approximation rank, The method according to claim 1 or 2.

7. each subsequent quantum circuit (28a-p) in the concatenated sequence of quantum circuits (28, 28a-p) operates on the set of ancilla qubits (42) and the set of computation qubits operated on by the previous quantum circuit; The method according to claim 6.

8. Each time the quantum circuit (28, 28a-p) is applied, the method further comprises: measuring a set of post-selected qubits (44) on which each quantum circuit (28, 28a-p) acts and which are not acted upon by a subsequent application of the quantum circuits (28, 28a-p) in a concatenated sequence of the quantum circuits (28, 28a-p); Each set of post-selected qubits (44) specifically includes a number of qubits equal to the logarithm of the approximate rank, The method according to claim 1 or 2.

9. preparing the initial quantum state includes applying a Hadamard gate to each qubit of the set of initial computation qubits; The method according to claim 1 or 2.

10. applying the quantum circuits (28, 28a-p) in a concatenated sequence of quantum circuits (28, 28a-p) p times to an input quantum state includes providing p sets of ancilla qubits (42), each set of ancilla qubits (42) being acted upon by a respective one of the quantum circuits (28, 28a-p) and including a number of ancilla qubits (42) equal to the logarithm of the approximate rank; The method according to claim 1 or 2.

11. A processing system (20) for solving an optimization problem given a function (24, 24a), said system (10) comprising: obtaining an approximation of said function (24, 24a) in the form of an approximated matrix multiplication operator MPO representation (22) of said function (24, 24a); selecting an approximate rank based on said approximated MPO representation (22); determining a rank-based orthogonal approximation (30, 34) of said approximated MPO representation (22) in the form of an orthogonal MPO having an isometric subtensor (32, 32a, 32b) of said approximation rank; determining an implementation of the orthogonal approximation (30, 34) in a quantum circuit (28, 28a-p) based on the encoding of the isometric subtensors (32, 32a, 32b) into quantum gates acting on computational qubits and an initialized set of ancilla qubits (42), in which the quantum circuit (28, 28a-p) is applied to an input quantum state multiple times in a concatenated sequence of quantum circuits (28, 28a-p); configured to implement A processing system (20).

12. To obtain the orthogonal approximation (30, 34) in the form of the orthogonal MPO, the system (10) determining an initial estimate of an orthogonal approximation (30, 34) of said function (24, 24a) in the form of a tensor network having isometric subtensors (32, 32a, 32b) of said approximation rank; starting from the initial estimate, iteratively optimizing the orthogonal approximation (30, 34) based on an optimization algorithm that minimizes a cost function subject to an isometry constraint of the isometric subtensors (32, 32a, 32b), the cost function attributing a cost to the orthogonal approximation (30, 34) based on a quality of the orthogonal approximation (30, 34) with respect to the approximated MPO representation (22); configured to implement The system (10) of claim 11.

13. the system (10) is further configured to communicate the implementation to a quantum computing system for executing the implementation on quantum hardware. A system (10) according to claim 11 or 12.

14. A hybrid quantum-classical computing system comprising the system (10) of claim 11 or 12 and quantum computing hardware, the hybrid quantum-classical computing system comprising: receiving said embodiment from said processing system (20); Implementing the embodiment on the quantum computing hardware; and receiving a computation result from the quantum computing hardware; configured to implement A hybrid quantum-classical computing system.

15. When executed by a processing system (20), it implements the method according to claim 1 or 2 or the system according to claim 11 or 12. A non-transitory medium containing machine-readable instructions.

Images