Cross-pollinating Quantum Computing and Machine Learning for the Solution of Parameterized Problems in Computational Science and Engineering

1. Introduction

Advancements in quantum computing (QC) are paving the way for a new age in simulations, promising to revolutionize several industry verticals, including communications, finance, medicine, materials, artificial intelligence, cryptography, etc., over the coming decades. The enormous potential for QC to tackle extremely complex issues is garnering ever greater interest as the capabilities of conventional computer systems are rapidly approaching their limits, and efforts globally to mount up the quantum race are intensifying. This disruptive and rapidly emerging new technology represents both a radical change in the laws of physics in computing as we know them and a huge, unimaginable opportunity, which is expected to result in significant changes to many aspects of modern life. Since its first conceptualization in the 1980s and early proof of principles for hardware in the 2000s, quantum computers can now be built with hundreds of qubits thanks to the progress made by companies, start-ups and corporations.

In contrast to classical computing, QC harnesses the laws of quantum mechanics to solve problems that are considered intractable for classical computers by exploiting counter-intuitive quantum physical phenomena like superposition and entanglement. Despite the fundamental differences between classical and quantum computing, any computational problem solvable by a classical computer can also be solved by a quantum computer, and vice versa, as stated by the Church-Turing thesis. However, certain quantum algorithms are demonstrated to have significantly lower time and/or memory complexity than corresponding known classical algorithms and are even expected to solve problems that no classical computer could in any feasible amount of time, a feat known as “quantum supremacy”. Examples include the famous Schor factoring algorithm [1] and the Grover search algorithm [2]. Therefore, the potential of quantum computing for scientific and industrial breakthroughs is immense, which justifies the considerable research devoted to its application in real-world problems.

Even though quantum technology is still in its infancy, the rapid development of quantum hardware and the enormous monetary investments made worldwide have led many to claim that the so-called Noisy-Intermediate Scale Quantum (NISQ) devices already in existence may soon outperform conventional computers. In short, NISQ devices are near-term quantum computers with insufficient physical qubits to execute effective error-correcting techniques. On a handful of problems created specifically for assessing the capabilities of quantum computers, it already has been demonstrated that existing NISQ computers outperform classical computers. The algorithms that run on these constrained devices may only need a few qubits, exhibit some noise resilience, and are frequently implemented as hybrid algorithms, where some operations are carried out on a quantum device and some on a classical computer. Particularly, the number of operations, or quantum gates, must be kept to a minimum since the more errors that are introduced into the quantum state during implementation, the greater the likelihood that the quantum state may decohere. Due to these restrictions, there are limits on the scope of algorithms that can be considered. Therefore, the crucial technological question is how to best utilize current NISQ technology to get a quantum advantage. Any such strategy must account for a limited number of qubits, limited connectivity of the qubits, and coherent and incoherent errors that limit quantum circuit depth.

In parallel, the field of computational science and engineering (CSE) is also exponentially growing, being inextricably linked to the quest of solving problems of rising complexity. In this environment, quantum computing is expected to bring tremendous breakthroughs to this field. Even to this day, many critical problems lay beyond our current computational capabilities and are deemed intractable even for the world’s most advanced supercomputers. Some examples are non-convex optimization problems, modeling physical systems at multiple length and/or time scales, molecular modeling, robust design problems, and others. With quantum computing coming to the fore, there is a silver lining in the computing world that one can finally solve such complex problems. Towards this direction, quantum variants of classical algorithms have been recently proposed for solving linear systems of algebraic equations [3] and differential equations [4], as well as for performing optimization [5,6] and machine learning [7,8,9] in quantum computers. All these problems lie at the heart of CSE and incorporating new quantum technologies may lead to unforeseeable advances in the field. Nevertheless, we are still in the NISQ era and there are still several obstacles and limitations that hinder the application of quantum algorithms in the field of CSE. Current state-of-the-art quantum technologies solve at best a very limited set of problems and it’s unlikely that there will be available general-purpose quantum computers in the near future that can solve problems more quickly than classical computers do. In addition, even quantum algorithms that are known to be faster than their classical counterparts most often require a fair amount of preparation done on a classical computer, either to set up the problem and/or interpret the answer.

Based on these, it becomes evident that quantum computing utilization requires efficient synergy with classical computing. In this regard, variational quantum algorithms (VQAs) have emerged as the leading strategy to obtain a quantum advantage on NISQ devices. VQAs are hybrid quantum-classical algorithms that consist of a parameterized quantum circuit, that is an ansatz of gates containing tunable parameters, and an optimization part performed on a classical computer to find the optimal values of the parameters in the circuit. Most prominent examples in this category include the VQLS algorithm [10,11] for solving linear systems of equations, the Quantum Approximate Optimization algorithm (QAOA) [5] and the Quantum Analytic Descent algorithm (QAD) [12] for optimization problems, as well as the Variational Quantum Eigensolver (VQE) [13], tasked with finding the lowest eigenvalue of a given matrix.

This work proposes a methodology that fuses hybrid quantum algorithms and in particular the VQLS algorithm, with machine learning to accelerate the solution of parameterized linear systems. Such systems arise frequently in CSE applications in the context of uncertainty propagation, parameter inference, optimization, or sensitivity analysis, since these problems are governed by parameterized partial differential equations that are converted to linear systems using numerical discretization methods (typically, finite element or finite difference methods). A general approach to treat these types of problems is to evaluate the system response for multiple parameter instances and post-process the results to extract relevant information on quantities of interest. However, this can be a computationally intensive process, especially when the linear systems under consideration are of high dimensionality. In this regard, this work explores the use of the VQLS algorithm as the main tool to solve parametrized linear systems and investigates ways to accelerate its convergence for a given snapshot of parameters by employing machine learning algorithms to efficiently deliver accurate predictions on the optimal circuit parameters of the corresponding linear system. The main intuition behind this approach is that successful initial prediction on the circuit parameters will reduce the number of iterations of the classical optimizer employed by the VQLS algorithm.

The paper is structured as follows: In Section 2 the fundamental principles of quantum computing are presented. Section 3 provides the mathematical and algorithmic background for the VQLS algorithm. Section 4 presents a methodology that utilizes machine learning to accelerate the convergence of the VQLS algorithm when used for solving parameterized linear systems. Section 5 demonstrates the proposed methodology in a set of numerical applications and, lastly, Section 5 summarizes the conclusions drawn from this work and discusses future research directions.

2. Basic Principles of Quantum Computing

Initially, we shall introduce the basic principles of quantum computing and specifically the well-established Circuit /or Gate paradigm, on which VQLS is based. Most famous algorithms, s.a. Shor’s for prime factoring [1] or HHL for linear systems [3], follow this paradigm.

2.1. Quantum Computation

A quantum computation is composed of three steps: input state preparation, transformations acting on this state and, finally, measurement of the output state. Due to the stochastic nature of quantum objects, the only way to gain part of the information hidden in a quantum state is by sampling, i.e. repeating the quantum algorithm several times and extracting probabilistic measures.

The elementary memory unit of QC is the qubit, a 2-level quantum system that can be prepared, manipulated and measured in a predefined way. Unlike a classical bit, which can take the distinct value of either 0 or 1, a qubit can be a linear combination — or superposition— of two orthogonal quantum states |0〉 and |1〉. Thus, using Dirac notation, a single-qubit system can be described as

$$|\psi \rangle =a|0\rangle +b|1\rangle ,$$

(1)

where the coefficients $a,b\in \mathbb{C}$ are called amplitudes of the state adhering to the constraint ${|a|}^2+{|b|}^2=1$.

Since not much can be achieved with a single-qubit system, multi-qubit systems can be produced by taking advantage of the peculiar quantum-mechanical properties of superposition and entanglement. Hence, an n-qubit state, which resides in a $2^n$-dimensional complex Hilbert space, can be written as

$$|\psi \rangle =\sum _{i=|0\rangle }^{2^n-1}{a}^i|i\rangle ,$$

(2)

where $a^i\in \mathbb{C},i=1,\cdots ,n$ the amplitudes following the normalization constraint $\langle \psi |\psi \rangle ={\sum }_i{|a^i|}^2=1$. The probability of measuring a particular outcome, $p^i$, is given by the squared magnitude of the corresponding amplitude, $|a^i{|}^2$. As mentioned, in order to calculate that probability, one has to repeatedly take measurements of output state.

Quantum computers use quantum gates, which are elementary operations that manipulate the quantum state of a system, to perform calculations. some of the most notable single-qubit gates are the Pauli gates $X,Y,Z$ and the Hadamard gate H. These gates are represented by unitary matrices, denoted as U and preserve the inner product between quantum states.

An example of a single-qubit circuit is illustrated in Figure 1. Starting from state |0〉, a Hadamard gate is applied to the qubit, followed by a Pauli-Y gate. The resulting state is $|\psi \rangle =YH|0\rangle =-i/\sqrt{2}(|0\rangle -|1\rangle )$. The measurement that determines whether the resulting state is $ket0$ occurs last. A measurement is a non-unitary operation which results with the qubit in either state |0〉 or |1〉. Due to the Hadamard gate, the probability of measuring the state |0〉 is $p(|0\rangle )={\left|-i/\sqrt{2}\right|}^2=0.5$.

In general, multi-qubit state manipulation is achieved by use of unitary operators which evolve an initial state to a target state.

$$|{\psi }_t\rangle =U|{\psi }_0\rangle ,$$

(3)

where $U\in {\mathbb{C}}^{4^n}$. Quantum computers also utilize entanglement, a phenomenon where two or more qubits become physically correlated. In circuit computing, entanglement is introduced via controlled operations. Every unitary operator $U\in {\mathbb{C}}^{4^n}$ has a controlled form

$$CU=|0\rangle \langle 0|\otimes I+|1\rangle \langle 1|\otimes U\in {\mathbb{C}}^{4^{n+1}},$$

(4)

which applies U conditionally on a target-qubit register, if a control-qubit has a |1〉 component. An example of a multi-qubit, controlled gate is illustrated in Figure 2, where the first qubit acts as control-qubit and the rest n qubits as targets.

2.2. Quantum Linear System Problem

When given a Linear System Problem (LSP) of N linear equations and N unknowns which can be expressed as $\mathbf{A}\mathbf{x}=\mathbf{b}$, where x is the vector of unknowns, A the matrix of coefficients and b is the right-hand side vector, one can construct the Quantum Linear System Problem (QLSP) [14,15]

$$A|x\rangle =|b\rangle ,$$

(5)

where $|b\rangle =b/{\parallel b\parallel }_2$ is the normalized right-hand side and A the normalized matrix satisfying $\parallel A\parallel \le 1$.

Although QLSP is practically a normalized version of the LSP, they are two distinct problems with latter being addressed mostly by quantum algorithms. The main difference is the requirement that both the input and output of a QLSP are given as quantum states. This means that any efficient algorithm for QLSP requires an efficient preparation of |b〉 and efficient readout of |x〉, both of which are open problems for the general case.

In the quantum linear systems algorithm literature, the term "efficient" usually refers to polylogarithmic complexity w.r.t the system size N. This means that Quantum algorithms like [16] solve the QLSP exponentially faster than even the fastest classical algorithms like Conjugate Gradient method. Solving and reading in quantum computing, though, are two completely distinct processes and as such, we must note though that in order to read out all the elements of |x〉, we would still require $\mathcal{O}\left(N\right)$ time, just as any efficient classical algorithm. Hence, a solution to QLSP should be used as a subroutine in an application where samples from vector x are useful, as discussed thoroughly in [17].

3. The VQLS algorithm

The Variational Quantum Linear Solver (VQLS) [10] is a Variational Hybrid Quantum-Classical (VHQC) algorithm used for solving QLSPs and designed with NISQ devices in mind. Given the QLSP $A|x\rangle =|b\rangle$, the quantum part of the algorithm uses a cost function which captures the similarity between vectors $A|x\rangle$ and |b〉, by means of either a local or global Hamiltonian. Classical optimization schemes can be then utilized in order to minimize the Hamiltonian loss and finally terminate when a desired precision $ϵ$ is reached.

Although variational algorithms are heuristic and complexity analysis of VQLS is a challenging open problem, the authors of [10] provided numerical simulations indicating efficient scaling in the linear system’s dimension N, condition number $\kappa$ and desired absolute error $ϵ$ w.r.t the exact solution. Specifically, they found evidence of (at worst) linear scaling in $\kappa$, logarithmic scaling in $1/ϵ$ and polylogarithmic scaling in N [10].

3.1. Overview

The algorithm requires the following as input: an efficient operator initializing a quantum state $|b\rangle =U|0\rangle$, a vector of initial parameters $\mathit{\alpha }$ and a linear decomposition of the normalized matrix $A={\sum }_{l=1}^L{c}_l{A}_l$, with $A_l$ being unitary matrices and $c_l\in \mathbb{C}$. Intuitively, this linear decomposition into unitaries is essential to construct the required circuits since all quantum gates are unitary operators themselves. VQLS requires the following assumptions be made: L scales as $\mathcal{O}\left(\mathrm{poly}\right(n\left)\right)$, A has a finite condition number $\kappa <\infty$, its spectral norm satisfies $\parallel A\parallel \le 1$. The $A_l$ unitaries must additionally be implemented with efficient circuits, which in quantum computing terms, this usually implies a worst case polynomial dependence on n. Information over efficient decomposition method for sparse A can be found at [10]. The output is the quantum state |x〉 that approximately satisfies the system $\mathbf{A}\mathbf{x}=\mathbf{b}$ up to a constant factor. Parameters $\mathit{\alpha }$ in the ansatz $V\left(\mathit{\alpha }\right)$ are adjusted in a hybrid, quantum and classical optimization scheme until the cost $C\left(\mathit{\alpha }\right)$ (local or global) reaches a predefined threshold. Lastly, the optimized ansatz is applied to construct the solution state which may be measured in order to compute observable quantities of interest. An overview of the algorithm can be seen in Figure 3.

As mentioned, in order to solve the QLSP one must prepare a state |x〉 s.t. $A|x\rangle$ is as close as possible to |b〉. In order to achieve that, VQLS employs an ansatz $V\left(\mathit{\alpha }\right)$ that prepares a candidate solution $|x(\mathit{\alpha })\rangle =V\left(\mathit{\alpha }\right)|0\rangle$. A parameter vector $\mathit{\alpha }$ is given as input to the quantum computer, which prepares $|x(\mathit{\alpha })\rangle$ and runs an efficient quantum circuit that estimates a cost function $C\left(\mathit{\alpha }\right)$. The cost function’s value quantifies how far the proposed state $A|x\rangle$ deviates from the desired state |b〉. This value $C\left(\mathit{\alpha }\right)$ is then returned to the classical computer which adjusts parameters $\mathit{\alpha }$ according to a classical optimization algorithm, attempting to reduce the cost. The process is repeated until a termination condition of the form $C\left(\mathit{\alpha }\right)\le \gamma$ is achieved, at which point VQLS outputs the solution parameters ${\mathit{\alpha }}_{\mathrm{opt}}$.

Finally, the optimal parameters ${\mathit{\alpha }}_{\mathrm{opt}}$ are given as input to the ansatz which computes the solution state $|x({\mathit{\alpha }}_{\mathrm{opt}})\rangle =V\left({\mathit{\alpha }}_{\mathrm{opt}}\right)|0\rangle$. Once the solution state is achieved we can finally measure quantities of interest in order to extract valuable information for the solution vector. The deviation of observable expectation values from those of the exact solution can be upper bounded based on the cost function.

3.1.1. VQLS Cost functions

According to [10], the cost functions can be described as expectation values of either global or local Hamiltonians. To avoid cumbersome notation, we will thereafter write $|\psi \rangle =A|x\rangle =A|x(\mathit{\alpha })\rangle$. As it is noted in [18], global functions may introduce barren plateaus in the energy field, which lead to exponentially vanishing gradients with respect to the number of qubits n. For improved trainability as n grows larger, the authors introduce local Hamiltonian, the normalized version of which is

$$C_L=\frac{\langle \psi |H_L|\psi \rangle }{\langle \psi |\psi \rangle }=\langle \Psi \left(\mathit{\alpha }\right)|H_L|\Psi \left(\mathit{\alpha }\right)\rangle ,$$

(6)

where $|\Psi \rangle =\frac{|\psi \rangle }{\parallel |\psi \rangle \parallel }$ is normalized. The effective Hamiltonian is

$$H_L=U\left(I-\frac{1}{n}\sum _{j=1}^n|0_j\rangle \langle 0_j|\otimes I_{\overline{j}}\right)U^{†},$$

(7)

and $|0_j\rangle$ is the zero state on qubit j and $I_{\overline{j}}$ the identity on all qubits but qubit j.

3.2. Cost evaluation

In order to estimate the cost $C_L$ one must evaluate each part of the fraction in Equation 6 separately, with the first and simpler term to evaluate being

$$\langle \psi |\psi \rangle =\sum _{l{l}^{\prime }}^L{c}_l{c}_{l^{\prime }}^{*}{\beta }_{l{l}^{\prime }}$$

(8)

where

$${\beta }_{l{l}^{\prime }}=\langle \mathbf{0}|V^{†}{A}_{l^{\prime }}^{†}{A}_{l}V|\mathbf{0}\rangle .$$

(9)

Also, for $C_L$ one needs to estimate the terms

$${\delta }_{l{l}^{\prime }}^j=\langle \mathbf{0}|V^{†}{A}_{l^{\prime }}^{†}U(|0_j\rangle \langle 0_j|\otimes I_{\overline{j}})U^{†}{A}_{l}V|\mathbf{0}\rangle ,$$

(10)

where index j implies application of corresponding operator to the j-th qubit, while $\overline{j}$ refers to every other qubit but j. By exploiting the fact that $|0_j\rangle \langle 0_j|=\frac{1}{2}(Z_j+I_j)$, 10 can be written as

$${\delta }_{l{l}^{\prime }}^j=\frac{1}{2}\left(\langle \mathbf{0}|V^{†}{A}_{l^{\prime }}^{†}U\left((I_j+Z_j)\otimes I_{\overline{j}}\right)U^{†}{A}_{l}V|\mathbf{0}\rangle \right)=\frac{1}{2}({\beta }_{l{l}^{\prime }}+{\zeta }_{l{l}^{\prime }}^j),$$

(11)

where

$${\zeta }_{l{l}^{\prime }}^j=\langle \mathbf{0}|V^{†}{A}_{l^{\prime }}^{†}U\left(Z_j\otimes I_{\overline{j}}\right)U^{†}{A}_{l}V|\mathbf{0}\rangle$$

(12)

Taking the previous into account, the cost function finally becomes

$$C_L=\frac{1}{2}-\frac{1}{2n}\frac{{\sum }_j^n{\sum }_{l{l}^{\prime }}^L{c}_l{c}_{l^{\prime }}^{*}{\zeta }_{l{l}^{\prime }}^j}{{\sum }_{l{l}^{\prime }}^L{c}_l{c}_{l^{\prime }}^{*}{\beta }_{l{l}^{\prime }}},$$

(13)

where each of the terms ${\beta }_{l{l}^{\prime }},{\zeta }_{l{l}^{\prime }}^i$ are computed by postprocessing the output of the circuits shown in Figure 4 and Figure 5 respectively.

The circuit in Figure 4 applies operator V on |0〉 and then, conditioned on the first qubit, applies two controlled operations, $C{A}_l$ followed by $C{A}_{l^{\prime }}^{†}$. Finally, a measurement takes place in the ancilla circuit. The $S^{†}$ gate is added to the circuit when the imaginary part of the expectation value is needed. In a similar manner, the circuit of Figure 10 additionally applies the $U^{†}$, $C{Z}_j$ and U operators. The case shown corresponds to $j=1$, when controlled-Z gate is applied to first qubit of working register. As before the $S^{†}$ gate is added to the circuit when the imaginary part of the expectation value is needed.

3.3. Ansatz

VQLS works by preparing the state $|x\rangle =V(\mathit{\alpha })|\mathbf{0}\rangle$. The unitary operator $V\left(\mathit{\alpha }\right)$ consists of a sequence of gates with trainable continuous parameters $\mathit{\alpha }=({\alpha }_1,\cdots ,{\alpha }_L)$ and can more or less be expressed as

$$V\left(\mathit{\alpha }\right)=G_{k_L}\left({\alpha }_L\right)\cdots G_{k^i}\left({\alpha }_i\right)\cdots G_{k_1}\left({\alpha }_1\right).$$

(14)

The discrete parameters $\mathbf{k}=(k_L,\cdots ,k_1)$ encode the types of gates and their placement in the circuit. Such empirical gate constructs are called ansatzes and are common in variational quantum algorithms in general.

If k is constant, we are referring to "fixed-structure" ansatz, since one optimizes only over $\mathit{\alpha }$. An example of fixed ansatz can be seen Figure 6. Although less versatile, fixed ansatzes are simpler than "variable-structure" ones, where one optimizes with respect to k as well. Lastly, another type that can be employed is the Quantum Alternating Operator Ansatz (QAOA) [19,20], which is known to be universal as the number of its layers tend to infinity [21,22].

4. AI-VQLS-solver: coupling Quantum Computing with machine learning to accelerate the solution of parameterized linear systems

In a multitude of applications, such as uncertainty propagation, parameter inference, optimization, sensitivity analysis etc. it is required to produce sufficient number of samples or solution datasets. These, in turn, require the solution of parametric linear systems, and besides the often unavoidable, immense matrix sizes, one has to deal with an overwhelming number of LSP instances which often exceed available computational resources. While quantum algorithms, and specifically VQLS, are expected to deal with the former challenge by exponentially reducing the computational complexity of individual LSPs, there still remains space for further acceleration of the overall parametric problem.

In this regard, it is essential to derive efficient methods that exploit the continuity of the LSP solution space with the respect to the parameters that produce the linear system, hereafter called physical parameters $\mathbf{p}\in {\mathbb{R}}^d$. Starting from the linear system with equation $\mathbf{A}\left(\mathbf{p}\right)·\mathbf{x}=\mathbf{b}$, with $\mathbf{A}\left(\mathbf{p}\right)$ being the parametrized system matrix, x the solution sought and b the right-hand side vector. The equivalent QLSP is the system with equation $A\left(\mathbf{p}\right)|x\rangle =|b\rangle$. Finally, when a variational quantum algorithm such as VQLS is used as solver, the system takes the form

$$A\left(\mathbf{p}\right)|x(\mathit{\alpha })\rangle =|b\rangle ,$$

(15)

where $\mathit{\alpha }=({\alpha }_1,\cdots ,{\alpha }_r)\in {\mathbb{R}}^r$ the VQLS real parameter-vector and $|x(\mathit{\alpha })\rangle$ the proposed solution of the QLSP. Using the cost function from Equation 6 the optimal parameters are obtained by the following minimization problem

$${\mathit{\alpha }}_{\mathrm{opt}}=arg\underset{\mathit{\alpha }}{min}C(\mathbf{p},\mathit{\alpha })=arg\underset{\mathit{\alpha }}{min}\langle \Psi (\mathbf{p},\mathit{\alpha })H_L|\Psi (\mathbf{p},\mathit{\alpha })\rangle ,$$

(16)

It is evident from Equation 16 that for every matrix $A\left(\mathbf{p}\right)$, a search over parametric space ${\mathbb{R}}^r$ is required in order to reach the solution $|x({\mathit{\alpha }}_{\mathrm{opt}})\rangle$ which minimizes the Hamiltonian cost $C_L$.

Our proposed scheme uses a machine-learning model $\mathcal{M}:{\mathbb{R}}^d\to {\mathbb{R}}^r$ in order to propose refined initial guesses to the VQLS algorithm. Instead of letting the initial VQLS parameters be random or constant, we take advantage of the assumed continuity of the cost function $C(\mathbf{p},\mathit{\alpha })$, so that for every physical parameter p, a refined initial guess ${\mathit{\alpha }}_0$ is proposed that is as close as possible to the optimal ones ${\mathit{\alpha }}_0-{\mathit{\alpha }}_{\mathrm{opt}}\le ϵ$ with some probability $P_{\mathcal{M}}$ depending on the chosen model. The acceptable error $ϵ$ between the proposed and optimal parameters is a hyperparameter for the overall scheme and depends on the value of the training loss.

Since, in general, no training dataset will be available, one has to initially rely on the usual workflow without refined initial parameter proposition i.e. ${\mathit{\alpha }}_0^i=\mathcal{S}\left({\mathbf{p}}^i\right)$ for $i=1,\cdots ,N_{\mathrm{train}}$ using some generic strategy $\mathcal{S}$. Once the first optimal solutions are produced, the model can be trained on the respective dataset $({\mathbf{p}}^i,{\mathit{\alpha }}_{\mathrm{opt}}^i),i=1,\cdots ,N_{\mathrm{train}}$. The trained model thereafter proposes refined initial parameters ${\mathit{\alpha }}_0^i=\mathcal{M}\left({\mathbf{p}}^i\right)$ for $i=N_{\mathrm{train}}+1,\cdots ,N_{\mathrm{total}}$. Again, $N_{\mathrm{train}}$ is considered a hyperparameter which must be chosen to balance the trade-off between accuracy and training-time.

The overall scheme is presented in Figure 7. Initially, as noted, the scheme follows the dashed arrows generating the training data by solving the first $N_{\mathrm{train}}$ QLSPs with the base-VQLS. Once the training dataset has been generated, the machine-learning model is trained with $x^i$ as input and ${\alpha }^i$ as output. The computational scheme thereafter changes as indicated by the solid arrows, with the model now proposing refined initial VQLS parameters ${\alpha }^i$, which are closer to the optimal, thus reducing the optimization cost for a QLSP instance.

In the model architecture context, in order to minimize training time without adding to the overall execution time, it is desirable to choose a ML model $\mathcal{M}$ with low complexity and efficient scaling with respect to its output. In our scheme, the only hyperparameter that can be controlled in this regard is the complexity of the ansatz V and more specifically the dimension r of the ansatz parameter vector $\mathit{\alpha }$. An ansatz that scales polynomially with the number of qubits is desirable because it enables one to construct models that remain efficient as the size of the system increases. Specifically, in our case r should scale at most as $\mathcal{O}\left(\mathrm{poly}\right(n\left)\right)$, where n the number of qubits required for the system, so that low-complexity, efficient ML models are sustainable. Some known ansatzes with polynomial qubit scaling include the Quantum Approximate Optimization Algorithm (QAOA) ansatz [5] and the Hardware-Efficient ansatz [23].

5. Numerical examples

To evaluate the performance of the proposed scheme we conducted the comparative experiments described in the subsequent section. The quantum circuit simulations were driven by the Qiskit open-source framework and specifically the statevector simulator [24]. The matrix generator for the parametric $2^n\times 2^n$ QLSP is chosen to be

$$A\left(\mathbf{p}\right)=\frac{1}{\xi }\left(I+\frac{1}{C}\sum _{j,k=1}^n{c}_{ij}\left(\mathbf{p}\right)(X_j\otimes Z_k)\right),$$

(17)

where n is the number of effective qubits used in the circuit, $c_{ij}\left(\mathbf{p}\right)$ are coefficients that depend on the physical parameters p, while $C,\xi$ are normalization constants. In both examples, the physical parameters $\mathbf{p}={(p_0,p_1,p_2)}^T$ are uniformly distributed in a 3-dimensional unit cube $D:[0,1]\times [1,2]\times [3,4]$. At each example we compare pure VQLS schemes with ones that fall under our paradigm. Since the circuits are simulated, the comparison is conducted in relation to the required number of iterations to solve each instance of the parametric system. The parametric QLSP is discretized into a total of $N_{\mathrm{total}}$ regular QLSPs, with equations

$$A\left({\mathbf{p}}^i\right)|x({\mathit{\alpha }}^i)\rangle =b,$$

(18)

and $i=1,\cdots ,N_{\mathrm{total}}$. For all examples the right-hand side is set to $|b\rangle =|\mathbf{0}\rangle \in {\mathbb{C}}^{2^n}$ and the total number of solutions $N_{\mathrm{total}}=1000$.

5.1. Convergence and continuity

A main problem when using variational algorithms is that the Hamiltonian loss-functions and especially their local versions are inherently non-convex which increases the probability for different optimization paths with neighboring initial parameters to land astray at completely different local-minima. Furthermore, quantum ansatzes are inherently periodic, since most parametric quantum gates themselves are not bijective with respect to their parameters. In fact, quantum gates are essentially rotations with period of at most $4\pi$ and the parameters themselves are rather angles of rotation in some complex Hilbert space. Thus, when an ansatz is deployed, in general there is often a multitude of parameter values that can produce the same state. Besides periodicity, an ansatz contains many different paths through which the quantum states are prepared, resulting in different quantum states that may possess the same properties and more specifically the same probability distribution when measured under a specific basis. Lastly, when dealing with real quantum processors, the unavoidable hardware noise will most definitely propagate through the optimization path resulting again in discontinuous optimal values.

Another source of discontinuities which may arise in low-qubit systems is over-parametrization of the quantum ansatz. The required number of parameters to fully span an n-qubit system’s Hilbert space is $N_{\mathrm{span}}=\mathcal{O}\left(2^n\right)$, whereas the number of parameters of an efficient ansatz usually scales as $N_{\mathrm{ans}}=\mathcal{O}\left(\mathrm{poly}\left(n\right)\right)$. Hence, as $n\to \infty$ more often than not the optimization problem will be under-parametrized, since $N_{\mathrm{ans}}\ll N_{\mathrm{span}}$. In smaller systems, though, where the previous inequality is not guaranteed, one might eventually use more parameters than required. This will lead to an underdetermined system with infinite solutions w.r.t a portion of its parameters, giving optimizers more paths to explore, but also more ways to diverge from previous strategies. A the simple illustrative example of this behaviour can be seen in Figure 9. The discontinuities in Figure 9a are due to the periodicity of the non-convex Hamiltonian loss which allow jumps between different local minima. In Figure 9b which has more than required ansatz parameters, the system becomes underdetermined resulting in more discontinuities.

From the above, it is apparent that in order to ensure precision and suffice to simple, scalable model architectures, it is imperative to take into account any potential hindrances posed by aforementioned discontinuities, which can obstruct training and prediction to a substantial extent. Coherently, various classical optimizations schemes used in VQLS may introduce different convergence behaviours with respect to the utilized algorithm. When dealing with convergence alone, the only optimization property that matters is the execution time, where the fastest algorithm to reach a minimum loss, can be considered best for the given problem. In our case though, continuity is of equal importance, since we seek to not only reach a solution, but to force systems with close physical parameters ${\mathbf{p}}^i,{\mathbf{p}}^j\in {\mathbb{R}}^d$ to produce output in the ansatz parameter space ${\mathit{\alpha }}_{\mathrm{opt}}^i,{\mathit{\alpha }}_{\mathrm{opt}}^j\in {\mathbb{R}}^r$ in the same neighbourhood. An optimization algorithm that favours explorability instead of exploitation has a higher chance to visit different local minima during the optimization of concurrent, similar QLSPs and thus resulting in optimal parameter discontinuities with higher probability. In the subsequent experiments, we were confined to gradient-free optimizers in order to evaluate their sufficiency with respect to the potential discontinuities.

5.2. Example 1: two-qubit system

The first example consists of the QLSP in Equation 17 for $n=2$, which results in a parametric system of size $4\times 4$, with the coefficients $c_{ij}$ chosen as trigonometric functions of the physical parameters $p_1,p_2,p_3$ e.g. $c_{00}=cos(p_1+p_2)$ . For the solution of the QLSPs a BFGS optimizer was used with numerical first and second order derivatives approximation. The ansatz used in the example is shown in Figure 10

The Hamiltonian cost function for a given QLSP instance can be seen in Figure 11. We observe an intensive non-convex behariour and periodicity resulting in different local minima which may have hindered optimization algorithm to consistently extract continuous optimal solutions ${\mathit{\alpha }}_{\mathrm{opt}}^i$ across all QLSP instances.

In this example, three strategies were implemented: in the first two the basic VQLS scheme was used with constant, ${\mathit{\alpha }}_0^i=0$, and uniformly random, ${\mathit{\alpha }}_0^i\sim {\mathcal{U}}_{[0,1]}$, proposed initial parameters respectively, while the third strategy consists of a naive nearby scheme, where the proposed initial parameters for each linear system were the optimal parameters of the previous one ${\mathit{\alpha }}_0^i={\mathit{\alpha }}_{\mathrm{opt}}^{i-1}$. The nearby strategy requires an initial reordering of the QLSPs w.r.t the dimensions of the physical parameters. The latter strategy takes advantage of the continuity of the QLSP and requires that ${\mathbf{p}}^i$ are mostly sorted at each dimension in order to avoid large gaps between the sought optimal ansatz parameters. This scheme is quite efficient in the sense that it does not require any training, $N_{\mathrm{train}}=1$, prior to its application and can be compared in a direct manner with the base VQLS.

In Figure 12 one can see the joint probability distributions of the optimal parameters ${\alpha }_{\mathrm{opt},k}^i$ for $k=0,1,2$, $i=1,\cdots ,1000$ for the $2^2\times 2^2$ example of this section. In this example the dataset had no discontinuities which led to a trainable ML model.

The comparative results of the three strategies can be summarized in Figure 13, where the required iterations distribution for solving one instance of the parametric QLSP $A\left({\mathbf{p}}^i\right)|x({\mathit{\alpha }}^i)\rangle =|b\rangle$ are shown. On the left is the pure-VQLS scheme with constant initial ansatz parameters, in the middle a similiar approach but with random initialization and on the right the naive nearby approach, where the proposed ansatz parameters for each system is taken as the optimal of the previous. The speedup factor of the latter over first two schemes based on their median values is $S\approx 1.9$.

5.3. Example 2: three-qubit system

For the second example, a three-qubit system was chosen, which results in a parametric system of size $8\times 8$, with the coefficients $c_{ij}$ being polynomial functions of the physical parameters $p_i$ e.g $c_{00}=p_1^2+p_0^2$. For the solution of the QLSPs a BFGS optimizer was used with numerical first and second order derivatives approximation. The 2-layered ansatz used for this example takes as input a vector $\mathit{\alpha }\in {\mathbb{R}}^{11}$. Each angle $a_k,k=0,\cdots ,10$ is used as a parameter for the $R_y$ gates which are structured as shown in Figure 14. For the solution of the systems a BFGS optimizer was used with numerical first and second order derivatives approximation.

Three distinct strategies were compared, the first two are the pure VQLS scheme with constant initial parameters and the naive nearby scheme, as before. Additionally, as a third approach a simple MLP with 2 hidden layers was trained in the first $N_{\mathrm{train}}=100$ samples using standard MSE loss and an Adam optimizer. In contrast to the nearby scheme, the optimal parameters used for training must cover the whole solution space in order to better approximate the manifold they reside in. In Figure 15, one can observe the non-convex Hamiltonian cost function $C_L$ for an instance of the QLSP. As in the first example, the Hamiltonian is non-convex resulting in different local minima, each one with similar loss values. This behaviour is expected when using local cost functions $C_L$.

In Figure 12 one can see the joint probability distributions of the optimal parameters ${\alpha }_{\mathrm{opt},k}^i$ for $k=0,1,2$ and $i=1,\cdots ,1000$ for the $2^3\times 2^3$ example of this section. In this example the dataset had no discontinuities which led to a trainable ML model.

Figure 16. Pairwise joint probability distributions of optimal parameters ${\alpha }_{\mathrm{opt},k}^i,k=0,1,2,i=1,\cdots ,1000$ for example 5.3. In this example there are no discontinuities resulting in easier training of ML models.

Figure 16. Pairwise joint probability distributions of optimal parameters ${\alpha }_{\mathrm{opt},k}^i,k=0,1,2,i=1,\cdots ,1000$ for example 5.3. In this example there are no discontinuities resulting in easier training of ML models.

The numerical results of the three strategies Constant, Nearby and MLP, are summarized in Figure 17. For each strategy, the graph shows the quartile distribution of the number of BFGS iterations required to reach a minimum Hamiltonian loss of at least $C_L^{min}={10}^{-12}$ for a given QLSP instance. The last two strategies seemingly provide the same level of acceleration with speedup factor over the original algorithm being $S\approx 2.9$ based on the median values. However, the MLP model didn’t outperform the naive Nearby method, since the mean square error, after an initial drop, usually was trapped in plateaus which were difficult to overcome without hyperparameter fine-tuning. Due to the fact that our experiments were simulated and did not include QPU runtimes, which are directly comparable to the time spend in fine-tuning, we were forced to be quite conservative with the latter. To our knowledge, unless one can use real quantum hardware and time metrics, it is not apparent how to take into account and compare the required MLP training time with respect to overall VQLS iterations. The assumption eventually made was that by using scalable models and as $N\to \infty$, which is the field quantum algorithms are supposed to excel, the model training time will be negligible with respect to the total solution time since $t_{\mathrm{train}}/t_{\mathrm{total}}\to 0$.

Regarding the low MLP performance, it seems that in order to outperform a naive nearby approach in such a small scale, a very delicate fine-tuning of the undermine the predictive ability of simple models played the most important role. Furthermore, the majority of the optimized parameters dataset ${\mathit{\alpha }}^i$ cannot be used for training, as is common practice, since increasing the training data means that more solutions will be generated using the base VQLS scheme which is slower.

6. Conclusions and future work

Variational algorithms — specifically in our context VQLS algorithm — are designed for the noisy qubits of the current NISQ era. In order to solve the QLSP, VQLS optimizes classically a parametric ansatz w.r.t a quantumly computed Hamiltonian cost function.

Our proposed machine-learning scheme aims to accelerate the convergence of variational quantum algorithms when they are utilized as solvers in parametric problems by taking advantage of continuity. As a proof of concept, we showed that there is potential in that direction, since we managed to achieve a speed-up factor of up to 2.9. In the studied cases, the nearby shceme which combines both speed-up and simplicity seems to be preferable. We also demonstrated that when using gradient-free solvers, and especially ones which favour explorability, discontinuities in the optimal ansatz parameters emerge for various reasons, most prominently due to non-convexity of Hamiltonian loss and periodicity of quantum ansatzes w.r.t. their paramaters. This considerably handicaps the performance of machine learning models such as MLPs, a problem that should be addressed in the future in order to take full advantage of modern ML models.

A natural next step in order to overcome those hindrances is to use exact gradient-informed optimizers, which shall favour exploitation over exploration of parameter-space. Additionally, using domain constraint optimizers will help alleviate discontinuities related to the periodicity of the Hamiltonian. In particular, the use of natural gradient methods initially proposed in [25] for classical neural networks and their quantum counterparts [26] may provide sufficient robustness to the optimation process and in this way addressing the discontinuities originating from non-convex loss functions. An extensive evaluation of transformation-preprocessing techniques may also help significantly in this regard. Without doubt, increasing the scale of the experiments will allow bigger training datasets, which will in turn enable more complex, deeper machine learning models to be trained. Finally, experiments with real quantum hardware noise must be conducted so that similar obstructive phenomena are highlighted and addressed.