SYSTEM AND METHOD FOR SOLVING A DIFFERENTIAL PROBLEM USING A QUANTUM COMPUTER

Abstract

 Methods and systems are disclosed for solving a differential problem using a hybrid computer system. The differential problem defines a differential equation and one or more constraints, and solving the differential problem comprises finding a target function defining an approximate solution of the differential equation and obeying the constraints. The target function is represented by a vector of coefficients with respect to a given function basis. The method comprises receiving or determining information on one or more quantum circuits encoding a Hamiltonian representing the differential problem with respect to the basis functions; and determining the solution to the differential problem. The determining comprises encoding the quantum circuits using the quantum processor; and determining an eigenvector of the Hamiltonian with a minimal eigenvalue, preferably the eigenvalue being equal to zero, the eigenvector representing the vector of coefficients of the target function.

Description

Technical field

[0001] The disclosure relates to solving differential problems, and in particular, though not exclusively, to methods and systems for solving differential problems using a hybrid computer system, and to a computer program product enabling a hybrid computer system to perform such methods.

Background

[0002] Solving differential equations (DEs) is a crucial task in many fields of science, engineering and economics. Differential equations provide a concise and powerful mathematical framework for formulating the dynamics of real-world phenomena. Because closed forms of the solution are often unavailable for practical problems, several classical techniques (e.g., finite-difference methods (FDM), finite-element methods (FEM) and spectral methods) have been developed for numerically solving differential equations to obtain approximate solutions. Many variations of these kinds of algorithms are known in the art.

[0003] Finite-difference methods break down a problem domain into a regular gird of discrete pieces, approximating differential terms by finite differences at each grid point. A numerical solution is typically found by collecting the differential equations and data constraints (e.g., boundary conditions, initial conditions, or fixed points) into a single matrix equation that can be solved by matrix inversion techniques. Finite-difference methods need to balance the resolution of the solution with the grid size, which can grow very quickly with the dimension of the problem, and which is practically speaking often limited by hardware constraints such as memory size, cache memory, et cetera.

[0004] Instead of a regular grid, finite-element methods use a mesh of elements to model the problem domain. Various solving techniques may be used; some popular methods employ so-called 'weak' formulations that describe the differential equations in integral form. The approximate solution can then be expressed as a linear combination of shape (basis) functions that are applied independently across each element. A system matrix equation is produced (and solved) after assembling individual local equations along with any data constraints.

[0005] In spectral methods, the solution is approximated by a finite set of orthogonal basis function (e.g., Fourier series or trigonometric functions for periodic problems, or Chebyshev or Legendre polynomials for non-periodic problems). This converts the problem into finding optimal coefficients of the truncated expansion of the sought solution through solving a matrix equation. Because basis function spans the entire domain, spectral methods provide a global (non-local) approximation and discretization is not required, making it distinct from finite-difference methods and finite-element methods.

[0006] Typically, these techniques can handle nonlinear differential equations (NDEs) using iterative methods. The order of accuracy can be increased by increasing the number of elements and/or higher-order shape function for finite-difference and finite-element methods, or by increasing the number of expansion terms for spectral methods. However, these introduce computational difficulties when tackling complex problems in high-dimensional spaces.

[0007] Alternatively, the increasingly popular field of machine learning has also been turned to the task of solving differential equations. This is part of the field known as SciML (scientific machine learning), where machine-learning models are tasked with learning relationships from scientific data. The solving of differential equations generally utilises PINNs (physics-informed neural networks) - deep neural networks that include differential constraints introduced with automatic differentiation. PINNs have been applied to a wide range of problems. However, it is still a developing field with issues left to address such as training guarantees and cost per epoch.

[0008] Quantum computing offers a fundamentally different way to compute. Correspondingly, certain processes and techniques can efficiently be implemented on quantum devices (or hybrid quantum computers), but not on classical computing devices. The expectation is that by utilizing these areas of efficiency, algorithms for problems of interest can be found and implemented which offer a speed up over classical methods. The solving of (certain classes of) differential equations offers one potential area for advantage, and consequently, there has been effort in creating quantum algorithms for this use using a variety of approaches. Many of these approaches are analogous to the approaches for classical algorithms described above.

[0009] For example, in quantum finite-differencing approaches, the problem is discretised into a system of equations and then solved. However, the absence of, e.g., register shifts makes implementing these algorithms on current hardware very challenging. Also proposed are spectral methods, where the solution is fit by finding coefficients for a given basis. Other proposed works can be categorized as SciML approaches using quantum machine learning, where a quantum model is variationally taught to respect the known differential constraints - similar to PINNs. However, variational methods tend to require a large number of circuit evaluations, which quickly leads to (very) long computation times given the repetition rate of current quantum hardware.

[0010] It is classically a hard problem to generate random samples from a given probability distribution. This is especially the case for multivariate problems. Because quantum systems are inherently probabilistic, generative modelling is one of the areas where quantum computers can provide an advantage. This is particularly true for multivariate problems, which can be associated with measurements on entangled quantum states.

[0011] Quantum generative modelling (QGM) aims to exploit trainable circuits that can prepare distributions as quantum states. Being a subject of emerging quantum machine learning (QML), QGM utilizes the Born rule inherent to quantum mechanics. The goal is to represent a parametrized probability distribution p_θ(x) as a probability to measure a bit string x from a variational state parametrized by a vector gate parameters θ. By using so-called quantum circuit Born machines (QCBMs) models can be learned directly from samples of a target distribution using various loss functions. Here, directly means that a sample can be determined directly from a measured wavefunction, so that for one sample in principle only one quantum measurement is needed. However, while making a step towards sampling-based advantage on quantum computers, QGM performance is largely limited, both by poor trainability and due to undefined or poorly defined derivatives. This limits in particular applications such as solving stochastic differential equations (SDE). Alternative solutions that have been proposed in the art rely on computationally-expensive amplitude encoding schemes, complex or deep quantum circuits, or other incompatibilities with realistic NISQ hardware devices available today and in the expected near future.

[0012] WO 2023/158404 A1 describes a way to sample probability distribution functions based on a variational approach, e.g., to solve a stochastic differential equation. While this approach addresses many of the above-mentioned drawbacks, such a variational approach is still relatively sensitive to noise.

[0013] Hence, there is a need in the art to quickly and efficiently solve differential problems using hybrid quantum computers.

Summary

[0014] It is an aim of the embodiments in this disclosure to provide a system and method for solving differential equations that avoids, or at least reduces the drawbacks of the prior art.

[0015] In a first aspect, embodiments in this disclosure relate to a method for solving a differential problem using a hybrid computer system comprising a classical computer system and a quantum processor. The quantum processor may comprise a plurality of quantum elements, e.g., qubits, qutrits, or other qudits. The differential problem defines one or more differential equations and, optionally, one or more constraints. Solving the differential problem comprises finding a target function defining an approximate solution of the differential equation, the target function being represented by a vector of coefficients with respect to a given function basis. The method comprises receiving or determining information on one or more quantum circuits encoding a Hamiltonian representing the differential problem with respect to the basis functions, and determining the solution to the differential problem. The Hamiltonian comprises a first contribution representing the differential equation and, optionally, a second contribution representing the one or more constraints. The determination of the solution to the differential problem comprises: encoding the one or more quantum circuits using the quantum processor; and determining an eigenvector of the Hamiltonian by finding a ground state, or zero-eigenvalue state, of the Hamiltonian, which eigenvector represents the vector of coefficients of the target function. More in particular, the determination of the solution to the differential problem may comprise: initialising at least a first part of the quantum processor (e.g., a first plurality of quantum elements) in an initial state, e.g., an initial fiducial state; evolving the at least a first part of the quantum processor into a final state, the final state targeting the ground state of the Hamiltonian; and receiving measurement data from the quantum processor, e.g., from measuring at least some of the plurality of quantum elements, to obtain the solution to the differential problem. The evolving the plurality of quantum elements may comprise applying a sequence of pulses to at least part of the plurality of quantum circuits, the sequence of pulses encoding the one or more quantum circuits.

[0016] Such a Hamiltonian effectively encoding a differential problem may also be referred to as an effective Hamiltonian.

[0017] The solution may be constrained by the constraints, e.g., boundary constraints, initial constraints, internal fixed points, and/or data driven constraints. In some cases, the solution may also be considered a fitted function to the constraints, wherein the fit is constrained by the differential equation. The described method allows to encode the solution in a quantum register of the quantum processor. The thus encoded solution may be evaluated at one or more points, or may be used as a starting point for further operations, as will be discussed in more detail below.

[0018] An advantage of this method is that the representation of the target function is a ground state of the Hamiltonian, and is therefore inherently stable: the quantum elements (typically qubits) cannot make an unwanted transition to a lower-energy state. Moreover, small deviations around the ground state remain close to the desired solution. This makes a function encoded in a ground state suitable for further computations; this includes not only determining function values of the target function at one or more points of interest, but also analysis of properties like convexity of the target function, using the target function as input for a further model, et cetera.

[0019] As a further example, the described method are particularly suitable for modelling ongoing processes, wherein the Hamiltonian can be updated in a relatively straightforward manner by adding or replacing one or more data constraints, and allowing the system to transition to the new ground state.

[0020] The target function may be defined as an overlap of a basis state representing the basis functions of the function basis, and the eigenvector representing the coefficients of the target function. The basis state is dependent on the one or more independent variables of the differential problem, whereas the coefficients vector is independent of the one or more independent variables. That is, the Hamiltonian may be formulated with respect to some function basis 〈τ(x)|; then, if the Hamiltonian's eigenvector corresponding to the ground state is represented by |ψ〉, then the solution to the differential problem is given by 〈τ(x)|ψ〉. In general, for a system using n encoding qubits, N = 2ⁿ basis functions can be used to represent the solution.

[0021] In brief, the method employs collecting individual latent contributions of the differential equation and constraints, to form an effective Hamiltonian operator

. An energy function ε(η) can be constructed as an expectation value of

with respect to an unknown quantum state |ψ〉 to yield a scalar energy. In principle, there are various options to find a solution to such a problem. For example, the problem can be treated as a conventional machine-learning task with a cost function corresponding to the energy, using known quantum machine-learning techniques, e.g., by variationally tuning N = 2ⁿ independent parameters {c_k} for k = 0, 1, ... , N - 1, related to each component of |ψ〉 (for a system with n encoding qubits).

[0022] Alternatively, a solution can be found by exploring the eigenproperties of the effective Hamiltonian operator

, which can be constructed as a function of a single variable η > 0. This approach does not require a variational Ansatz, thus avoiding common training issues, such as barren plateaus and local minima trapping, that are generally applicable to variational protocols. Moreover, as a spectral method, there is no dependence on a grid of collation points, making the method scalable to higher dimensions.

[0023] In general, the Hamiltonian

can be represented by a real symmetric matrix with distinct real non-negative eigenvalues and mutually orthonormal eigenvectors. Thus, the effective Hamiltonian operator can be written as

, where Q is a matrix where the k^th column is the k^th eigenvector q_k, and Λ is a diagonal matrix containing all eigenvalues λ_k. This is also known as the eigenspace representation. In this eigenspace representation, ε(η) can be expressed as a sum of N squared projection lengths

weighted by the eigenvalue λ_k associated with eigenvector q_k. The diagonal entries of A control the definiteness of ε(η), and the total energy is always nonnegative and (at least approximately) equal to zero for solving the problem. This implies that

must be a positive-semidefinite matrix for certain

with k = 0,1, ..., N. A positive-semidefinite

with exactly one eigenvalue being zero can be achieved when

, where

is the largest value of η in

. The energy minimization requires that |ψ〉 aligns with

, the eigenvector (ground state) associated with the zero eigenvalue (

) of

, and is perpendicular to the other eigenvectors (

for k = j), thus solving the differential problem (up to a normalisation factor) without the use of variational methods.

[0024] This solution can be obtained in a fully analogue manner, avoiding some of the issues associated with limited qubit correlation times and lack of error correction on current-era and near-term quantum systems. As the effective Hamiltonian can be encoded using digital gates (e.g., specific programmable pulse sequences), a general-purpose quantum processor can be used to solve the differential problem. Alternatively, special-purpose quantum processors may be designed to solve specific differential problems in an efficient manner.

[0025] Depending on the context, the value of the normalisation factor may or may not be required. For instance, if the solution is a probability density function, the normalisation factor is equal to 1 by definition. If one is only interested in the 'shape' of the solution (e.g., convexity, increasing or decreasing parts, linear and non-linear parts, et cetera), or if the problem is defined in arbitrary units, the value of the normalisation factor is immaterial.

[0026] However, if computation of the normalisation factor is required, there are several potential solutions. For instance, if a non-zero function value is known for some input value, one may evaluate the function for that input value and derive the normalisation factor from the ratio between the known function value and the evaluated value. The normalisation factor can also be computed without evaluating the quantum solution (i.e., using only a classical computer), by determining the eigenvalue corresponding to the ground state of the matrix representation of the effective Hamiltonian. In that case, the task of searching for

corresponds to solving a set of linear equations

, indicating that the scalar value of

can be specifically located without the need of any optimization procedures. It is noted that numerical methods are known in the art for finding a smallest eigenvalue using a classical computer, without needing to determine all eigenvalues. By substituting

and the corresponding ground state

into the quantum model, the solution to the differential problem is obtained. However, it is noted that in general, it is not needed to (explicitly) determine the smallest eigenvalue and corresponding eigenvector.

[0027] Thus, the method may be considered a method for encoding a solution to a differential equation in an eigenstate of a Hamiltonian of a quantum processor. Such an encoding in a ground state is relatively stable even for entangled states, and resilient to perturbations. Excited states will exhibit similar behaviour, and will tend to be at least close to solutions of the differential problem, obeying constraints posed by the differential equation and the boundary constraints. This is in contrast to solutions of variational problems, where a small change in the encoding state can correspond to a large difference in the encoded function. Moreover, such an encoding is generally well scalable to, e.g., larger qubit registers.

[0028] The encoded solution can then be used to study various properties of the solution, e.g., as described in C.A. Williams et al., 'Addressing the Readout Problem in Quantum Differential Equation Algorithms with Quantum Scientific Machine Learning', arXiv:2411.15259v1 [quant-ph], which is hereby incorporated by reference in its entirety. For example, the convexity of the solution may be studied; if the solution represents a state of a physical system, that state may be classified (e.g., running properly or being at fault), et cetera. In general, if the problem is cast as an optimisation problem, a quantum advantage may be expected in cases where the model is trainable, i.e., if the resulting optimisation problem is convex and the optimisation landscape sufficiently smooth.

[0029] Thus, as an example, a probability distribution (e.g., a probability density function) may be loaded into the quantum register, where the probability distribution is a solution to a differential problem (for example, a stochastic differential equation). Hence, in an embodiment, the solution to the differential problem represents a probability density function and the method further comprises generating, by the quantum processor, samples based on the probability density function. The generation of samples may comprise applying a unitary transformation associated with the function basis, and measuring the quantum register in a computational basis associated with the function basis. This measurement may yield bitstrings as measurement results, each bitstring representing a sample based on the probability density function. Methods to generate samples from a probability distribution encoded in a plurality of quantum elements is described in more detail in WO 2023/158404 A1, which is hereby incorporated by reference in its entirety.

[0030] It is noted that once the solution is encoded in at least some of the quantum elements, the solution may be scaled to a larger number of quantum elements using known methods. This is e.g. the procedure used in qubit-wise learning for O. Kyriienko et al., 'Protocols for trainable and differentiable quantum generative modelling', Phys. Rev. Research 6 (2024) 033291, but also reminiscent of classical methodologies for super-resolution, as mentioned e.g. in Ren, Pu, et al. "PhySR: Physics-informed deep super-resolution for spatiotemporal data." Journal of Computational Physics 492 (2023): 112438. This way, for example, the resolution of generated samples may be increased. Thus, in an embodiment, the solution (the final state) is encoded in a first selection of quantum elements from the plurality of quantum elements. In such an embodiment, the method may further comprise transferring, e.g., scaling, the solution to a second selection of quantum elements from the plurality of quantum elements, the second selection of quantum elements containing a different number, typically a larger number, of quantum elements than the first selection of quantum elements.

[0031] In an embodiment, the method further comprises evaluating the target function. Evaluating the target function may comprise encoding an independent variable of the target function in a latent space of the quantum processor state by initialising a second part of the quantum processor in a fiducial state and evolving the fiducial state with a feature map associated with the given function basis to obtain an encoded state, the feature map being based on the independent variable; and measuring an overlap of the encoded state with the final state representing the vector of coefficients with respect to the given function basis.

[0032] In other words, the independent variable may be encoded using a feature map, using known feature map encoding techniques, and an overlap with the ground state representing the vector of coefficients may be determined to obtain the target function value (dependent variable) for the given value of the independent variable.

[0033] In an embodiment, determining the eigenvector of the Hamiltonian with the minimal eigenvalue comprises encoding the Hamiltonian and determining a ground state of the Hamiltonian, the determining the ground state preferably comprising imaginary time evolution of the Hamiltonian. Other embodiments may use other annealing techniques to determine the ground state of the Hamiltonian, e.g., as reviewed in T. Albash et al., 'Adiabatic quantum computing' arXiv:1611.04471v2 [quant-ph] and Y. Suba

et al., 'Quantum algorithms for systems of linear equations inspired by adiabatic quantum computing' arXiv:1805.12549v2 [quant-ph].

[0034] In an embodiment, the basis functions comprise a constant function. This allows easy encoding of generic constraints, such as boundary conditions or data constraints for arbitrary values of the independent variable(s).

[0035] In an embodiment, the basis functions are real-valued functions. In such an embodiment, the differential problem may also be real-valued. Using real-valued basis functions for solving a real-valued differential problem is more efficient, as a solution is sought in a smaller space.

[0036] In an embodiment, the basis functions define a complete function basis. This ensures that derivatives and products of basis functions can be written in terms of basis functions.

[0037] In an embodiment, the basis functions define an orthogonal function basis.

[0038] In an embodiment, the basis functions are the Chebyshev polynomials of the first kind and/or the Chebyshev polynomials of the second kind. The Chebyshev basis functions fulfil all the above-mentioned requirements, and are therefore a natural choice for the basis. Additionally, it is known how to prepare a feature map for the Chebyshev basis, facilitating read-out of the solution. In such an embodiment, the unitary transformation associated with the function basis may be an inverse Chebyshev transformation.

[0039] In other embodiments, a different function basis may be used. The function basis may be selected based on the differential problem. For example, a polynomial function basis is suitable for solving polynomial differential problems.

[0040] In an embodiment, the Hamiltonian is represented by a sum of a first matrix representing the differential equation and one or more second matrices representing the one or more constraints. The first matrix may be defined as a quadratic form of a third matrix, the third matrix being defined as a sum of matrices corresponding to terms of the differential equation with respect to the basis functions. Each of the one or more second matrices may comprise a quadratic form of a fourth matrix, the fourth matrix corresponding to one of the one or more constraints with respect to the basis functions.

[0041] In some embodiments, the second matrices may be provided in pairs, and each pair of second matrices may comprise a quadratic form of the fourth matrix and a product of a normalisation factor with an identity matrix. Such pairs of second matrices may be used to encode inhomogeneous constraints, e.g., a constraint of the form f⁽ⁿ⁾(x₀) = y₀ with y₀ ≠ 0, wherein f⁽ⁿ⁾(x₀) represents an n^th order derivative (with n ≥ 0) of a function f evaluated at a point x₀.

[0042] A quadratic form of a matrix A is defined as the product of a transpose of that matrix with the matrix itself, i.e., A^TA. A quadratic form is always positive (semi)-definite, and hence, all its eigenvalues are larger than (or equal to) zero. A sum of positive (semi)-definite matrices is itself also a positive (semi)-definite matrix, and hence, a Hamiltonian defined this way is guaranteed to be at least positive semi-definite.

[0043] Typically, each term in the differential equation can be mapped to a matrix operation acting on the vector of basis functions.

[0044] In some embodiments, the Hamiltonian may be based on a weighted sum. For example, the constraints may be weighted such that all constraints together have a weight equal to that of the differential equation. In embodiments wherein the constraints are based on measurements, the constraints may have a weight based on an accuracy of the measurement.

[0045] In some embodiments, the first matrix representing the differential equation may project the solution on the null vector. In other words, the differential equation may be represented by a matrix

(the third matrix as defined above), such that

, wherein |ψ_sol〉 represents the eigenvector representing the vector of coefficients of the target function.

[0046] Similarly, each of the fourth matrices representing one of the one or more data constraints may project the solution on the null vector. In other words, each of the one or more boundary conditions may be represented by a matrix

, such that

, wherein |ψ_sol〉 again represents the eigenvector representing the vector of coefficients of the target function. For example, the data constraints may be provided in the form f⁽ⁿ⁾(x₀) = 0, wherein f⁽ⁿ⁾(x₀) represents an n^th order derivative (with n ≥ 0) of a function f evaluated at a point x₀. In other words, in some embodiment, the one or more constraints are homogeneous data constraints.

[0047] In other embodiments, the data constraints may correspond to arbitrary values, e.g., f⁽ⁿ⁾(x₀) = y₀ with typically y₀ ≠ 0. In such cases, it may not be known a priori how many data constraints are needed to obtain the solution; hence, in such cases, additional data constraints may be added until a solution to the differential problem is obtained.

[0048] In an embodiment, the one or more quantum circuits further define operations for encoding a value of an independent variable of the differential equation in a first quantum state using a feature map associated with the basis functions; encoding the eigenvector of the Hamiltonian with the minimal eigenvalue in a second quantum state; and measuring an overlap between the first and second quantum states.

[0049] For example, if the chosen basis functions are the Chebyshev basis functions, the feature map can be the Chebyshev feature map or a derivative Chebyshev feature map.

[0050] In an embodiment, the Hamiltonian comprises a scaling factor. In such an embodiment, determining the eigenvector of the Hamiltonian with the minimal eigenvalue may comprise computing a set of eigenvalues of a part of the Hamiltonian that is scaled by the scaling factor, and selecting a smallest eigenvalue of the set of eigenvalues.

[0051] The determination of the set of eigenvalues may be performed by solving a linear system of equations. This may be done using, e.g., using the quantum processor or using the classical computer system, depending, inter alia, on the size of the problem.

[0052] In an embodiment, the differential equation is a homogeneous differential equation.

[0053] In an embodiment, the differential equation is a polynomial differential equation.

[0054] In an embodiment, applying one or more quantum circuits to the quantum register comprises translating each of the one or more quantum circuits into a sequence of signals (e.g., pulses, typically electric or electromagnetic pulses such as voltage pulses, optical pulses, or microwave pulses) and using the sequence of signals to operate the quantum elements in a quantum register of the quantum processor.

[0055] In an embodiment, receiving hardware measurement data comprises applying a read-out signal to quantum elements of the quantum processor and, in response to the read-out signal, measuring quantum hardware measurement data.

[0056] In an embodiment, the pulses in the sequence of pulses are electromagnetic pulses, e.g., laser pulses or microwave pulses. In such an embodiment, encoding the quantum circuits in the sequence of pulses may comprise determining, for each pulse in the sequence of pulses, one or more of: a pulse duration, a pulse timing, a pulse amplitude, and a pulse phase.

[0057] In general, the embodiments described herein can be implemented as sequences of operations that can be executed on, e.g., a neutral atom quantum computer, either in digital mode of operation, in analog mode, or in digital-analog mode. For example, the Hamiltonian may be encoded in digital mode, whereas the determination of the ground state of the Hamiltonian may be done in analog mode. In particular, a hybrid quantum computer as described herein may include one or more of: a gate-based qubit device, a digital/analog quantum device, a purely analog quantum device, a continuous-variable operating quantum device and/or a gaussian boson sampling device. The hybrid quantum computer may comprise a plurality of intertwined layers, wherein each of the layers is implemented by one of the aforementioned quantum devices. Each of the one or more quantum devices, or each of the plurality of intertwined layers, may be implemented as a neutral atom quantum device, a photonic qubit quantum device, a superconducting qubits quantum device, and/or a trapped ion device.

[0058] One aspect of this disclosure relates to a system for solving a differential problem, the differential problem defining a differential equation and one or more constraints, the solving comprising finding a target function defining an approximate solution of the differential equation, the target function being represented by a vector of coefficients with respect to a given function basis, the system comprising a hybrid computer system comprising a classical computer system and a quantum processor, the hybrid quantum computer system being configured for:

receiving or determining information on one or more quantum circuits encoding a Hamiltonian representing the differential problem with respect to the basis functions; and

determining the solution to the differential problem, the determining comprising:

encoding the quantum circuits using the quantum processor; determining an eigenvector of the Hamiltonian with a minimal eigenvalue, preferably the eigenvalue being equal to zero, the eigenvector representing the vector of coefficients of the target function; and

optionally, measuring the eigenvector to obtain the solution to the differential problem.

[0059] In an embodiment, the hybrid computer system is further configured for executing any one or more of the method steps described above.

[0060] One aspect of this disclosure relates to a hybrid quantum computer system comprising a computer-readable storage medium having computer-readable program code stored thereon, and a processor, preferably a microprocessor, coupled to the computer-readable storage medium, wherein responsive to executing the computer-readable program code, the processor is configured to perform any of the methods described herein.

[0061] One aspect of this disclosure relates to a computer program or suite of computer programs comprising instructions (e.g., at least one software code portion) or a computer program product comprising instructions, which instructions, when executed by a hybrid quantum computer system, configure the hybrid quantum computer to execute any of the methods described herein.

[0062] One aspect of this disclosure relates to a computer-readable medium, e.g., a non-transitory computer-readable storage medium or a computer-readable carrier signal, comprising instructions (e.g., at least one software code portion), which instructions, when executed by a hybrid quantum computer system, configure the hybrid quantum computer to execute any of the methods described herein.

[0063] As will be appreciated by one skilled in the art, aspects of the present invention may be embodied as a system, method or computer program product. Accordingly, aspects of the present embodiments may take the form of an entirely hardware embodiment, an entirely software embodiment (including firmware, resident software, micro-code, etc.) or an embodiment combining software and hardware aspects that may all generally be referred to herein as a "circuit," "module" or "system". Functions described in this disclosure may be implemented as an algorithm executed by a microprocessor of a computer. Furthermore, aspects of the present embodiments may take the form of a computer program product embodied in one or more computer readable medium(s) having computer readable program code embodied, e.g., stored, thereon.

[0064] Any combination of one or more computer readable medium(s) may be utilized. The computer readable medium may be a computer readable signal medium or a computer readable storage medium. A computer readable storage medium may be, for example, but not limited to, an electronic, magnetic, optical, electromagnetic, infrared, or semiconductor system, apparatus, or device, or any suitable combination of the foregoing. More specific examples (a non- exhaustive list) of the computer readable storage medium would include the following: an electrical connection having one or more wires, a portable computer diskette, a hard disk, a random access memory (RAM), a read-only memory (ROM), an erasable programmable read-only memory (EPROM or Flash memory), an optical fibre, a portable compact disc read-only memory (CD-ROM), an optical storage device, a magnetic storage device, or any suitable combination of the foregoing. In the context of this document, a computer readable storage medium may be any tangible medium that can contain, or store a program for use by or in connection with an instruction execution system, apparatus, or device.

[0065] A computer readable signal medium may include a propagated data signal with computer readable program code embodied therein, for example, in baseband or as part of a carrier wave. Such a propagated signal may take any of a variety of forms, including, but not limited to, electro-magnetic, optical, or any suitable combination thereof. A computer readable signal medium may be any computer readable medium that is not a computer readable storage medium and that can communicate, propagate, or transport a program for use by or in connection with an instruction execution system, apparatus, or device.

[0066] Program code embodied on a computer readable medium may be transmitted using any appropriate medium, including but not limited to wireless, wireline, optical fibre, cable, RF, etc., or any suitable combination of the foregoing. Computer program code for carrying out operations for aspects of the present embodiments may be written in any combination of one or more programming languages, including a functional or an object oriented programming language such as Java, Scala, C++, Python or the like and conventional procedural programming languages, such as the "C" programming language or similar programming languages. The program code may execute entirely on the user's computer, partly on the user's computer, as a stand-alone software package, partly on the user's computer and partly on a remote computer, or entirely on the remote computer, server or virtualized server. In the latter scenario, the remote computer may be connected to the user's computer through any type of network, including a local area network (LAN) or a wide area network (WAN), or the connection may be made to an external computer (for example, through the Internet using an Internet Service Provider).

[0067] Aspects of the present embodiments are described below with reference to flowchart illustrations and/or block diagrams of methods, apparatus (systems), and computer program products according to embodiments of the embodiments. It will be understood that each block of the flowchart illustrations and/or block diagrams, and combinations of blocks in the flowchart illustrations and/or block diagrams, can be implemented by computer program instructions.

[0068] These computer program instructions may be provided to a processor, in particular a microprocessor or central processing unit (CPU), or graphics processing unit (GPU), of a general purpose computer, special purpose computer, or other programmable data processing apparatus to produce a machine, such that the instructions, which execute via the processor of the computer, other programmable data processing apparatus, or other devices create means for implementing the functions/acts specified in the flowchart and/or block diagram block or blocks.

[0069] These computer program instructions may also be stored in a computer readable medium that can direct a computer, other programmable data processing apparatus, or other devices to function in a particular manner, such that the instructions stored in the computer readable medium produce an article of manufacture including instructions which implement the function/act specified in the flowchart and/or block diagram block or blocks.

[0070] The computer program instructions may also be loaded onto a computer, other programmable data processing apparatus, or other devices to cause a series of operational steps to be performed on the computer, other programmable apparatus or other devices to produce a computer implemented process such that the instructions which execute on the computer or other programmable apparatus provide processes for implementing the functions/acts specified in the flowchart and/or block diagram block or blocks.

[0071] The flowchart and block diagrams in the figures illustrate the architecture, functionality, and operation of possible implementations of systems, methods and computer program products according to various embodiments of the present embodiments. In this regard, each block in the flowchart or block diagrams may represent a module, segment, or portion of code, which comprises one or more executable instructions for implementing the specified logical function(s). It should also be noted that, in some alternative implementations, the functions noted in the blocks may occur out of the order noted in the figures. For example, two blocks shown in succession may, in fact, be executed substantially concurrently, or the blocks may sometimes be executed in the reverse order, depending upon the functionality involved. It will also be noted that each block of the block diagrams and/or flowchart illustrations, and combinations of blocks in the block diagrams and/or flowchart illustrations, can be implemented by special purpose hardware-based systems that perform the specified functions or acts, or combinations of special purpose hardware and computer instructions.

[0072] The embodiments will be further illustrated with reference to the attached drawings, which schematically will show embodiments according to the invention. It will be understood that the invention is not in any way restricted to these specific embodiments. Identical reference signs refer to identical, or at least similar elements.

Brief description of the figures

[0073] 

Fig. 1 is a schematic illustration of a method according to an embodiment;

Fig. 2A is a flowchart of a method according to an embodiment and Fig. 2B is a flowchart of a corresponding example;

Fig. 3A is a quantum circuit for measuring an overlap, and Fig. 3B-D are graphs showing results of a second-order differential equation solved with a method according to an embodiment;

Fig. 4A-D are graphs showing results of differential equations solved with a method according to an embodiment;

Fig. 5A-D are graphs showing results of a second-order differential equation solved with a method according to an embodiment;

Fig. 6A and 6B are quantum circuits for measuring an overlap, and Fig. 6C-E are graphs showing results of a differential equation solved with a method according to an embodiment;

Fig. 7A and 7B are quantum circuits for measuring an overlap, Fig. 7C and 7D show results of solving nonlinear differential equations using a method according to an embodiment, and Fig. 7E is quantum circuit for encoding a data constraint;

FIG. 8A-C are quantum circuits for quantum imaginary time evolution for use in a method according to an embodiment;

Fig. 9A shows a quantum circuit for a derivative operator, Fig. 9B shows a quantum circuit of a R-shift operator, and Fig. 9C shows an alternative quantum circuit for a derivative operator;

Fig. 10A-D depict circuit diagrams for evaluating an overlap, respectively a derivative thereof, according to various embodiments;

Fig. 11A-C schematically show quantum circuits for an amplitude Chebyshev feature map, a Fourier feature map, and a Chebyshev-grid-based spectral feature map, respectively;

Fig. 12A and 12B schematically show quantum circuits for a quantum Fourier transform and a quantum Chebyshev transform, respectively;

Fig. 13A-E schematically depict hardware-level schematics illustrating the application of logical operations to qubits using a quantum circuit;

Fig. 14A and 14B depict systems describing a Digital-Analog implementation of a quantum algorithm and a quantum feature map, respectively; Fig. 14C illustrates pulses corresponding to the quantum feature map;

Fig. 15A and 15B present quantum circuits based on a Digital-Analog implementation of quantum circuits;

Fig. 16A and 16B present quantum circuits based on a Digital-Analog implementation of quantum circuits; and

Fig. 17A and 17B are block diagrams illustrating, respectively, an exemplary hybrid data processing system and an exemplary classical data processing system that may be used for executing methods and software products described in this disclosure.

Detailed description

[0074] This disclosure describes quantum solvers for differential problems (comprising a differential equation and one or more data constraints) that combine an effective Hamiltonian approach with spectral methods to address various types of differential equation in a latent space. Execution of such a quantum solver by a hybrid quantum computer results in an encoding of a solution to the differential problem (up to a normalisation factor) in a quantum register of the hybrid quantum computer. This encoded solution may be used to, e.g., evaluate and/or plot the solution; more typically, however, the encoded solution is used as a starting point for further computational steps using the hybrid quantum computer. For instance, the encoded solution may be used to generate samples according to a distribution implicitly defined by the differential problem. Thus, the method could be considered a method for encoding a solution to a differential problem in a quantum register. A further advantage of the method is that, as the solution is encoded in a ground state, the solution is relatively stable over time.

[0075] The differential equation can be a linear or a nonlinear differential equation, and can be an ordinary or partial differential equation. The differential equation may be a homogenous differential equation. The differential equation may be a polynomial differential equation, or a polynomial approximation of a differential equation. A polynomial differential equation is a differential equation containing only polynomial terms of the dependent and independent variables.

[0076] The data constraints can be, e.g., boundary conditions, initial conditions, final conditions, (internal) fixed points, et cetera. In some embodiments, the data constraints may be based on measurements or observations on a system that is being modelled.

[0077] The function basis that is used for the spectral analysis can any suitable function basis, e.g., a Chebyshev basis comprising Chebyshev polynomials of the first and/or second kind. The formulation of quantum models in terms of the Chebyshev basis allows for the latent-space construction of an effective Hamiltonian operator of the underlying differential equations accompanied by given data constraints within a physics-informed SciML paradigm. In other words, it recasts the differential equations into a total energy function minimization problem, which can be mapped on an eigenvalue problem. Consequently, this method does not require variational optimization (although methods using variational optimization are not excluded), but enables finding a ground state (i.e., the zero-energy eigenstate) solution finding through algebraically solving a system of equations.

[0078] The choice of the latent-space basis or function basis to construct the quantum model directly affects the formulation of the effective Hamiltonian. If a linearly-independent and purely-real function basis with a purely-real differentiation matrix is used, such as the aforementioned Chebyshev basis, both linear and nonlinear differential equation can be transformed into an eigenvalue problem of a real-valued effective Hamiltonian in a finite-dimensional Hilbert space. This eigenvalue problem may be described by a concise form of the total energy function (as explained below with reference to, e.g., eq. (23), (38), and (46)).

[0079] By contrast, when a Fourier basis is utilized, the solution (overlap model) is formulated as a Fourier-series, and the Fourier differentiation matrix is purely imaginary. Due to the complex-valued quantum state, the resulting total energy function is of a complicated form, rendering it difficult to find a ground state in this scenario. Similarly, if Hartley basis is selected, the top-left entry of Hartley differentiation matrix is not a constant, being a differential operator, and thus it does not allow us to easily access to the derivative of the overlap model. Therefore, in the following example, the quantum models will be built based on the Chebyshev expansion. In addition, the x-dependent Chebyshev basis state can be readily prepared from available Chebyshev feature map circuit.

[0080] In principle, the ground state may be prepared using any known method, including fully analogue systems and digital solutions such as quantum imaginary-time evolution algorithms based on quantum singular value transformation (QSVT). The solution can be evaluated by determining an overlap between the ground state (representing the solution in latent space) and a feature map associated with the function basis. E.g., if the Chebyshev basis is used, a quantum Chebyshev feature map circuit may be used to encode independent variables. The Chebyshev feature map circuit is described in WO 2025/002636 A1 and C.A. Williams et al., 'Quantum Chebyshev Transform: mapping, embedding, learning and sampling distributions', arXiv:2306.17026 [quant-ph], which are hereby incorporated by reference in its entirety.

[0081] In the figures, a few examples are provided showcasing various types of differential equations. Nevertheless, it is understood that the current method is not limited to these examples.

[0082] Fig. 1 is a schematic illustration of a method according to an embodiment. In this particular example, a general solution to a second-order differential equation (DE) is sought,

. A quantum model f_Q(x) of an (independent) scalar variable x is constructed as a sum of three weighted Chebyshev polynomials of the first kind, T_k(x) ≡ cos(k arccos(x)), to approximately represent the unknown solution function f(x) by a polynomial of degree 2 (for the purpose of schematic illustration shown in Fig. 1, panel (a)). In other words, the solution has the form

, or more specifically, a quantum model is constructed such that

, wherein 〈τ_k(_x)| represents the Chebyshev basis,

is a scale factor, and |ψ〉 represents a vector of coefficients with respect to a given function basis, in this case the Chebyshev basis.

[0083] With this quantum model, each term in the differential equation is associated with a quantum state. For instance, f(x) ≃ f_Q(x) is mapped to a scaled quantum state, represented as a column vector of three undetermined coefficients, and scaled by a factor of

in a latent (Chebyshev) space. The derivative term

is mapped to a transformed quantum state in the same latent space. As will be shown below, a matrix

may be constructed that maps a function in the Chebyshev basis to its derivative, so that



. Higher-order derivatives can be obtained by repeated application of such matrix. Thus, the example differential equation may be encoded as

, wherein

represents the identity matrix:

[0084] More in general, every term of the differential equation, and each constraint of the differential problem (e.g., boundary/initial conditions) has its own latent representation.

[0085] Advantageously, the constraints are chosen such that the function or a derivative thereof is invariant under scaling, which is typically achieved by the function value or derivative being equal to zero (possibly after a function transformation). In the present example, as shown in panel (b), these data constraints may be given by:

and/or

[0086] By combining (e.g., summing) these individual latent representations, an effective Hamiltonian operator

may be constructed such that the solution of the differential problem is a ground state of the Hamiltonian, e.g.,



, or

. An energy function ε(η) yielding a scalar energy can then be defined as an expectation value of

with respect to the (so far) unknown quantum state |ψ〉. This can be considered the equivalent representation of the differential problem in the latent space. A more detailed description of the construction of the effective Hamiltonian

is provided below with reference to Figs. 2 and 3.

[0087] The minimization of the energy function ε(η), and thus the solving of the differential problem, can be treated as a machine-learning task by either varying three independent parameters {c_k} for k = 0, 1, 2, related to each component of |ψ〉. Alternatively, the eigenproperties of the effective Hamiltonian operator can be explored. The effective Hamiltonian operator may be constructed as a function of a single variable η > 0. An advantage of the latter approach is that it does not require a variational Ansatz, thus avoiding common training issues (such as barren plateaus and local minima trapping) that are generally applicable to all variational protocols.

[0088] In general, the effective Hamiltonian operator

can be represented by a real symmetric matrix with distinct real eigenvalues and mutually orthonormal eigenvectors. In the eigen-space representation, ε(η) can be expressed as a sum of three squared projection length

weighted by the eigenvalue λ_k associated with eigenvector q_k. It is known that (i) the diagonal entries of Λ control the definiteness of ε(η) and that (ii) the total energy is always non-negative, and by construction equal to zero for the solution to the differential problem. These facts imply that

is a positive-semidefinite matrix with exactly one eigenvalue equal to 0 for a certain value η ∈ {η_i} = argmin_η ε(η) with i = 0, 1, 2. In particular, a positive-semidefinite

with one eigenvalue being zero is achieved when η = η (corresponding to λ = λ_min), where η_max is the largest value in {η_i}. Thus, the energy minimization leads to the solution |ψ_sol〉 (up to normalisation), which is the lowest-energy eigenstate of the (effective) Hamiltonian

. The normalisation factor or scaling factor

may be obtained from a non-zero data constraint (in this example, boundary condition f(-1) = 1).

[0089] Thus, in an exemplary embodiment, the method may comprise the steps of: formulating the differential equation and, optionally, the data constraints as overlap models in a quantum Chebyshev space; summing the latent contributions to form an effective Hamiltonian; preparing a ground state of the effective Hamiltonian, e.g., by thermalization and/or QSVT; and applying a Chebyshev feature map to evaluate the overlap model for a given input value.

[0090] Fig. 2A is a flow chart of a method according to an embodiment and Fig. 2B provides an explicit example. In a first step 202, a differential problem is obtained (e.g., received or generated) as input. The differential problem comprises a differential equation DE and zero or more data constraints DC, e.g., boundary conditions, an initial value, measurement data, et cetera. Typically, the number of data constraints should be equal to or larger than the number of degrees of freedom of the differential equation minus one, but the number may be (substantially) higher, e.g., when fitting a function to measurements.

[0091] In the current example, the differential equation reads

[0092] Two different kinds of data constraints are considered: (i) a homogeneous data constraint

[0093] and (ii), a set of s inhomogeneous data constraints:

[0094] This differential problem can be solved analytically, and its solution is (for the homogeneous case)

with C a scaling factor.

[0095] A step 204 comprises selecting a function basis. The function basis may be used to express the solution as a vector of coefficients with respect to that basis. The function basis may be selected based on the differential problem. In some cases, it may be necessary to rescale the independent variable, as is known in the art. A useful function basis, especially for polynomial, real-valued differential problems, is the Chebyshev basis. In general, a system with n encoding qubits gives access to N = 2ⁿ basis functions.

[0096] A Chebyshev polynomial of the first kind of order k may be defined as:

[0097] A Chebyshev polynomial of the second kind of order k may be defined as:

[0098] Depending on the computational problem being solved, Chebyshev polynomial of the first or second kind, or a combination thereof, may be best suited to encode the solution function. The corresponding state |x〉 is

[0099] The solution is then obtained as a vector of coefficients with respect to the selected function basis.

[0100] The problem can then be recast into finding a vector |ψ〉 = [c₀, c₁, ... , c_N-1]^T such that

. is a solution to the differential problem. This vector |ψ〉 is x-independent.

[0101] For the current example, it is assumed a two-qubit system is used, and that the function basis comprises the zeroth through third Chebyshev polynomials of the first kind:

[0102] Based on the selected model, a step 206 comprises writing each term of the differential equation as a matrix operator in a latent space.

[0103] The derivative of a polynomial is a lower-order polynomial; hence, when a polynomial basis is used, the derivative can be written on the same polynomial basis.

Higher-order derivatives can be obtained by repeatedly applying

:

[0104] And hence, the differential equation can be written as:

[0105] With

, where

is the identity matrix. As will be shown in more detail below, matrix operations can also be defined for several other operations, such as multiplication by x^p; as this increase the maximum order of the polynomial basis, one or more additional qubits may be used to increase the set of available basis functions.

[0106] As all x-dependency is contained in |τ(x)〉, and the differential equation must be true for all x, the differential equation DE can be written as:

[0107] This equation is automatically obeyed when |ψ〉 is an eigenvector of

with eigenvalue zero.

[0108] For the data constraints (both the homogeneous and inhomogeneous ones), the function at independent variable value x_j is given by:

[0109] (using T₀(x) = 1) and with

[0110] this can be written as

[0111] where in the homogeneous case, y₀ = 0. From this, it follows that:

using 〈ψ|ψ〉 = 1 in the last step.

[0112] A step 208 comprises constructing an effective Hamiltonian

based on the matrix operators. The effective Hamiltonian may be parametrised by a scaling factor η. In some embodiments, the effective Hamiltonian may comprise distinct weights for the contribution of the differential equation and the data constraint(s). However, in the current example, all weights are identical.

[0113] Thus, taking the above matrix operators together (in this example, with equal weights), leads to:

[0114] for the homogenous case and to:

[0115] for the inhomogeneous case. In either case, this can be written as:

[0116] Which is solved by a zero-eigenvalue eigenstate of effective Hamiltonian

.

[0117] It is noted that by defining a reduced Hamiltonian

as

[0118] the inhomogeneous problem can be redefined as finding the (unique) eigenvector |ψ_sol〉 with eigenvalue zero of the effective Hamiltonian

. This is equivalent to finding a normalised eigenvector v_sol with eigenvalue λ_sol of W such that |ψ_sol〉 = v_sol and

.

[0119] A subsequent step (not shown) comprises receiving or determining one or more quantum circuits defining operations of the effective Hamiltonian

.

[0120] A step 210 comprises determining an eigenvector of the effective Hamiltonian

using a quantum processor, e.g., by determining the ground state of the effective Hamiltonian. This is discussed in more detail below. The solution in latent space is then given by

, and the general solution by

. If necessary, the scaling factor

may be determined by evaluating the function for a known point (e.g., a non-zero data constraint).

[0121] An optional step 212 comprises evaluating the determined solution, using the quantum model solution, at one or more input points. In particular, a known feature map may be used to encode |x_j〉; for example, |τ(x_j)〉 =

(x_j)|Ø〉 wherein

(x) is the Chebyshev feature map. Thus,

[0122] Additionally or alternatively, an optional step 214 comprises generating samples using the encoded solution, e.g., by applying a quantum transform associated with the used feature map (e.g., a quantum Chebyshev transform for the current example) and performing a projective measurement.

[0123] Similarly, a function evaluation of a derivative of the target function can be obtained by computing an overlap with a derivative Chebyshev feature map.

[0124] As noted before, once the solution has been encoded as a ground state of the effective Hamiltonian, it may be used for various kinds of computations, not limited to mere function evaluations. For example, the solution may be used to generate random samples that follow a distribution encoded by the solution, as input for a quantum machine learning algorithm (e.g., a quantum neural network, QNN) such as a classification algorithm, and so on.

[0125] It is noted that in some cases, the ground state of the effective Hamiltonian may not correspond to a solution of the differential problem. For example, the data constraints essentially appear in squared form, so that sign information may be lost in the process. In such cases, the method may further comprise adding additional data constraints (and/or replacing the data constraints already used), and repeating steps 206-210.

[0126] More in general, the following observations may be made. Although in the following, reference is made to the Chebyshev basis using Chebyshev polynomials of the first kind, the same or similar considerations apply to other function bases. Depending on the specific context, the function basis may fulfil one or more of the following requirements: real-valued, comprising a constant function, complete, orthogonal, and/or orthonormal.

[0127] Similarly, in the following x is used as a general independent variable, which can be either one- or multidimensional. Some specificities related to multidimensional variables (and/or several one-dimensional variables) are discussed in more detail further below.

[0128] Since the basis functions T_k(x) satisfy a discrete orthogonality condition and form a complete orthogonal basis over the interval [-1, 1], any arbitrary function that is smooth and continuous on the same interval can be expressed as

, with c_k =

for k > 1 and

. This is known as a Chebyshev expansion. In practical applications, a good approximation can generally be obtained using only a finite set of k < N basis functions. To prevent overfitting, it may be beneficial if N is not too large. In practice, it has been found that for many applications, good results are obtained if N is chosen as, e.g., 16, 32, or 64, which generally requires, respectively, 4, 5, or 6 qubits for encoding the coefficients; thus, an economical use may be made of the available computational resources. Functions or problems that are defined on a different (preferably compact) domain, can be scaled to the interval [-1,1].

[0129] Similar to the Chebyshev expansion, the quantum model f_Q(x) is built with finite n qubits as a state overlap, represented by the scaled Hermitian inner product of an x-dependent Chebyshev basis state |τ(x)〉_n and a quantum state |ψ〉_n (generally consisting of unknown coefficients to be determined), living in N = 2ⁿ bra and ket spaces, respectively, given by

[0130] where

denotes a scaling factor, η > 0. The quantum model may also be referred to as an overlap model, as a function evaluation typically involves an overlap of two quantum states, one quantum state representing the function coefficients (x-independent) and the other quantum state representing the function basis for a given value of the independent variable (x-dependent).

[0131] The purely real Chebyshev basis state reads

[0132] where |k〉, k = 0, ...,N - 1 is the n-qubit computational basis. The Chebyshev basis states in eq. (5) are orthonormal when evaluated at the so-called Chebyshev nodes x_j := cos(π(j + ½)/N) ∈ (-1, 1), defined as the roots or zeros of T_N(x), and also when evaluated at the Chebyshev extrema

. These Chebyshev states

or

form an orthonormal basis set.

[0133] Eq. (4) indicates that the unknown |ψ〉_n is naturally separable from the x-dependent basis, i.e., all x-dependency of f_Q(x) is contained in the function basis and the coefficients vector |ψ〉 is x-independent. Therefore, the quantum model f_Q(x) can be symbolically differentiated with respect to x via a Chebyshev differentiation matrix,

, as follows:

[0134] where the fact that

is employed and

is a constant upper triangular matrix with T in the superscript denoting a matrix transpose. An example of a quantum circuit for implementing derivative operator

is shown in Fig. 9. The double-struck capital letter (also called blackboard bold) is utilized to signify a non-unitary matrix throughout the application. This constant differentiation operator allows easy access to the m^th-order derivative of the quantum model by matrix multiplication of

m times, i.e.:

where

denotes the m^th derivative of f_Q(x).

[0135] The Chebyshev polynomials and their derivatives are related by:

where U_k(x) denotes a Chebyshev polynomial of the second kind of order k. This can then be further expanded to:





This allows the derivation of the matrix elements for

. This also means that the derivative of a Chebyshev polynomial is a linear combination of lower order Chebyshev polynomials. Therefore, an effective operator

can be extracted such that:

(where the subscript n has been dropped).

[0136] If one wants to evaluate a derivative of the solution function, the operator

may be used to encode a derivate Chebyshev feature map. The operator

is non-unitary, but can be implemented using known techniques to implement non-unitary operators, such as a linear combination of unitaries or embedding in a higher-order space.

[0137] In addition, using the product-to-sum identity of Chebyshev polynomials, it can be shown that the product of a variable x to the p^th power and 〈τ(x)|_n equates to the product of 〈τ(x)|_n+1 and an x-independent matrix

:

with the integer p ∈ [0, N]. This requires an additional qubit; i.e., 〈τ(x)|_n denotes an n-qubit basis, and 〈τ(x)|_n denotes an (n + 1)-qubit basis. If higher powers are required (p > N), more additional qubits may be used; each additional qubit doubles the set of basis functions. When necessary, differential equations that include a non-polynomial function of the independent variable may be approximated using, e.g., a Taylor series expansion or other suitable polynomial expansion of the non-polynomial function, or, in some cases, they can be converted into a polynomial differential equation using a suitable coordinate transformation.

[0138] It is generally desirable to encode all elements of the differential problem (e.g., terms of the differential equation, and data constraints) with respect to the same basis. Thus, if the (n + 1)-qubit basis 〈τ(x)|_n is used for some elements, other elements of the differential problem may be defined with respect to the same basis. To that end, projections matrices

and

may be defined, such that:

[0139] These non-square projection/embedding matrices are unique constant matrices that allow for the projection of the n-qubit quantum state to a higher dimensional Chebyshev basis (or vice versa). Such projection is useful because , 〈τ(x)|_n and 〈τ(x)|_n+1 are linearly-independent basis sets, while both x^p 〈τ(x)|_n and (〈τ(x)|_n ⊗ 〈τ(x)|_n) are not linearly-independent basis sets.

[0140] The data constraints can be included in the same formalism by formulating them as state overlaps with the basis functions, to obtain their latent representations. Usually, the term "boundary conditions" is used for constraints imposed at one or more boundaries of the problem domain (e.g., x₀ = -1 and x₁ = 1), whereas the more general term "data constraints" also covers internal constraints. Such constraints can be written as, for example:



where

and

are obtained from

and 〈τ(x₁)|_n =

, respectively. This leverages the fact that the zeroth-order Chebyshev polynomial is just a constant function, and consequently,



where

is an x-dependent (data-related) matrix with nonzero elements only for the first row.

[0141] It is useful to distinguish between "invariant" or "homogeneous" constraints, and "regular" or "inhomogeneous" constraints. The former are constraints of the form

(in which case r_x0 and s_x1 as used in eqs. (16) and (17), respectively, are equal to 0), whereas the latter are of the form

(and hence r_x0 ≠ 0 and s_x1 ≠ 0). Here, f^(m)(x) represents the m^th derivative of f(x).

[0142] In some cases, a regular constraint may be converted into an invariant constraint by a suitable transformation, e.g.,

et cetera. In general, regular data constraints can be used to determine a normalisation factor or scale of the (inherently normalised) solution in the latent space.

[0143] Additionally, regular data constraints can be used to encode polynomials of the independent variable, as can be seen from the following equations:

and

so that

Non-polynomial functions may be approximated using a polynomial expansion, e.g., using a Maclaurin series:

In the following,

will be used for notational simplicity.

[0144] Thus, both a function value f(x₀) and a derivative function value f'(x₁) can be expressed as their respective state overlaps in the same basis as well and the corresponding weighting factors r_x0 and s_x1 are then computed, subject to 〈c_r|_n|c_r〉_n = 〈c_s|_n|c_s〉_n = 1. |c_r〉_n = |c_s〉_n = |Ø〉 always holds in the case of constant boundary conditions (f(x₀) = y₀ and

. Eqs. (16) and (17) imply that available data can be supplied in the form of boundary conditions and serve as the physics-informed constraints, playing a crucial role to guide the model toward the desired outcome. Unlike boundary/initial conditions, data constraints are not necessarily restricted to be located at the boundaries (e.g., x₂ and x₃ ∈ [-1,1]).

[0145] Therefore, BC = {f(x₀), f'(x₁)} may be regarded as regular boundary conditions used to determine a unique analytical solution to a physical problem and DC = {f(x₂), f'(x₃)} as given data-related constraints used to minimize the total energy, which make two situations distinguishable by symbols. Note that the elements of DC can be exactly the same as those of BC, partially inclusive of boundary terms (DC = {f(x₀), f(x₂), f'(x₃)}) or exclusive when the boundary values are unavailable in reality. Both boundary and initial conditions can be incorporated into the effective Hamiltonian as part of data constraints to capture the local behaviour of the solution to DEs.

[0146] As

is a sparse matrix with non-zero elements in its first row, it can be easily implemented by block encoding techniques such as Fast Approximate Block-Encoding algorithm (FABLE). With FABLE, the implementation of unitary

, a block encoding of

requires 2n + 1 qubits (including n + 1 ancilla qubits), which typically contains 2n Hadamard gates and n SWAP gates in general and 2²ⁿ R_Y and 2²ⁿ CNOT gates at most. The same holds for (right-side) products of

, such as

.

[0147] Fig. 3A-D depict a solution of a differential problem comprising an ordinary differential equations with constant coefficients using a method according to an embodiment.

[0148] In particular, Fig. 3A schematically shows a circuit for evaluation of a function value for a given input value x. The circuit is initialized with a zero-product state |Ø〉 = |0〉^⊗n and a state preparation circuit

is responsible for preparing the desired state of a problem of interest, followed by an inverse Chebyshev feature map

to obtain the overlap probabilities in the Chebyshev basis. In this example, a generic preparatory circuit

is used to prepare the ground state of the effective Hamiltonian, represented by |0_a〉 |ψ〉_n. Some examples of such preparatory circuits (based, e.g., on quantum imaginary-time evolution) are discussed in more detail below with reference to Fig. 8A-C. This preparatory circuit is followed by an inverse Chebyshev feature map

to perform an overlap measurement with respect to the continuous variable x, given by 〈0_aτ(x)|0_aψ〉.

[0149] It is stressed that the function evaluation is given by the state overlaps themselves, and not by the absolute value square of the state overlaps. Note that f_Q(x) = Re[f_Q(x)] as the resulting quantum states are purely real under Chebyshev basis

, where the subscript a refers to the ancilla qubit(s) and |Ø〉 ≡ |0〉^⊗n denotes an n-qubit zero product state. In general, such overlap measurement can be carried out through the Hadamard test, or any other suitable method, as discussed in more detail below with reference to Fig. 10A-D. However, the Hadamard test relies on the use of a set of controlled operations via an extra ancillary qubit to control each quantum gate in the circuit. As a result, the implementation of the Hadamard test for evaluations of overlaps tends to increase circuit complexity and decrease noise resiliency. By contrast, interferometric measurements, which involve both ground and combined states can be implemented more easily on current hardware. It turns out that a single overlap measurement can be obtained from two separate evaluations of the expectation value of an observable

= |0_aØ〉 〈0_aØ| using the same circuit architecture; this will be discussed in more detail below with reference to Fig. 6A and 7A, and eqs. (24), (41), and (54).

[0150] As an example, a second-order ordinary differential equation (ODE) is considered with constant coefficients (a, b and

):

[0151] The state-overlap form of this ODE is

where

is an x-independent matrix. Because 〈τ(x)|_n is a linearly independent basis set, the x-dependent part can be removed (as the equation has to be valid for all x), keeping only the latent representation of the DE, that is,

. This leads to

by multiplying its transpose from the left, i.e., computing the quadratic form. By the same token, according to eqs. (16) and (17), the quadratic forms of given data constraints DC = {f(x₀ = 0), f'(x₁ = 0)} are

and

, respectively. Next, the latent differential equation may be combined with the data constraints to obtain a total energy function ε(η) as a sum of individual latent contributions given by

where

is the effective Hamiltonian matrix that consists of the symmetric (quadratic) forms of the latent representation of the x-independent differential equation and the x-dependent (data-related) constraints. In general, different weights may be used for different terms in the Hamiltonian, but for simplicity's sake, all weights are taken equal to 1.

[0152] Since eq. (23) is the latent representation of eq. (22), the differential problem becomes an eigenvalue problem

such that

. Considering the total energy is always non-negative, this means that if one of the positive eigenvalues of

goes to zero, the corresponding eigenvector is the ground state of the solution to eq. (22). Simply put,

must be a positive-semidefinite matrix for certain η ∈

with the integer i ∈ [0, N - 1].

[0153] In the case that the data constraints are inhomogeneous, the effective Hamiltonian can instead be written as:

, wherein Î_n is the n-qubit identity matrix, and wherein

and

are defined in eqs. (16) and (17). In that case, the requirement of the zero eigenvalue (λ = 0) enables the determination of the scaling factor

. Given that w_i is the (i + 1)^th positive eigenvalue of the symmetric positive-definite matrix

can be specifically found by solving a set of linear equations

. The positive-semidefinite

with (exactly) one eigenvalue being zero can be accomplished when

, where

is the largest value in the set

, if appropriate datasets are supplied. Since the smallest eigenvalue w_min is the eigenvalue closest to zero, the selection of

can take the minimum energy consumption to reach the ground state of the overall system based on their reciprocal relation.

[0154] In Fig. 3B-D, three possible types of solution to eq. (22) are investigated, with coefficients C = {a, b, c} and data constraints BC = {f(0), f'(0)}. The overlap model f_Q(x) is evaluated under different data constraints DC_I. Plots of f(x) and f_Q(x) are given for

3B C = {1, 4, 4}, BC = {0.5, -0.5} and DC_I = { f(-1) = 0}, η_e = 1.29.

3C C = {1, -2, -3}, BC = {0, 1} and DC_I = { f(0) = 0}, η_e = 32.47.

3D C = {1, 5,400}, BC = {1,0} and DC_I = {f'(0) = 0}, η_e = 451.714.

[0155] In each panel, n and η_e represent the number of qubits and the scaling factor of the quantum model, respectively.

[0156] Fig. 3B and 3C, show exemplary differential equations having real repeated and distinct roots, whose analytical solutions are f(x) = 0.5 (exp(-2x) + x exp(-2x)) and f(x) = 0.25 (exp(3x) - exp(-x)), respectively; this represents a lognormal-like function (dashed curve). For the supplied data constraint, the lowest-energy eigenstate of



is denoted as |ψ_g〉_n and the resulting state overlap solution

(gray solid curve) is in direct proportion to the exact solution. Here and in the following, T denotes the Gram operator:

.

[0157] For complex pairs of roots, as in the case of Fig. 3D, the differential equation has an analytical solution

, which represents a damped harmonic oscillator. Both η_e and ground state |ψ_g〉_n are then calculated according to the aforementioned methods. The resultant overlap model is displayed as the solid curve in Fig. 3D.

[0158] Similar results were obtained in simulations using inhomogeneous data constraints (not shown). As data constraints can manipulate the matrix elements of W and control the definiteness of

, the selection of inhomogeneous data constraints can determine the eigen-behaviour of the overlap model. Adequate inhomogeneous data constraints can help the model to rapidly probe the ground state of the problem.

[0159] Eq. (24) or (25) can be used to evaluate the overlap models of ordinary differential equations and non-linear differential equations, and to verify the solutions. This is also known as an interferometric measurement,

where

is the overlap probabilities of measuring the Chebyshev basis state |0_a〉τ(x)〉_n given (n + 1)-qubit ground state |ψ_G〉_n+1 = |0_a〉|ψ_g〉_n where |ψ_g〉_n is prepared by the imaginary-time evolution based on the obtained

, and

is the overlap probabilities of measuring the Chebyshev basis state |0_a〉|τ(x))_n given the (n + 1)-qubit combined state |ψ_C〉_n+1 = (|ψ_G〉_n+1 + |ψ_r〉_n+1)/

where

= ∥ψ_G〉_n+1 + |ψ_r〉_n+1|² is a normalization factor. The reference state |ψ_r〉_n+1 is set as |0〉^⊗(n+1) so as to easily decouple the reference function when computing the real part of the product term, and N = 2ⁿ is used. In the special cases of either f(x) ≥ 0 (for instance Fig. 3B) or f(x) ≤ 0 ∀x ∈ [-1,1], it is sufficient to evaluate the |Ô_G(x)|² term.

[0160] Fig. 4A-D show results of solving Legendre differential equations and associated Legendre differential equations using a method according to an embodiment. These are examples of ordinary differential equations with variable coefficients.

[0161] As a further example, the general Legendre differential equation with 0 ≤ m ≤ l ∈ N is discussed, which is given by:

[0162] The complete solution to this ordinary differential equation is

where

and

are two linearly independent functions called associated Legendre functions of the first and second kind, respectively. They are frequently used when solving Laplace's equation and scattering problem in spherical coordinates. For the azimuthally symmetric case (m = 0), the state-overlap form of this differential equation is

where

is an x-independent matrix. The multiplication of x² with the n-qubit Chebevshev basis indicates that the leading power of x is raised from N - 1 to N - 1 + 2 in the Chebyshev polynomials, so that n + 1 qubits are needed to ensure that each term in the equation shares the same linearly independent basis set (see eqs. (13) and (14); in eq. (27), the subscript n has been dropped); consequently, each term in A comprises a matrix

, with 0 ≤ p < N (using

for brevity). This allows to remove the x-dependent part and keep intact the latent representation of the ODE.

[0163] Data constraints including function and derivative values are also expressed in terms of the (n + 1)-qubit Chebyshev basis, and thus eqs. (16) and (17) are modified as follows:



where

and

are obtained from

and

, respectively. They are x-dependent (data-related) matrices with nonzero elements on their first row as usual, and may be encoded using sparese matrix encoding techniques as described above. With l-dependent boundary conditions, this equation has an analytical solution

, the Legendre polynomials (shown in black dashed curves). The data constraints are chosen such that

and

The total energy function ε(η) is the same as eq. (23).

[0164] In the case of m ≠ 0, multiplying both sides of eq. (26) by (1 - x²) yields the latent governing equation

In the case of m = 1 with fixed boundary conditions, this differential equation has an analytical solution

, the associated Legendre polynomials, again shown as black dashed curves. The effective Hamiltonian has effectively the same shape as for m = 0, but with a different matrix A as described in eq. (30). Following the same procedures as described above, the evaluated quantum models (shown as solid curves) closely follow the target functions for m = 0 and exhibit best-fit behaviour even though singular points occur at x = -1 and x = 1 for m = 1.

[0165] More specifically, results are shown for the following simulations. In each panel, n represents the number of qubits of the quantum model. As before, η_e represents the scaling factor.

[0166] For m = 0, eq. (26) with l-dependent boundary conditions BC = {f(-1) = (-1)^l, f(1) = 1} is used to determine a unique analytical solution

. The overlap model f_Q(x) is evaluated for different integers of l under the same data constraints DC = {f(0), f(±1)}. The plots show

and f_Q(x) for (Fig. 4A) l = 4 and η_e = 1.75 and (Fig. 4B) l = 5 and η_e = 1.49.

[0167] For m = 1, eq. (26) with fixed boundary conditions BC = {f(-1) = 0, f(1) = 0} is used to determine a unique analytical solution

. The overlap model f_Q(x) is evaluated for different integers of l under the different data constraints DC. The plots show

and f_Q(x) for (Fig. 4C) l = 5 and DC = { f'(0) = 0}, η_e = 108.51 and (Fig. 4D) l = 6 and DC = {f(0) = 0}, η_e = 135.38.

[0168] Again, similar results were obtained (not shown) using inhomogeneous data constraints (for m = 1, l = 5 and DC = { f(±x₀), f(±x₁₂)}, and for m = 1, l = 5 and DC = {f(±x₀), f(±x₈)}). Since, as noted, eq. (26) has two regular singular points at -1 and 1 for m ≠ 0, the 0^th and j^th Chebyshev nodes of opposite signs located between two singular points are chosen as DC. Finally, we have

. The right choice of appropriate and meaningful datasets is essential in today's data-driven model to facilitate physics-informed machine learning in order to unravel insights and discover patterns as well as to extract valuable knowledge and underlying information. It is found that the overlap models also in these cases closely follow the target functions for m = 0, and exhibit small relative errors nearby the singular points of -1 and 1 for m = 1.

[0169] Fig. 5A-D show results of solving second-order inhomogeneous differential equations using a method according to an embodiment. The second-order inhomogeneous differential equations are given by

with coefficients and a source function on either side of the equal sign, which may be denoted as

= {o(x), p(x), q(x);

.

[0170] Fig. 5A-D shows solution for eq. (31) with coefficients and a source function denoted as

= {o(x), p(x), q(x); r(x)} and data constraints BC = {f(0), f'(0)} (which are used to determine a unique analytical solution f(x)), as follows:

Fig. 5A:

= {(x - 1), -x, 1; (x - 1)²}, BC = {1/2 , -1/8} and DC_I = {f'(-0.2) = 0}, η_e = 1.35.

Fig. 5B: C = {1, 4, 4; e^-2x}, BC = {-1, 1} and DC_I = {f(-0.73) = 0}, η_e = 15.74.

Fig. 5C:

= {1, -5, 6; x e^2x}, BC = {1/3, 1/3} and DC_I = {f(0.111) = 0}, η_e = 0.82.

Fig. 5D: C = {1, 4, 20; xe^-2x sin 4x}, BC = {-1, -1} and DC_I = {f'(-0.75) = 0} , η_e = 77.51.

Here, η_e is again the scaling factor used to obtain the scaled solution. Plots are shown for f(x) and

. In each panel, n represents the number of qubits of the quantum model.

[0171] For

= {(x - 1), -x, 1; (x - 1)²} and data constraints DC = {f(-0.2) = 0, f(0) = 0.5}, i.e., the parameters used for Fig. 5A, this differential equation eq. (31) has an analytical solution f(x) = 3/2 e^x - x (x + 13/8) - 1. Following the same procedures as described above, eq. (21) is used to find the latent representation of r(x). The quantum model of eq. (31) is

, where

is a data-encoded matrix, where again

. With the supplied homogeneous data constraints, the lowest-energy eigenstate of



is calculated, and then the state overlap solution is found. The scaling factor can be determined by f(x_s = 0) = 0.5, the same as the one appearing in A. Eventually, the evaluated quantum model agrees with the exact solution, as shown in Fig. 5A. For complicated (e.g., non-polynomial) source functions, a polynomial expansion such as the Maclaurin expansion of r(x) is applied and the first t + 1 terms are taken for approximation eq. (21). For simplicity, eq. (22) with r(x) on the right hand side may be considered. The corresponding latent governing equation is

Again, three exemplary differential equations with real repeated, real distinct, and complex pairs of roots are solved (shown in Fig. 5B-D, respectively), whose analytical solutions are

and f(x) =

, respectively. Given the abovementioned data constraints, the effective Hamiltonian operators are

for the first two cases and

for the last. The resultant quantum models evaluated with t = 3, 5, 7 reveal excellent agreement with the exact solutions, as shown in Fig. 5B-D, respectively.

[0172] Fig. 6A and 6B show quantum circuits and Fig. 6C-E show results of solving partial differential equations using a method according to an embodiment. In particular, Fig. 6A shows a quantum circuit that can be used to evaluate the two-dimensional state overlaps, 〈0_a τ(x) 0_a τ(y)|0_a ψ_x 0_a ψ_y〉, through an interferometric measurement where a single overlap measurement can be obtained from two separate measurements of overlap probabilities. The circuit is initialized with a zero-product state and a state preparation circuit

is employed to prepare the desired state of a problem of interest. A permutation circuit P̂ is responsible to shuffle the ancillary and prepared states in a specific order and generate separable input states for parallel inverse Chebyshev feature maps to obtain the overlap probabilities in the Chebyshev basis. An alternative circuit, not requiring a permutation subcircuit, is discussed in more detail below, around eq. (49).

[0173] Fig. 6B shows an alternative quantum circuit, where the effect of the permutation circuit has been absorbed in the state preparation circuit.

[0174] Fig. 6C-E show, respectively, (Top) intensity-scale surface plots of quantum effective Hamiltonian model-solved (Fig. 6C) Laplace's equation (eq. (33)) with n = 3 and η_e = 5.21559, (Fig. 6D) heat equation (eq. (39)) with n = 4 and η_e = 603.863 and (Fig. 6E) wave equation (eq. (40)) with n = 5 and

. (Bottom) the corresponding error plots ε defined as the difference between the estimated and true values. All plots share the same intensity scale bar on the right.

[0175] The methods described herein can be applied to both univariate and multivariate cases. The effective Hamiltonian allows to express the solution of differential problem as the ground state of a multidimensional system, without the need of a quantum variational Ansatz and costly grid evaluations, as used in known methods. For simplicity, we first demonstrate how to tackle a partial differential equation (PDE) involving partial derivatives of an unknown function (or a dependent variable) f(x, y) with respect to two independent variables x and y. For encoding of two independent variables (as shown in Fig. 6A and 6B), the quantum model f_Q(x, y) can be built in parallel registers of equal qubit numbers n_x = n_y = n as the two-dimensional state overlap given by

[0176] In the following, 〈τ(x, y)|_2n is used as shorthand notation for the two-dimensional linearly-independent Chebyshev basis set 〈τ(x)|_n ⊗ 〈τ(y)|_n. As an example, Laplace's equation (which is an elliptic PDE that arises in many applications, such as electrostatics and fluid dynamics, to describe steady-state phenomena), is solved s=using the following parameters:



Eq. (33) has an analytical solution f(x, y) = csch(π) cosh(πx/2) sinh((y + 1) π/2). The state-overlap form of this PDE is

, where

is an x- and y-independent matrix with the symbol ⊕ denoting the Kronecker sum operation. Boundary and initial conditions including function and derivative values are formulated as the two-dimensional state overlap as well. Homogeneous Dirichlet and Neumann boundary conditions respectively take the form

[0177] Similar to eqs. (16) and (17),

and

are obtained from 〈τ(y₀)|_n =

and

, respectively. They are y-dependent (data-related) matrices with nonzero elements on their first row as usual.

[0178] In contrast to 2D collocation points in PINNs, the regions of minimum variation of the function value and/or the first derivative are chosen as homogeneous data constraints, DC = {f(±1,y) = 0,f(x, -1) = 0}, and/or

, such that

with at least two weights w_i ≠ 0; or more in general, at least one non-zero weight for each independent variable.

[0179] Again, inhomogeneous data constraints may be represented by:



wherein u_y0 and v_y1 are the corresponding weighting factors associated with the function value and the derivative value, subject to 〈c_u|_n|c_u〉_n = 〈c_v|_2n|c_v〉_2n = 1 and the same for the other two.

[0180] Also in that case, the regions of minimum variation of the function value and the first derivative are chosen as data constraints,

, such that



is obtained, where the (symmetric) Gram operator

(·) = (·)^T(·) is used to concisely represent the data constraints (operands) in the formula.

[0181] In either case (homogeneous or inhomogeneous data constraints), the total energy function for PDE is written as

Both η_e and ground state |ψ_g〉_2n are then calculated according to the aforementioned methods; typically by minimizing eq. (38) to obtain the lowest-energy eigenstate (ground state) of the effective Hamiltonian

. The exact value of the scaling factor can be obtained via

using an inhomogeneous data constraint {f(x_s, y_s) ≠ 0} . The quantum model can be measured with the quantum circuit shown in Fig. 6A and eq. (43).The resultant overlap model is displayed as a grey-level surface plot in Fig. 6C and nicely follows the hyperbolic sine and cosine function with small errors nearby its maximum along the boundary.

[0182] As a further example, a time-dependent PDEs in one spatial dimension is shown in Fig. 6D and 6E. The dependent variable is a function of time t and (one-dimensional) space x. Two representative examples, the heat and wave equations, are examined as these parabolic and hyperbolic PDEs are widely used to describe many physical and engineering problems such as diffusion processes/heat flow and vibrating systems/wave motion, respectively.



where IBC is the set of initial and boundary conditions (i.e., data constraints) to specify the distribution at the initial time and domain boundaries to pin down an unique solution. Eqs. (39) and (40) have unique analytical solutions f(t, x) = exp(-4 π²k t) cos(2πx) and f(t, x) = cos(4πt) sin(2πx), respectively. The corresponding quantum model is f_Q(t, x) =

revised from eq. (32), with 〈τ(t, x)|_2n representing the two-dimensional linearly-independent Chebyshev basis set 〈τ(t)|_n ⊗ 〈τ(x)|_n. Both PDEs have the same state-overlap form

, with their respective



and

. The independent variables x and y in eqs. (34) and (35) are respectively replaced by t and x. Given

and DC = { f(1/8,x) = 0, f(7/8,x) = 0, f(t, ±1) = 0}, this ultimate leads to effective Hamiltonians

for eq. (39) and

for eq. (40). The results are shown as intensity-scale surface plots in Fig. 6D and 6E, respectively. Overall, the quantum models closely follow mixed decaying and oscillating (heat) solution and highly oscillatory (wave) solution, with small errors appearing at their maxima and minima. Again, very similar results may also be obtained using (a larger number of) inhomogeneous data constraints (results not shown).

[0183] Eq. (41) can be used to evaluate the overlap models of partial differential equations, and to verify the solutions.

where

is the overlap probability of measuring the two-dimensional Chebyshev basis state |0_a〉τ(x)〉_n|0_a〉τ(y)〉_n given (2n + 2)-qubit ground state |ψ_G〉_2n+2 = |0_a〉|ψ_gx〉_n|0_a〉|ψ_gy〉_n = P̂(|0_a〉|0_a〉|ψ_g〉_2n) where |ψ_g〉_2n is prepared by the imaginary-time evolution based on the obtained

and P̂ is a permutation circuit responsible to shuffle the ancillary and prepared states in a specific order and generate separable input states for parallel registers, and |Ô_C(x, y)|² =

is the overlap probability of measuring the two-dimensional Chebyshev basis state |0_a〉|τ(x)〉_n|0_a〉τ(y)〉_n given (2n + 2)-qubit combined state |Ψ_C〉_2n+2 = (|ψ_G〉_2n+2 + |ψ_r〉_2n+2)/

where

= ||ψ_G〉_2n+2 + |ψ_r〉_2n+2|² is a normalization factor. The reference state |ψ_r〉_2n+2 is set as |0〉^⊕(2n+2) so as to easily decouple the reference function when computing the real part of the product term, and N = 2ⁿ is used.

[0184] Fig. 7A and 7B show quantum circuits and Fig. 7C and 7D show results of solving nonlinear differential equations using a method according to an embodiment. In particular, Fig. 7A shows a quantum circuit that can be used to evaluate the state overlaps, 〈0_a τ(x)|0_a ψ〉, through the interferometric measurement where a single overlap measurement can be obtained from two separate measurements of overlap probabilities. The circuit acts on two (n + 1)-qubit registers, each qubit register comprising an ancilla qubit and n data qubits. The circuit is initialized with a zero-product state and a state preparation circuit

is responsible to prepare the desired state of a problem of interest. This is followed by a permutation circuit P̂. A data constraint circuit D̂ is applied to the first (n + 1)-qubit register, followed by a layer of Hadamards H. An inverse Chebyshev feature map is then applied on the second (n + 1)-qubit register to obtain the overlap probabilities in the Chebyshev basis.

[0185] Fig. 7B shows an alternative quantum circuit, where the effect of the permutation circuit has been absorbed in the state preparation circuit.

[0186] Fig. 7C shows a plot of f(x) and f_Q(x) for eq. (42) with DC = {f'(0) = 0,}, η_e = 3.75367. Fig. 7D shows a plot of f(x) and f_Q(x) for eq. (47) and DC = {f(0.02615) = 0}, η_e = 0.7984. In both panels, n represents the number of qubits of the quantum model.

[0187] In this final example, the developed strategy is applied to nonlinear differential equations (NDEs). We first consider a nonlinear differential equation whose analytical solution is an even function, 1 - x²/8, written as

The problem is solved following the same workflow as in the previously described linear cases, while noticing that the quantum state is encoded twice as |ψ〉_n ⊗ |ψ〉_n in the nonlinear (function squared) term while the quantum state for the linear term is only encoded once. To maintain the same dimension of Hilbert space, the linear terms are multiplied with the identity (e.g., as expressed in eq. (19)) such that |ψ〉_n ⊗ |ψ〉_n appears in all linear terms (at the expense of a 'doubled' Hilbert space).

[0188] Thus, for non-linear differential equations, the quantum model of eq. (4) and the homogeneous boundary conditions of eqs. (28) and (29) become, after introducing the identity:





where eq. (15) is employed to create a linearly independent basis set. For conciseness, (〈ψ|_n ⊗ 〈ψ|_n) may also be written as 〈Ψ|_2n (and analogously for analogous expressions). The equations contain f(x_s) encoded in the expressions. It is noted that eq. (43) is essentially identical to eq. (4) and both equations are applied to the linear and nonlinear terms of non-linear differential equations, respectively. The quantum model, represented by the state-overlap form of the nonlinear differential equations can be written as

where

is an x-independent matrix. Given the data constraint DC = {f'(0) = 0}, the effective Hamiltonian can be written as

. The total energy function ε(η) for the non-linear differential equation is then given by:

[0189] Like linear DEs, solving the second-order NDEs requires one invariant constraint to obtain the state overlap solution and one regular constraint to determine the scaling factor. The regular constraint is selected as DC_R = {f(x_s = 0) = 1}, the same as the one appearing in

.

[0190] The consequence of the doubled Hilbert space is that

possesses at least 2^2n-1 degenerate states associated with zero eigenvalue λ_deg = 0. Most importantly, the task of solving an NDE is converted to that of seeking a degenerate state |Ψ_g〉_2n = |ψ_g〉_n ⊗ |ψ_g〉_n that closely matches the analytical solution, which can be verified by substituting the chosen ground state into eq. (43) and determining the scaling factor with the given DC_R. On the other hand, the quantum model in eq. (43)indicates that the same independent variable x is encoded to both registers. In this case, based on the identical quantum model in eq. (4), one can simply replace the inverse Chebyshev feature map in the first register with a layer of Hadamards on each qubit and then measure the overlap in the Hadamard-Chebyshev basis, as illustrated in Fig. 7A. The quantum model is measured with the described quantum circuit along with eq. (54), and the result is consistent with the analytically exact solution, as shown in Fig. 7B.

[0191] Another example shows a non-linear differential equation with a quadratic non-linearity on the unknown function, given by

[0192] This non-linear differential equation has no analytical solution, so instead a numerical solution to the non-linear differential equation is found. The latent governing equation is

Given DC_I = {f(0.02615) = 0}, the resultant effective Hamiltonian operator is



. DC_R = {f(x_s = 0.5) = 0} is utilized and the choice of DC_R only affects the accuracy of the solution.

[0193] After performing measurements, the result is shown in Fig. 7D, demonstrating excellent agreement with the numerical solution. Overall, the representation of quantum states in a higher dimensional space is the key to solve NDEs in the degenerate eigenspace of the effective Hamiltonian

for the methods described herein.

[0194] In general, boundary (initial) conditions specify the value of the dependent variable at the domain borders (initial time) of the independent variable, which impose constraints on the solution to restrict the solution space and guarantee a unique solution to the differential equation. In PINNs (or quantum PINNs), the optimized solution is found by (variationally) minimizing a total loss function, e.g., a mean squared error (MSE), defined as the sum of individual contributions that match the differential equation and satisfy the boundary/initial conditions and data-based constraints. In general, PINNs requires a costly loss function evaluation at uniformly-spaced collocation points (

×

) for solving PDEs like those shown in Fig. 6D and 6E. By contrast, the methods described herein do not require evaluation at a plurality of collocation points.

[0195] In the methods described herein, however, the aim is to minimize the total (global) energy function, the expectation value of an effective Hamiltonian

, which is the sum of the symmetric forms of latent PDE and data constraint matrices whose size is determined by the number of qubits (n_t + n_x) used to encode the temporal and spatial quantum states.

[0196] Different from PINNs (or quantum PINNs), the formulated energy function is regarded as an eigenvalue problem and thus the solution naturally corresponds to the lowest energy eigenstate that resides in the eigenspace of

. Furthermore,

is a real symmetric matrix, which guarantees that the eigensolution exists unlike the variational or quantum matrix inversion methods. This also bypasses variational optimization and avoids barren plateaus and trainability issues in initialization and throughout optimization, offering a novel means to reach the solution of DEs in a variation-free manner. The method can thus be regarded as a novel means to reach the solution of DEs in a machine-learning amenable and yet variation-free manner. Notwithstanding the natural target of (early) fault-tolerant quantum devices, the present algorithm may be run on current, noisy quantum hardware leveraging the plethora of methods existent in the literature to (at least approximately) solve eigenvalue problems. This is demonstrated for the key step in solving the Legendre DE example (Eq. 17), using a realistic emulation of neutral atom QPUs to prepare the relevant ground states, using the pulser library. This library is compatible with cloud-accessible devices, and allows to design a pulse sequence preparing the targeted states with high-fidelity, while incorporating hardware constraints like the modulation bandwidth of the pulse.

[0197] In addition, two data constraints suffice to solve the second-order linear/nonlinear DEs (as shown in Figs. 3-5 and 7) and three to four data lines (Fig. 6) or discrete points along lines for PDE problems. In general, for higher-order linear and nonlinear DEs, d - 1 invariant constraints are needed and one regular constraint, where d is the number of degrees of freedom (e.g., the order of the differential equation). Consequently, the methods described herein do not require costly evaluations on dense collocation points, which typically raises computational costs when the problem dimensionality increases, making the current approach distinct from PINNs.

[0198] The invariant constraint plays a crucial role in separating the overlap model and scaling factor, enabling us to independently tackle them. It enforces a physics-informed restriction on the null space of

to ensure the state overlap solution respecting the laws of physics governed by the differential equation, while meeting the specifics of data and improving the model toward the desired outcome. Simple rescaling of a state overlap solution is more efficient than directly solving a differential equation with extra cost for a large number of collocation points. It is worth mentioning that

is positive-definite for linear differential equations, and thus the ground state corresponds to the lowest-energy state. However, for non-linear differential equations,

is (only) positive semi-definite, so the ground state can be one of a number of degenerate states associated with zero energy.

[0199] Alternatively, several non-homogeneous (or regular) data constraints may be used; e.g., for the shown examples, it is enough to solve linear/nonlinear DEs with 3-5 data points (Figs. 4, 5, and 7) and PDEs with 5-7 data lines (Fig. 6) or with discrete points along lines. In that case, the matrix elements of

may be deliberately manipulated through supplying appropriate data constraints such that the resultant

turns out to be a positive-semidefinite matrix. Note that x-independent

itself is a positive-semidefinite matrix and acts as a global (regularization) term to ensure the solution respects the laws of physics governed by the differential equation, while x (or t)-dependent

is a rank-one matrix and behaves as a local term responsible to meet the specifics of the solution and improve the model toward the desired outcome. Unlike PINNs (or quantum PINNs), the formulated total energy function can be regarded as an eigenvalue problem and thus the solution naturally corresponds to the scaled zero-energy eigenstates that resides in the null-space of

. This bypasses the variational optimization and avoids the barren plateaus and trainability issues in the initialization and throughout the optimization, offering a novel mean to reach the solution of DEs in a variation-free manner.

[0200] The former A is composed of qubit-dependent constant matrices

and the later

is dependent on the given local temporal and spatial datasets. These data constraints can be provided at anywhere within the domain rather than limited at the domain extreme ends like boundary/initial conditions.

[0201] An advantage of the (Chebyshev-enabled) Hamiltonian differential eigensolver described herein concerns the capability of systemically building an effective Hamiltonian and expressing the ground state of a multidimensional system. Moreover, the described effective Hamiltonian method offers the flexibility to adapt to the inputs of a multivariate system in a cost-effective manner. For example, the previously included permutation circuit P̂ in Fig. 6A can be removed if eq. (32) is replaced with the equivalent quantum model

where

is a non-square constant matrix that has only one 1 in each row. In this alternative model, the unknown quantum state |ψ〉_2n+1 is implicitly engineered to directly follow the desired input pattern of the inverse Chebyshev feature maps, |0_a〉ψ〉_2n+1 = |0_a〉|ψ_x〉_n|0_a〉|ψ_y〉_n, and therefore P̂ is no longer needed.

[0202] Dirichlet and Neumann boundary conditions (as described in eqs. (34) and (35)) accordingly become

[0203] For homogenous data constraints, r_x0, u_y0, s_x1, and/or v_y0 are equal to zero. For inhomogeneous data constraints, the corresponding weighting factors associated with the function and derivative values are calculated, subject to



.

[0204] For the case of Laplace's equation, the state-overlap quantum model is

with

. The (2n + 1)-qubit ground state can be found by eigen-solving the total energy function

using the same approach as described previously. After discarding N² degenerate zero-energy states, the ground state |ψ_g〉_2n+1 associated with the lowest eigenvalue is found. Both approaches return the same values of η_e under the same data constraints for different PDE problems, but this alternative approach provides the ground state exactly matching the inverse Chebyshev feature map circuits without the need of the complex permutation circuit (because the effect of P̂ has been absorbed into the effective Hamiltonian operator). Therefore, the resulting quantum circuit (Fig. 6B) used to evaluate the two-dimensional state overlaps is similar to that in ODEs (Fig. 3A).

[0205] Fig. 7B shows an alternative quantum circuit for applying the methods described herein to PDEs. Quantum circuit used to evaluate the two-dimensional state overlaps, (0_a τ(x) 0_a τ(y)|0_a ψ_x 0_a ψ_y〉, through the interferometric measurement where a single overlap measurement can be obtained from two separate measurements of overlap probabilities. The circuit is initialized with a zero-product state and a state preparation circuit is utilized to prepare the desired state of a problem of interest for parallel inverse Chebyshev feature maps to obtain the overlap probabilities in the Chebyshev basis.

[0206] To verify if the same procedure is applicable to NDEs, Eqs. (28), (29) and (30) are rewritten as:





where

is a non-square constant matrix with the matrix elements like

. For the case of Eq. (39), the corresponding quantum model is

, with



, and the total energy function is

, with

. The (2n + 1)-qubit ground state can be found from the degenerate zero-energy space of

with the number of degeneracy increased by N². Again, both approaches return the same value of η_e. The resulting quantum circuit, which can be used to measure the state overlap related to the quadratic nonlinearity, is depicted in Fig. 7B.

[0207] It is noted that if the solution function has neither zero crossings nor zero slopes in the Chebyshev domain, invariant constraints do not exist. For instance, consider this situation with a NDE whose analytical solution is a cubic function -(x - 3)³/27, written as:

[0208] Apparently, invariant constraints are not available for this analytical solution in the interval [-1,1]. A new dependent variable may be defined as f̃(x) = f(x) - f(0) = f(x) - 1 so that Eq. (39) becomes:

[0209] Now, eq. (53) has DC_I = {f̃(0) = 0} and the methods described above for homogeneous or invariant data constraints are applicable. The corresponding effective Hamiltonian operator is

with



, where DC_R = { f̃(-1/16) = f(-1/16) - 1} is utilized. As soon as

is obtained, the original quantum model solution can be retrieved by



.

[0210] Eq. (54) allows evaluation of the quantum model that contains the quadratic nonlinearity.

where

is the overlap probability of measuring the Hadamard-Chebyshev basis state |+_a〉|+〉^⊗n|0_a〉|τ(x)〉_n given 2(n + 1)-qubit ground state |Ψ_G〉_2n+2 = (Î₁ ⊗ D̂_n ⊗ Î_1+n)(Î₁ ⊗ P̂)|0_a〉|0_a〉|Ψ_g〉_2n with |Ψ_g〉_2n prepared by, e.g., imaginary-time evolution (as described below) based on the obtained

. P̂ is a permutation circuit and D̂ is a data constraint circuit with

being a subnormalization factor.

is the 2-norm of |Ψ_G〉_2n+2.

is the overlap probability of measuring the Hadamard-Chebyshev basis state |+_a〉|+〉^⊗n|0_a〉|τ(x)〉_n given 2(n + 1)-qubit ground state |Ψ_C〉_2n+2 = (|Ψ_G〉_2n+2 + |Ψ_r〉_2n+2)/

with

being the 2-norm of (|Ψ_G〉_2n+2 + |Ψ_r〉_2n+2). The reference state |Ψ_r〉_2n+2 is set as |0〉^⊗(2n+2) so as to easily decouple the reference function when computing the real part of the product term, and N = 2ⁿ is used.

[0211] Fig. 7E is an example of a data constraint circuit, for example, circuit D̂ in Fig. 7A and 7B. More in particular, Fig. 7E shows an implementation of

via linear combinations of unitary operations, where

.

is an inverse Chebyshev feature map, and Ŝ_n+1 is a special type of the Grover diffusion operator for performing a reflection about the all zeros state.

[0212] Fig 8A-C show various quantum circuits that can be used in embodiments.

[0213] Fig. 8A shows quantum circuits for QSVT-based quantum imaginary time evolution. Panel (a) shows a

(ϕ, x) circuit used to generate a polynomial approximation of a mixed parity function based on the sum of even and odd polynomials via LCU. Relevant QSP subroutine circuits for even and odd parity are shown in Fig. 8B. Panel (b) shows a Re[

(ϕ, x)] circuit used to find the phase angles generating a polynomial approximation to e^-tx, where t is the evolution time. Panel (c) shows a

circuit used to generate a mixed parity polynomial transformation of a normalized effective Hamiltonian operator

given a set of optimal phase angles obtained from (b). Relevant QSVT subroutine circuits for even and odd parity are shown in Fig. 8C. Panel (d) shows a

circuit producing

in top left block used to find the ground state of

in the second n-qubit register after measuring the first (a + 3)-qubit register and postselecting on the |0〉^⊗(a+3) outcome.

is a symmetric operator, in this context typically the effective Hamiltonian. H denotes the Hadamard gate.

[0214] Fig. 8B shows quantum circuits for even and odd parity QSP circuits. Quantum signal processing (QSP) is performed by the repeated applications of the single-qubit reflection and parameterized Z-rotation gates, R̂(x) and Ŝ(ϕ). The phase angles ϕ = {ϕ^e, ϕ^o} comprise even and odd parity groups,

and

where d^e and d^o are the maximum degrees of even and odd polynomial components of p(x), respectively. Note that R̂(x) is a symmetric operator.

[0215] Fig. 8C shows quantum circuits for even and odd parity QSVT circuits. Quantum singular value transformation (QSVT) operation is composed of repeated applications of block-encoding unitary and projector-controlled phase gates,

and Π̂(ϕ), where



is a normalized effective Hamiltonian (symmetric) operator. ϕ_opt =

is a set of optimal phase angles.

is an a-qubit block-encoding of

, which is constructed using the combination of PREPARE and SELECT subroutines based on Pauli decomposition of

.

[0216] As soon as the effective Hamiltonian operator

is ready, the zero-energy eigenstate (ground state) can be prepared by the quantum imaginary-time evolution (QITE). Known subroutines, such as quantum signal processing (QSP) and quantum singular value transformation (QSVT), may be used to realize QITE. First, phase angles are searched for generating a polynomial approximation of the target function, p(x) ≈ e^{-t x}. Since the function e^{-t x} has indefinite parity, p(x) can be synthesized by the sum of even and odd polynomials. This sum operation is accomplished by a mixed parity polynomial approximation circuit

(ϕ, x) consisting of opposite parity QSP circuits

and

via a linear combination of unitaries (LCU), as shown in Fig. 8A, panel (a), where the phase angles ϕ = {ϕ^e, ϕ^o} are split into even and odd parity groups,

and ϕ^o =

for each QSP circuit (Fig. 8B), with d^e and d^o being the maximum degrees of even and odd polynomial components of p(x), respectively. Note that p(x) is a purely real polynomial corresponding to the top left entry of Re[

(ϕ, x)] circuit which can be similarly constructed using a linear combination of unitaries, as shown in Fig. 8A, panel (b). The phase angles are then adjusted until p(x) gets closer to the target function e^{-t x} and then a set of optimal phase angles

is obtained. We then upload the optimal phase angles to a mixed parity polynomial transformation circuit

composed of opposite parity QSVT circuits

and

via LCU, as shown in Fig. 8A, panel (c). The argument

is a normalized effective Hamiltonian operator,

=

. Inside each QSVT circuit (Fig. 8C), the unitary

is an a-qubit block-encoding of

, which can be constructed using PREPARE and SELECT subroutines based on Pauli decomposition of

, with the numbers of ancillary qubits proportional to the numbers of Pauli string operators. Finally,

is a purely real matrix polynomial corresponding to an (a + 3)-qubit encoding of

circuit which can be similarly constructed via LCU, as shown in Fig. 8A, panel (d). The acceptable value of time evolution t is determined by the maximum eigenvalue of

. For the DEs shown in Fig. 4A and 4B, t = 120 is enough to recover the n-qubit ground state in the second register.

[0217] Fig. 9A shows a quantum circuit for a derivative operator and Fig. 9B shows a quantum circuit of a R-shift operator. More in particular, Fig. 9A shows a quantum circuit

for a derivative operator

for a Chebyshev basis for 3 qubits (n = 3). The quantum circuit uses a first register 902 with n = 3 qubits (shown at the bottom) for encoding the Chebyshev basis functions, a single-qubit temporary register 904 (shown in the middle) to temporarily store qubit values for the R-shift circuits in block 930, a second register 906 with m = n - 1 = 2 qubits, and an ancilla 908. The first register may comprise a non-zero input, while the temporary register, the second register, and the ancilla are initialised in the zero state. The qubits in the first register may be numbered 1 ... n from top to bottom, with qubit n considered the least significant bit. For the sake of controlling the ancilla (or other qubits), the qubits in the first register may be interpreted as a bit string; i.e., if a qubit is controlled by bitstring 101 = 5, it is controlled by the first and third qubits, but not by the second.

[0218] The second register 906 is equipped with Hadamard gates H at the beginning and end of the quantum circuit. Two initial R_Y gates are applied to the ancilla 908. Subsequently, a first block 910 comprises n R_Y gates acting on the ancilla, shown as R_Y(θ_k), k = 1,...,n, which are controlled by the (n + 1 - k)^th qubit (or bitstring 2^k-1) in the first register 902. A second block comprises 2^m - 1 R_Y gates acting on the ancilla, shown as R_Y(ϕ_l), l = 1,..., 2^m - 1, which are controlled by the bitstring 2l + 1 in the in the first register 902 and the bitstring l in the second register 906.

[0219] The argument variables are as follows:

, θ_k = -2^k-1Δ, and φ_l = 2 arccos[α(2l + 1)] - π + (2l + 1) Δ, with

(N - 1), and N = 2ⁿ.

[0220] The third block 930 of the quantum circuit is composed of a series of R-shift circuits R̂_n, an example of which is shown in Fig. 9B.

[0221] Fig. 9C shows an alternative quantum circuit for a derivative operator, more in particular a Chebyshev differentiation circuit. The implementation of

via linear combinations of unitary operations, where

, is a normalised differentiation operator. The m-register is equipped with m Hadamard gates at the beginning and end of the circuit (where m = n). The R_Y gates acting on the top ancilla are controlled by the bit string |j〉_n = |2^k-1〉 in the n register of

, and bit strings |i〉_m = |0〉^⊗m and |j〉_n = |2^k-1〉 in the m and n registers of

, respectively. The argument variables are as follows: ω_c =

, θ_k = 2 arcsin(-2 ^k-1 α),

. The last part of each unitary subroutine is composed of n R-shift (R̂_n) circuits for

and one L-shift (L̂_n) circuit plus 2 n SWAP gates for

. Fig. 9C shows the case for n = 3 is shown here.

[0222] Fig. 10A-D depict circuit diagrams for evaluating an overlap, respectively a derivative thereof, according to various embodiments. Other embodiments may use different circuits to measure the overlap of two wave functions. In particular, Fig. 10A shows a circuit diagram for implementation of a squared overlap of the general form

[0223] In this example, a so-called SWAP-test (as defined in H. Buhrman et al. "Quantum fingerprinting", Physical Review Letters 107:16 (2001), page 167902, which is hereby incorporated by reference) is used to measure the expectation value |〈ψ_j|ψ_k〉|² of an overlap of two wave functions |ψ_j〉 and |ψ_k〉 directly. The wave functions |ψ_j〉 and |ψ_k〉 typically correspond to the solution wavefunction |ψ_sol〉 (representing function coefficients) and the basis functions of the chosen basis (e.g., a one-dimensional or multi-dimensional Chebyshev basis). The swap test may be performed by having two parallel registers of n qubits, one for a first wave function |ψ_j〉 1002 and one for a second wave function |ψ_k〉 1004, plus a single ancillary qubit 1006. Thus, the swap test requires 2n + 1 qubits where n is the number of qubits used to express |ψ_j〉.

[0224] A unitary operator

is applied to the registers for |ψ_j〉 and |ψ_k〉, respectively, in order to obtain |ψ_j〉 =

|0〉 1008 and |ψ_k〉 =

|0〉 1010. A Hadamard gate H 1012 is applied to the ancillary qubit, and then a controlled-SWAP operation 1014 is performed between all qubits of the other registers, conditioned on the ancilla qubit with a control operation 1016. The controlled swap onto the size 2n register is made up of controlled swap on the i^th qubit of each of the two size n registers, for i ∈ 1: n. Subsequently, another Hadamard gate 1018 is applied to the ancillary qubit and the qubit is then measured in the Z basis 1020 to read the squared overlap value |〈ψ_j|ψ_k〉|².

[0225] If ψ_j and ψ_k are equal, the SWAP operation will not do anything, and hence the two Hadamard gates will combine to an X gate and the measurement results always exactly 1. If ψ_j and ψ_k are different, the controlled-SWAP gate will incur a phase-kickback on the ancillary qubit, which is detected as an average value of the measurement between 0 and 1. This way, the squared overlap can be estimated using this circuit, which is physically implementable even on current-generation qubit-based quantum computers.

[0226] Fig. 10B shows a circuit diagram for an alternative implementation of an overlap measurement. This circuit may be used to determine an overlap

(or functions thereof). Using this circuit, the real and imaginary parts of the overlap 〈ψ_j|ψ_k〉 may be measured by means of a so-called 'Hadamard-test'. In this case, both x and y are mapped to one register of n qubits 1022. A single ancillary qubit 1024 may be used for measurements. Thus, the swap test requires n + 1 qubits where n is the number of qubits used to express |ψ_j〉.

[0227] First a Hadamard gate 1028 is applied to the ancillary qubit. Then, an (undaggered) unitary operator

is applied to the register to create a quantum state |ψ_k〉 =

|0〉, followed by a daggered unitary operator

, in order to obtain a quantum state |ψ_j,k〉 =

|0〉 1030. The application of the unitary operators is controlled by the ancillary qubit 1032. This circuit employs the identity

[0228] An S-gate S^b = exp(-i b π Z/4) is then applied 1034 to the ancillary qubit, where Z is the Pauli-Z operator. It is noted that if b = 0, this gate does nothing. This is followed by an Hadamard gate 1036 and then a measurement operation 1038. Steps 1034-1038 may be applied twice, once with b = 0 to measure the real part of the overlap Re(〈ψ_j|ψ_k〉) and once with b = 1 to measure the imaginary part of the overlap Im(〈ψ_j|ψ_k〉), but in many applications, only either the real or imaginary part is needed, so that a single application suffices.

[0229] Based on these measurements, various overlap-related functions may be constructed (typically using a classical computer), e.g.

or

[0230] Fig. 10C shows a quantum circuit diagram for an alternative implementation of an overlap measurement. The quantum circuit is similar to Fig. 10B, but instead of the controlled-unitary 1030, a control-free version 1042 of the (same) unitary is used. Finally, all n qubits are measured (rather than only the ancillary qubit), and the number of times a zero is observed divided by the total number of measurements is an estimate of the squared overlap value |〈ψ_j|ψ_k〉|².

[0231] Fig. 10D is a circuit diagrams of Hadamard tests to measure overlaps. Fig. 10D shows a Hadamard test measuring

and

for b = 0 and b = 1, respectively. S^b denotes the phase gate, exp(-π b Z/4). When b = 0 and b = 1 are used the real and imaginary part is evaluated, respectively. By summing over b = 0/1 the full overlap can be evaluated.

[0232] The construction of the circuit shown in Fig. 11A, which encodes a data point x with a Chebyshev feature map, is based on the observation that the Chebyshev polynomials are defined as T_k(x) = cos(k arccos(x)) and are therefore of the form of cosine; in particular,

. The cosine can be recovered from

, which can be implemented with the use of an ancillary qubit.

[0233] Thus, the Chebyshev polynomials can be prepared using a combination of exponents for some scaled variable x, where each amplitude is embedded via a phase feature map, as described in O. Kyriienko et al., 'Protocols for trainable and differentiable quantum generative modelling', Phys. Rev. Research 6 (2024) 033291, which is hereby incorporated by reference in its entirety. The positive and negative phase can be combined using a so-called Linear Combination of Unitaries (LCU) approach, as described in Childs et al., 'Quantum algorithm for systems of linear equations with exponentially improved dependence on precision', SIAM Journal on Computing 46 (2017) pages 1920-1950, which is hereby incorporated by reference in its entirety.

[0234] In this approach, maps are conditioned on the state of the (top) ancilla qubit, and effectively interfered, choosing the ancilla collapsed to 0 outcome.

[0235] In a first step 1102, a layer of Hadamard gates is applied to all qubits (including the ancilla) to prepare a superposition of all basis states. In a step 1104, a single-qubit scaled phase shift gate

is applied to the encoding qubits, which is defined as

where j is the qubit number. This is followed by a layer 1106 of controlled single-qubit scaled phase feature maps

that are conditioned on the ancilla. Conceptually, the first layer corresponds to exp(i x) and the second layer corresponds to exp(-i 2x). However, since exp(-i x) = exp(i x) exp(-i 2x), this results in equal-weight combinations of exponents exp(i x) + exp(-i x). Thus, the amplitudes for the basis functions are scaled, based on the variable x such that T_k(x) are recovered.

[0236] Finally, the weight of the constant (T₀(x)) term is adjusted. The depicted example uses a round of Grover's iterate circuit 1108, rotating the state around |Ø〉 for a fixed angle -π/2. As the measurement operation 1112 on the ancilla (including the preceding Hadamard gate) commutes with rotation, the measurement can be postponed to the end of the circuit. This allows post-selection of the desired state also on systems where only all qubits can be read out simultaneously.

[0237] This way, the resulting quantum Chebyshev feature map creates a Chebyshev state via an x-parametrized isometry.

[0238] The layers 1102-1108 may together be referred to as the Chebyshev feature map 1110, as in practice the measurement on the ancilla may often be postponed, with one or more other circuits in between. The Chebyshev feature map may be represented by a unitary operator

(x).

[0239] This results in a state such that at the set of the N = 2ⁿ Chebyshev nodes an orthonormal basis is formed. Therefore, there exists a bijection and corresponding transformation between this basis and any different orthonormal basis. Fig. 5A shows the Chebyshev transform which performs the transformation between the Chebyshev and computational bases.

[0240] Fig. 11B schematically shows a quantum circuits for a Fourier feature map (Fourier encoding or spectral encoding), for the sake of comparison. This circuit does not use an ancilla. The circuit comprises a layer of Hadamard gates 1122 to create a superposition of all basis states, followed by a layer 1124 of (unscaled) single-qubit phase feature maps

[0241] Fig. 11C schematically shows a quantum circuits for a phase-encoded Chebyshev-like feature map, for the sake of comparison. This circuit does not use an ancilla. The circuit comprises a layer of Hadamard gates 1132 to create a superposition of all basis states, followed by a layer 1134 of (scaled) single-qubit phase feature maps

[0242] For the sake of clarity, it is noted that the feature maps referred to as Chebyshev encoding feature maps in previous disclosures by the present inventors are fundamentally different from the Chebyshev feature map described herein. For example, the previously described feature maps have, in general, a non-zero imaginary part, and are not orthonormal. This prevents their conversion to amplitude encoding (which is standard in quantum information processing), and building distributions for generative modelling in particular. By contrast, the currently disclosed Chebyshev feature map is orthonormal and real (i.e., the imaginary part is zero), and may therefore be considered an amplitude encoding.

[0243] Fig. 12A schematically shows an example of a quantum Fourier transform circuit, as may be used in layer 1206 of the quantum Chebyshev transform shown in Fig. 12B.

[0244] Fig. 12B schematically shows an example of a Quantum Chebyshev transform circuit according to an embodiment. This circuit maps computational basis states

into Chebyshev states

. Similar to the Chebyshev amplitude encoding, this circuit uses an ancillary qubit. In the depicted example, an Hadamard gate is applied to the ancilla, followed by a CNOT-ladder controlled by the ancilla.

[0245] The transform involves a quantum Fourier transform (QFT) circuit (described above with reference to Fig. 12A) applied to n + 1 qubits (n encoding qubits plus one ancilla), followed by phase-adjusting and permutation circuits. Here, single-qubit rotations of the form

are used with the single-qubit phase shift gate P(ϕ) = diag{1, exp(i ϕ)}. The unitary operators U₁ and U₂ applied to the ancilla can be combined into a single-qubit gate U_a = U₂ U₁.

[0246] The quantum Chebyshev transform shown in Fig. 12B can be understood as a specific version of the (discrete) cosine transform. Considering the normalized output states evaluated at specific

with the integer j ∈ [0, N - 1], on the condition of ancilla bit measured in |0〉 state; in this case, each of them corresponds to the (j + 1)^th column of the type-II discrete cosine transform matrix, DCT^II_N, and they form a set of orthonormal basis states.

[0247] Therefore, the above-mentioned Chebyshev feature map circuit takes at least N times to construct a full DCT^II_N matrix. Fig. 12B shows an alternative approach similar to quantum Fourier transform (QFT) to obtain the quantum Chebyshev transform (QChT) equivalent to DCT^II_N.

[0248] The circuit coverts a complex-valued matrix into a purely real matrix through a series of unitary gates. The circuit 1200 starts with a Hadamard gate 1202 on the ancilla (being considered the most significant bit), followed by a sequence of controlled-NOT gates 1204 applied on the other n qubits to create a (n + 1)-qubit entangled state. This is followed by a (n + 1) QFT circuit 1206 applied to both the ancillary and encoding qubits. After the QFT, controlled-R_z gates 1208 are introduced to adjust the relative phase of the encoding qubits, and a unitary operator

is applied to the ancilla.

[0249] The resulting state (behind the U₁ gate) has purely real amplitude of |0〉_a|Φ〉 and a purely imaginary amplitude of |1〉_a|Φ〉 for any n-qubit intermediate state |Φ〉. A permutation circuit 1210 is used to reorder amplitudes of the conditioned states, i.e., |1〉_a |Φ〉 → |1〉_a |(Φ + N - 1) mod N〉, followed by a layer 1212 of n-controlled-NOT gates to ensure that the amplitude of |0〉_a|Φ〉 (|1〉_a|Φ〉) returns to purely real (imaginary). This is followed by a layer 1214 comprising a unitary operator

applied to the ancilla, and then multi-controlled R_X gates. Finally, on the condition of the input and output ancilla bits being respectively initialized with and measured in the |0〉 state (not shown), the quantum Chebyshev transform is obtained.

[0250] A representative example of the Chebyshev transform matrix for n = 2 is given by

where only the top-left corner of a three-qubit unitary is shown, and similar structure block with sine-valued coefficients appears in the right-bottom corner. For a general n-qubit quantum Chebyshev transformation, the unitary transformation can be expressed as

where |τ(x)〉 is given by Equation 11 and

are the Chebyshev nodes.

[0251] To the circuit illustrated in Fig. 12B, a permutation circuit identical to block 1210 can be appended, in order to obtain the type II discrete sine transform. In this case, the input and output ancilla bits are initialized with and measured in |1〉, respectively.

[0252] The quantum Chebyshev transform circuit 1200 can be run in reverse direction to perform the inverse quantum Chebyshev transform, which equivalently maps Chebyshev basis states

to computational basis states

.

[0253] Fig. 13A is a hardware-level schematic of the actions effectuating the logical operations shown in circuit diagrams described herein, e.g., with reference to Figs. 3A, 6A-B, 7A-B, 8-12, and 14-16. Unitary operators, e.g. those to encode the quantum kernel feature map and derivatives thereof, can be decomposed into a sequence of logical gate operations. These logical gate operations are transformations in a quantum Hilbert space over the qubits encoding an input parameter. In order to transform the internal states of these qubits, a classical control stack is used to send pulse information to a pulse controller that affects one or more qubits. The controller may send a sequence of such pulses in time and for each qubit independently. An initialization pulse is used to initialize the qubits into the |0〉 state 1302. Then, for example a series of single-qubit pulses is sent to the qubit array in 1304, which may apply a single-layer feature map. Additionally, two-qubit pulse sequences can be used to effectively entangle multiple qubits with a feature map 1306. The duration, type, strength and shape of these pulses determine the effectuated quantum logical operations. 1308 indicates a 'break' in the depicted timeline, which means the sequence of gates may be repeated in a similar fashion in the direction of the time axis 1312. At the end of the pulse sequences, one or more of the qubits are measured 1310.

[0254] Fig. 13B is a hardware-level schematic of the actions effectuating the logical operations shown in circuit diagrams described herein, specified for a photonic/optical quantum processor. Unitary operators, e.g. those used to encode the quantum kernel feature map and derivatives thereof, can be decomposed into a sequence of optical gate operations. These optical gate operations are transformations in the quantum Hilbert space over the optical modes. In order to transform the internal states of these modes, a classical control stack is used to send pulse information to a pulse controller that affects one or more modes. The controller may formulate the programmable unitary transformations in a parametrised way.

[0255] Initially the modes 1314 are all in the vacuum state 1316, which are then squeezed to produce single-mode squeezed vacuum states 1318. The duration, type, strength and shape of controlled-optical gate transformations determine the effectuated quantum logical operations 1320. At the end of the optical paths, one or more modes are measured with photon-number resolving, Fock basis measurement 1322, tomography or threshold detectors.

[0256] Fig. 13C is a hardware-level schematic of the actions effectuating the logical operations shown in circuit diagrams described herein, specified for a Gaussian boson sampling device. Unitary operators, e.g. those used to encode the quantum kernel feature map and derivatives thereof, can be decomposed into a sequence of optical gate operations. These optical gate operations are transformations in the quantum Hilbert space over the optical modes. In order to transform the internal states of these modes, a classical control stack is used to send information to optical switches and delay lines. The controller may formulate the programmable unitary transformations in a parametrised way.

[0257] Initially the modes 1326 are all in a weak coherent state, which is mostly a vacuum state with a chance of one or two photons and negligibly so for higher counts. Subsequently, the photons travel through optical waveguides 1328 through delay lines 1330 and two-mode couplers 1332 which can be tuned with a classical control stack, and which determines the effectuated quantum logical operations.

[0258] At the end of the optical paths, one or more modes are measured with photon-number resolving 1334, or threshold detectors.

[0259] Fig. 13D is a hardware-level schematic of the actions effectuating the logical operations shown in circuit diagrams described herein, which may be executed on a neutral-atom-based quantum computer. On this type of hardware, unitary operators, e.g. those used to encode the quantum feature map and derivatives thereof, can be decomposed in two different kinds of operations: digital or analog. Both of these operations are transformations in the quantum Hilbert space over atomic states.

[0260] Schematic (a) of Fig. 13D depicts a digital quantum circuit 1338, wherein local laser pulses may be used to individually address neutral atoms to effectuate transitions between atomic states which effectively implement sets of standardized or 'digital' rotations on computational states. These digital gates may include any single-qubit rotations, and a controlled-pauli-Z operation with arbitrary number of control qubits. Additionally, such digital gate operations may also include 2-qubit operations.

[0261] Schematic (b) of Fig. 13D depicts an analog mode 1346 of operation, wherein a global laser light pulse may be applied to groups of, or all, atoms at the same time, with certain properties like detuning, Rabi frequencies and Rydberg interactions to cause multi-qubit entanglement thereby effectively driving the evolution of a Hamiltonian 1344 of the atomic array in an analog way. The combined quantum wavefunction evolves according to Schrödinger's equation, and particular, unitary operators

=

, where

denotes the Hamiltonian and t the time, can be designed by pulse-shaping the parametrised coefficients of the Hamiltonian in time. This way, a parametric analog unitary block can be applied, which entangles the atoms and can act as a variational ansatz, or a feature map, or other entanglement operation.

[0262] The digital and analog modes can be combined or alternated, to yield a combination of the effects of each. Schematic (c) of Fig. 13D depicts an example of such digital-analog quantum circuit, including blocks 1346_1-3 of digital qubit operations (single or multi-qubit) and analog blocks 1348_1-3.

[0263] It can been proven that any computation can be decomposed into a finite set of digital gates, including always at least one multi-qubit digital gate (universality of digital gate sets). This includes being able to simulate general analog Hamiltonian evolutions, by using Trotterization or other simulation methods. However, the cost of Trotterization is expensive, and decomposing multi-qubit Hamiltonian evolution into digital gates is costly in terms of number of operations needed.

[0264] Digital-analog circuits define circuits which are decomposed into both explicitly-digital and explicitly-analog operations. While under the hood, both are implemented as evolutions over controlled system Hamiltonians, the digital ones form a small set of pre-compiled operations, typically but not exclusively on single-qubits, while analog ones are used to evolve the system over its natural Hamiltonian, for example in order to achieve complex entangling dynamics.

[0265] It can be shown that complex multi-qubit analog operations can be reproduced/simulated only with a relatively large number of digital gates, thus posing an advantage for devices that achieve good control of both digital and analog operations, such as neutral atom quantum computer. Entanglement can spread more quickly in terms of wall-clock runtime of a single analog block compared to a sequence of digital gates, especially when considering also the finite connectivity of purely digital devices.

[0266] Further, digital-analog quantum circuits for a neutral quantum processor that are based on Rydberg type of Hamiltonians can be differentiated analytically so that they can be used in variational and/or quantum machine learning schemes, including the differential quantum circuit (DQC) schemes as described in this application.

[0267] In order to transform the internal states of these modes, a classical control stack is used to send information to optical components and lasers. The controller may formulate the programmable unitary transformations in a parametrised way.

[0268] At the end of the unitary transformations, the states of one or more atoms may be read out by applying measurement laser pulses, and then observing the brightness using a camera to spot which atomic qubit is turned 'on' or 'off', 1 or 0. This bit information across the array is then processed further according to the embodiments.

[0269] Fig. 13E is a hardware-level schematic of the actions effectuating the logical operations shown in circuit diagrams described herein, specified for a photonic/optical quantum processor. The quantum model can be decomposed into a sequence of optical gate operations. These optical gate operations are transformations in the quantum Hilbert space of the photons. In order to transform the internal states of these photons, a classical control stack is used to send information to a universal multiport interferometer. The controller may formulate the programmable unitary transformations in a parameterized way.

[0270] Initially the photons 1355 are in Fock states, weak coherent states or coherent states. The duration, type, strength and shape of controlled-optical gate transformations determine the effectuated quantum logical operations 1356.

[0271] At the end of the optical paths, the modes are measured with photon-number resolving, Fock basis measurement 1357, tomography or threshold detectors.

[0272] Fig. 14A depicts a system describing a Digital-Analog implementation of a quantum algorithm. In particular, Fig. 14A represents a Digital-Analog implementation where the qubit interactions can be switched on and off and there exists single-qubit addressability for the application of single-qubit gates, a process described in 1400. In order to transform the internal states of these qubits, a classical control stack may be used to send information to a controller that affects one or more qubits. The controller may send multiple information carriers (e.g., pulses) in time and for each qubit independently. An initialization protocol is used to initialize the qubits, in this example into the |0〉 state 1402. The current example shows only four qubits, but the same principles may be extended in a straightforward manner to systems with more (or less) than four qubits.

[0273] Then, a quantum feature map is applied 1404 to encode quantum information, e.g., an input variable, into a Hilbert space associated with the quantum processor.

[0274] Following application of the feature map, a variational Ansatz is applied, implemented as a variational quantum circuit 1406. In the current example, the variational Ansatz comprises three single-qubit gates R_Y - R_X - R_Y with different rotational angles θ_i. These single-qubit gates are typically sets of standardized or 'digital' rotations on computational states applied to different qubits with different rotation angles/parameters. These digital gates include any single-qubit rotations according to the θ_i argument of the rotation. A single rotation/single-qubit gate is not sufficient to perform arbitrary rotation regardless of the parameter/angle since it can only rotate over a single axis. In order to accomplish a general rotation block, three gates in series with different rotation angles/parameters can be used, in this example represented by a block of single-qubit gates 1406.

[0275] Then, the entanglement in this Digital-Analog approach is generated by a wavefunction evolution 1408, described by the block

. During this evolution, the qubits are interacting amongst themselves, letting the system evolve for a specified amount of time. This process produces the necessary entanglement in the system. The combined quantum wavefunction evolves according to Schrödinger's equation, and particular unitary operators

=

in which

is the Hamiltonian of the system (for example, for neutral atoms the Hamiltonian that governs the system is

with n̂ being the state occupancy of the Rydberg atoms), and t is the evolution time. In this way, a parametric analog unitary block can be applied, which entangles the atoms and can act as a variational Ansatz.

[0276] After the evolution of the wavefunction 1408, another set of single-qubit gates may be applied, similar to the process described in the block of gates R_Y - R_X - R_Y 1406.

[0277] Then, the wavefunction may be evolved once more, and finally a measurement in the computational basis occurs as described in 1410.

[0278] Optionally, additional steps may be included, represented by the ellipses, e.g., before, after, or in between the blocks shown. For example, a different initial state than |0〉 may be prepared prior to application of the feature map. Furthermore, the blocks can be executed in a different order, for instance, in some embodiments, block 1408 might precede block 1406. One or more of the blocks 1404-1408, or variations thereof, may be repeated one or more times prior to the measurement 1410.

[0279] Fig. 14B represents a typical quantum feature map 1412 that starts by preparing the quantum register to the |0〉 quantum state 1402 and then applying single-qubit rotations in all of the qubits with (typically) different rotation angles 1404. In this example, single qubit rotations 1404 (here chosen as R_Y(f_i(x)) act on each qubit individually and are parametrized by a variable x. For a non-linear feature map encoding such as the Chebyshev feature map, the variable x may be used as an angle of rotation.

[0280] The application 1404 of one rotation operation to each qubit in a register may be referred to as a single layer of rotation operations. This type of encoding of classical information into quantum states is known as angle encoding, which means that a data vector is represented by the angles of a quantum state. Angle encoding can be found in many Quantum Machine Learning algorithms, with the main advantage being that it only requires n = log NM qubits to encode a dataset of M inputs with N features each. Thus, considering an algorithm that is polynomial in n, it has a logarithmic runtime dependency on the data size.

[0281] Fig. 14C illustrates different pulses (potentially differing in pulse duration and/or amplitude/height) that can be sent to the quantum circuit (different qubits) to encode the classical information. The angle of the rotation that is described by the f_i(x) dictates the duration and geometry of the pulse.

[0282] Fig. 15A presents a quantum circuit based on a Digital-Analog implementation 1500 that has no single-qubit addressability and uses a wavefunction evolution that can be turned on and off for the entanglement generation. Fig. 15A differs from Fig. 14A in that it shows a technique known as data re-uploading that alters the quantum feature map 1504 and 1510. Loading classical data into a quantum system is a non-trivial task, and uploading large amounts of data makes it an even more critical issue. However, when the case of big-data is not considered, techniques like data re-uploading can be utilized. As the name suggests, this entails uploading the classical data multiple times, and encode them via rotational angle/parameter encoding into single-qubit rotations throughout the quantum circuit. Such a data re-uploading technique may start by setting the quantum register to an initial state, e.g., the all-zero state 1502, followed by angle encoding that is based on a feature map 1504, in this case a tower feature map.

[0283] Following that, as already explained above with reference to 1406, three single-qubit gates in series are used to perform a general rotation block that acts as the Ansatz 1506. In this example the encoding is done through the θ_i variable of the R_Y - R_X - R_Y rotations; however, any type of combined rotations of the type (R_X, R_Y, R_Z) could be applied here. Since there is no single-qubit addressability in this example, each R_Y and R_X gate is applied at all qubits simultaneously with the same rotation angle θ_i. Subsequently, the wavefunction is allowed to evolve for a specified amount of time 1508.

[0284] Then the classical data are re-uploaded through angle encoding into single-qubit rotations 1510. However, each time the (same) data are uploaded, the information is encoded using different angles than in the previous data uploading steps. For example, the amount of rotation may be doubled in each data (re-)uploading step 1504,1510 (i.e., the rotational angle of the might be increasing as 1-2-4, etc.).

[0285] Once again, the rotations R_Y - R_X - R_Y are applied to all qubits simultaneously with the same rotation angle for each single qubit gate, but different than 1506. The classical information can be repeatedly encoded into the quantum feature map through the data re-uploading technique, while changing the angle of the single qubit rotational angle (1504, 1510, etc.) every time after the wavefunction evolution. The resulting quantum feature map can become more expressive as a tower feature map by doing this process of serial data re-uploading with different levels of rotational angles every time the re-uploading is done.

[0286] Fig. 15B illustrates a different (and at current technology levels more realistic) form of a Digital-Analog implementation. In this implementation, the system has single-qubit addressability, but the wavefunction evolution cannot be turned off, thus staying always on. This system requires less control over the qubit interactions (since they are always on) and only requires single-qubit addressability.

[0287] The single-qubit gates are applied while the wavefunction evolution (qubit interaction) is also happening. In the shown example, the quantum circuit 1518 starts by encoding information in the |0〉 quantum state 1502, followed by a quantum feature map 1504 (analogous to 1404).

[0288] Then, a variational Ansatz is applied, comprising both wavefunction evolution

1508 and simultaneous application of a block of single-qubit gates 1516 (similar to block 1406). After multiple repetitions of block 1516, each time with different rotational angles at the single-qubit gates, the circuit ends with a measurement in the computational basis 1514.

[0289] Fig. 16A illustrates a quantum circuit 1600 that implements a Digital-Analog implementation with single-qubit addressability, wherein the wavefunction evolution (e.g., atom interactions) can be turned on and off. This quantum circuit is equivalent to the quantum circuit in Fig. 15A, with the time of the wavefunction evolution 1608 being t = π (in dimensionless coordinates), while all other operations are the same as the ones described in 1500. The same technique of data re-uploading with a tower feature map 1604, 1610 and the same Ansatz comprising single-qubit gates 1606,1612 followed by wavefunction evolution 1608 are used (all similar to 1500).

[0290] The reason that such a specific time period for the wavefunction evolution 1608 is selected, is due to the fact that such evolution results in the application of Controlled-Z (CZ) gates between pairs of qubits. The benefit of describing the wavefunction evolution as a set of CZ gates lies in the fact that this Digital-Analog implementation of a quantum circuit resembles the structure of the equivalent Digital quantum circuit, which instead of the operation 1608 performs a set of CNOT gates between pairs of qubits included in the quantum circuit. This allows more straightforward comparisons to be drawn between the Digital and Analog implementation of such feature maps, and to evaluate their relative performance.

[0291] Fig. 16B illustrates a quantum circuit that implements a Digital-Analog implementation with single-qubit addressability 1618, similar to 1600, but with the entanglement generated via CZ gates between the qubits, resembling the atom interaction / entanglement generation.

[0292] Fig. 17A depicts a system for determining a solution for an computational problem using a hybrid quantum computer system according to an embodiment. The system 1702 may include a quantum computer system 1704 comprising one or more quantum processors 1708, e.g. a gate-based qubit quantum processor, and a controller system 1710 comprising input output (I/O) devices which form an interface between the quantum processor and the outside world, e.g., the one or more classical processors of a classical computer 1706. For example, the controller system may include a system for generating control signals for controlling the quantum processing elements. The control signals may include for example a sequence of pulses, e.g. microwave pulses, voltage pulses and/or optical pulses, which are used to manipulate qubits. Further, the controller may include output device, e.g. readout circuits, for readout of the qubits and the control signals for readout of the quantum processing elements, e.g. a readout pulse for reading a qubit. In some embodiments, at least a part such readout circuit may be located or integrated with the chip that includes the qubits.

[0293] The system may further comprise a (purely classical information) input 1712 and an (purely classical information) output 1714. The data processor systems may be configured to solve an optimization problem using the quantum computer. Input data may include information about the optimisation problem one wishes to solve. This information may include training data such as input-output pairs, a differential equation, boundary conditions, initial values, regularization values, etc. The input data may be used by the system to construct quantum circuits, in particular quantum feature maps and/or parametrized quantum circuits, and to classically calculate values, e.g. sequences of pulses, which may be used to initialize and control qubit operations according to the quantum circuit. To that end, the classical computer may include a quantum circuit generator 1707. The input data may be used by the system to classically calculate values, e.g. parameter settings, which may be used to initialize the quantum circuit that is implemented on the quantum processor. Similarly, output data may include loss function values, sampling results, correlator operator expectation values, optimisation convergence results, optimized quantum circuit parameters and hyperparameters, and other classical data.

[0294] Each of the one or more quantum processors may comprise a set of controllable quantum processing elements, e.g. a set of controllable two-level systems referred to as qubits. The two levels are |0〉 and |1〉 and the wave function of a N-qubit quantum processor may be regarded as a complex-valued superposition of 2^N of these (distinct) basis states. The embodiments in this application however are not limited to qubits but may include any multi-level quantum processing elements, e.g. qutrits, that is suitable for performing quantum computation. Examples of such quantum processors include noisy intermediate-scale quantum (NISQ) computing devices and fault tolerant quantum computing (FTQC) devices.

[0295] The quantum processor may be configured to execute a quantum algorithm in accordance with the gate operations of a quantum circuit. The quantum processor may be implemented as a gate-based qubit quantum device, which allows initialization of the qubits into an initial state, interactions between the qubits by sequentially applying quantum gates between different qubits and subsequent measurement of the qubits' states. To that end, the input devices may be configured to configure the quantum processor in an initial state and to control gates that realize interactions between the qubits. Similarly, the output devices may include readout circuitry for readout of the qubits which may be used to determine a measure of the energy associated with the expectation value of the Hamiltonian of the system taken over the prepared state.

[0296] Furthermore, the one or more quantum processors may comprise a set of continuous variable systems, such as optical or photonic quantum computers. Furthermore, the quantum processor may comprise a Gaussian Boson Sampling device as described above with reference to Fig. 13C or any of the other hardware described above with reference to Fig. 13A-E.

[0297] In some embodiments, the first data processor system may be implemented as a software program for simulating a quantum computer system 1704 comprising a quantum processor system 1708. Hence, in that case, the software program may be a classical software program that runs a classical computer 1706 so that quantum algorithms can be developed, executed and tested on a classical computer without requiring access to a hardware implementation of the quantum processor system.

[0298] Fig. 17B is a block diagram illustrating an exemplary classical data processing system described in this disclosure, for example classical computer 1706. Classical data processing system 1720 may include at least one processor 1722 coupled to memory elements 1724 through a system bus 1726. As such, the data processing system may store program code within memory elements 1724. Further, processor 1722 may execute the program code accessed from memory elements 1724 via system bus 1726. In one aspect, data processing system may be implemented as a computer that is suitable for storing and/or executing program code. It should be appreciated, however, that data processing system 1720 may be implemented in the form of any system including a processor and memory that is capable of performing the functions described within this specification.

[0299] Memory elements 1724 may include one or more physical memory devices such as, for example, local memory 1728 and one or more bulk storage devices 1730. Local memory may refer to random access memory or other non-persistent memory device(s) generally used during actual execution of the program code. A bulk storage device may be implemented as a hard drive or other persistent data storage device. The classical data processing system 1720 may also include one or more cache memories (not shown) that provide temporary storage of at least some program code in order to reduce the number of times program code must be retrieved from bulk storage device 1730 during execution.

[0300] Input/output (I/O) devices depicted as key device 1732 and output device 1734 optionally can be coupled to the data processing system. Examples of key device may include, but are not limited to, for example, a keyboard, a pointing device such as a mouse, or the like. Examples of output device may include, but are not limited to, for example, a monitor or display, speakers, or the like. Key device and/or output device may be coupled to data processing system either directly or through intervening I/O controllers. A network adapter 1736 may also be coupled to data processing system to enable it to become coupled to other systems, computer systems, remote network devices, and/or remote storage devices through intervening private or public networks. The network adapter may comprise a data receiver for receiving data that is transmitted by said systems, devices and/or networks to said data and a data transmitter for transmitting data to said systems, devices and/or networks. Operation modems, cable operation modems, and Ethernet cards are examples of different types of network adapter that may be used with classical data processing system 1720.

[0301] As pictured in FIG. 17B, memory elements 1724 may store an application 1738. It should be appreciated that classical data processing system 1720 may further execute an operating system (not shown) that can facilitate execution of the application. Application, being implemented in the form of executable program code, can be executed by classical data processing system 1720, e.g., by processor 1722. Responsive to executing application, data processing system may be configured to perform one or more operations to be described herein in further detail.

[0302] In one aspect, for example, classical data processing system 1720 may represent a client data processing system. In that case, application 1738 may represent a client application that, when executed, configures classical data processing system 1720 to perform the various functions described herein with reference to a "client". Examples of a client can include, but are not limited to, a personal computer, a portable computer, a mobile phone, or the like.

[0303] In another aspect, data processing system may represent a server. For example, data processing system may represent an (HTTP) server in which case application 1738, when executed, may configure data processing system to perform (HTTP) server operations. In another aspect, data processing system may represent a module, unit or function as referred to in this specification.

[0304] The terminology used herein is for the purpose of describing particular embodiments only and is not intended to be limiting of the invention. As used herein, the singular forms "a," "an," and "the" are intended to include the plural forms as well, unless the context clearly indicates otherwise. It will be further understood that the terms "comprises" and/or "comprising," when used in this specification, specify the presence of stated features, integers, steps, operations, elements, and/or components, but do not preclude the presence or addition of one or more other features, integers, steps, operations, elements, components, and/or groups thereof.

[0305] The corresponding structures, materials, acts, and equivalents of all means or step plus function elements in the claims below are intended to include any structure, material, or act for performing the function in combination with other claimed elements as specifically claimed.

[0306] The description of the present embodiments has been presented for purposes of illustration and description, but is not intended to be exhaustive or limited to the invention in the form disclosed. Many modifications and variations will be apparent to those of ordinary skill in the art without departing from the scope and spirit of the invention. The embodiment was chosen and described in order to best explain the principles of the invention and the practical application, and to enable others of ordinary skill in the art to understand the invention for various embodiments with various modifications as are suited to the particular use contemplated.

Claims

1. A method for determining a solution to a differential problem using a hybrid computer system comprising a classical computer system and a quantum processor, the differential problem defining a differential equation and one or more constraints, the solving comprising finding a target function defining an approximate solution of the differential equation, the target function being represented by a vector of coefficients with respect to a given function basis, the method comprising:

receiving or determining information on one or more quantum circuits encoding a Hamiltonian representing the differential problem with respect to the basis functions, the Hamiltonian comprising a first contribution representing the differential equation and a second contribution representing the one or more constraints; and

determining the solution to the differential problem, the determining comprising:

- initialising at least a first part of the quantum processor in an initial state;

- evolving the at least first part of the quantum processor into a final state, the final state targeting a zero-eigenvalue state of the Hamiltonian, the final state representing the vector of coefficients with respect to the given function basis; and

- receiving measurement data from the quantum processor, to obtain the solution to the differential problem.

2. The method as claimed in claim 1, further comprising evaluating the target function, wherein evaluating the target function comprises:

encoding an independent variable of the target function in a latent space of the quantum processor state by initialising a second part of the quantum processor in a fiducial state and evolving the fiducial state with a feature map associated with the given function basis to obtain an encoded state, the feature map being based on the independent variable; and

measuring an overlap of the encoded state with the final state representing the vector of coefficients with respect to the given function basis.

3. The method as claimed in claim 1, wherein the solution to the differential problem represents a probability density function, and wherein the method further comprises:
generating, by the quantum processor, samples according to the probability density function, the generating the samples comprising:

- applying a unitary transformation associated with the function basis; and

- measuring the quantum processor in a computational basis, a measurement yielding a bitstring, the bitstring representing a sample according to the probability density function.

4. The method as claimed in claim 1 or 2,

wherein the basis functions comprise a constant function; and/or

wherein the basis functions are real-valued functions; and/or

wherein the basis functions define a complete orthonormal function basis;

preferably the basis functions being the Chebyshev basis.

5. The method as claimed in any one of the preceding claims, wherein the Hamiltonian is represented by a sum of a first matrix representing the differential equation and one or more second matrices representing the one or more constraints,

the first matrix being defined as a quadratic form of a third matrix, the third matrix being defined as a sum of matrices corresponding to terms of the differential equation with respect to the basis functions, and

each of the one or more second matrices comprising a quadratic form of a fourth matrix, the fourth matrix corresponding to one of the one or more constraints with respect to the basis functions.

6. The method as claimed in any one of the preceding claims, wherein the one or more constraints are homogeneous data constraints.

7. The method as claimed in any one of the preceding claims, wherein the one or more quantum circuits further define operations for:

encoding a value of an independent variable of the differential equation in a first quantum state using a feature map associated with the basis functions;

encoding the eigenvector of the Hamiltonian with the minimal eigenvalue in a second quantum state; and

measuring an overlap between the first and second quantum states.

8. The method as claimed in any one of the preceding claims, wherein the Hamiltonian comprises a scaling factor, and wherein determining the eigenvector of the Hamiltonian with the minimal eigenvalue comprises:

computing a set of eigenvalues of a part of the Hamiltonian that is scaled by the scaling factor; and

selecting a smallest eigenvalue of the set of eigenvalues.

9. The method as claimed in any one of the preceding claims,
wherein the differential equation is a homogeneous differential equation; and/or wherein the differential equation is a polynomial differential equation.

10. The method as claimed in any one of the preceding claims,

wherein applying a quantum circuit to the quantum processor comprises translating the quantum circuit into a sequence of signals, e.g. optical or electric pulses, and using the sequence of signals to operate quantum elements, e.g. qubits, in a quantum register of the quantum processor; and/or,

wherein receiving measurement data comprises applying a read-out signal to the quantum elements, e.g. qubits, in a quantum register of the quantum processor and, in response to the read-out signal, obtaining data from the quantum hardware representing projective measurements.

11. The method as claimed in any of the preceding claims, wherein the pulses in the sequence of pulses are electromagnetic pulses, e.g., laser pulses, and wherein encoding the quantum circuits in the sequence of pulses comprises determining, for each pulse in the sequence of pulses, one or more of: a pulse duration, a pulse timing, a pulse amplitude, and a pulse phase.

12. The method as claimed in any of the preceding claims,

wherein the hybrid quantum computer includes at least one of: a gate-based qubit device, a digital/analog quantum device, a purely analog quantum device, a continuous-variable operating quantum device and/or a gaussian boson sampling device; and/or,

wherein the hybrid quantum computer comprises at least one of: a neutral atom quantum device, a photonic qubit quantum device, a superconducting qubits quantum device, and/or a trapped ion device.

13. A system for solving a differential problem, the differential problem defining a differential equation and one or more constraints, the solving comprising finding a target function defining an approximate solution of the differential equation, the target function being represented by a vector of coefficients with respect to a given function basis, the system comprising a hybrid computer system comprising a classical computer system and a quantum processor, the hybrid quantum computer system being configured for:

receiving or determining, by the classical computer system, information on one or more quantum circuits encoding a Hamiltonian representing the differential problem with respect to the basis functions; and

determining the solution to the differential problem, the determining comprising:

- initialising at least a first part of the quantum processor in an initial state;

- encoding the quantum circuits and, using the encoded quantum circuits, evolving the at least first part of the quantum processor into a final state, the final state targeting a zero-eigenvalue state of the Hamiltonian, the final state representing the vector of coefficients with respect to the given function basis; and

- receiving measurement data from the quantum processor, to obtain the solution to the differential problem.

14. A computer program product comprising instructions which, when the computer program is executed by a hybrid quantum computer, cause the hybrid quantum computer to execute the method according to any one of claims 1-12.

15. A computer-readable storage medium comprising instructions which, when executed by a hybrid quantum computer, cause the hybrid quantum computer to execute the method according to any one of claims 1-12.

Drawing