Variational Quantum Algorithm for Solving Second-Order Linear Differential Equations

Sugaya, Ryo; Sogabe, Tomohiro; Kemmochi, Tomoya; Zhang, Shao-Liang

doi:10.2478/qic-2025-0012

Full Article

1.

Introduction

Solving partial differential equations (PDEs) is key to describing various issues in condensed matter physics, fluid simulation, electromagnetic field analysis, and other fields as mathematical models [1]. Therefore, identifying efficient algorithms for solving PDEs is crucial. Thus, numerous studies has explored the numerical analysis of PDEs.

Quantum algorithms (QAs), which are computational algorithms employing quantum computers, have attracted attention owing to their ability to solve specific issues more efficiently compared with classical computers. Notably, QAs have been demonstrated to outperform classical computers in solving certain problems. Hence, QAs have been deployed to solve problems in various fields such as machine learning [2,3], integer factorization [4], and combinatorial optimization [5]. Further, the Harrow–Hassidim–Lloyd (HHL) algorithm [6], a QA for solving linear equations, can solve linear equations obtained via PDEs discretization. However, the HHL algorithm requires a sufficient number of qubits and error correction techniques on quantum computers, which may require a long development time. To resolve this shortcoming, quantum-classical hybrid algorithms, employing noisy intermediate-scale quantum (NISQ) devices [7] with intermediate numbers of qubits, combined with classical computers, have gained attention. Among these quantum-classical hybrid algorithms, the variational quantum algorithms (VQAs) can operate on NISQ devices.

VQAs [8] are algorithms for minimizing objective functions containing approximate solution vectors. These approximate solution vectors are generated by quantum gates (ansatz gates) depending on parameters. The parameter updates are performed using iterative algorithms to minimize the objective function represented by approximate solution vectors. Thus, VQA is highly versatile and has been applied to various tasks such as solving eigenvalue problems using the variational quantum eigenvalue solver [9,10], quantum machine learning [2], and quantum approximate optimization algorithms [5].

To obtain the value of the objective function using VQA, the computation is performed by measuring the expectations. Furthermore, two conventional methods for calculating expectations on a quantum computer are projective measurements and the Hadamard test [11]. To compute expectations using projective measurement and the Hadamard test, the matrices involved must be measurement operators and unitary, respectively. Therefore, to implement VQAs, decomposing a matrix into a linear combination of unitary matrices and measurement operators represents a crucial challenge.

Many researchers have proposed methods for solving linear equations using VQAs [12–18]. Notably, Sato et al. proposed VQAs to solve linear equations obtained by the discretization of the Poisson equation [18]. To solve Ax = b for a symmetric positive definite matrix A, the minimizing problem for x^†Ax/2 – x^†f was solved using VQAs in the previous study, which was aimed at solving the Poisson equation, where x^†denotes the conjugate transpose of x. Furthermore, the number of terms involving expectations in the objective function value was accelerated to O(1) regardless of the matrix size by decomposing the coefficient matrix into a linear combination of unitary matrices and measurement operators. However, A must be a symmetric positive definite to ensure its formulation as a convex quadratic programming problem. Additionally, as A in the Poisson equation exhibits a specific structure, extending the methods employed in previous studies to other PDEs may be challenging. For example, for differential equations containing advection terms in fluid simulations, A becomes non-symmetric, thereby complicating the use of existing VQAs to solve the problem.

Thus, we propose VQAs for solving general second-order linear differential equations containing constant coefficients, e.g., the Poisson equation. Generally, second-order linear differential equations containing constant coefficients include advection terms, thus yielding non-symmetric coefficient matrices when discretized. For example, A obtained from the discretization of the Poisson equation with the Dirichlet boundary condition represents a special tridiagonal Toeplitz matrix that is symmetric and positive definite. Conversely the A obtained by discretizing a second-order linear differential equation containing constant coefficients and the Dirichlet boundary condition is a general tridiagonal Toeplitz matrix. To handle the non-symmetric A, we considered the normal equation A^†Ax = A^†b for the original linear equation (Ax = b) because A^†A is always a symmetric positive definite when A is a nonsingular, making it reducible to a convex quadratic programming problem¹. Furthermore, solving the normal equation is equivalent to solving the least squares problem for minimizing the norm ∥Ax – b∥. Similar to previous study [18], as matrices A^†A and A in the normal equation are not unitary matrices, they must be decomposed into linear combinations of unitary matrices and measurement operators. Thus, for A in second-order linear differential equations containing constant coefficients with the periodic boundary condition, we decomposed them into a linear combination of unitary matrices and measurement operators. Furthermore, by adding a correction term to A with the periodic boundary condition, we attempted the decomposition of the A for the case of the Dirichlet boundary condition. However, as the matrices employed for this correction do not possess unitarity and Hermiticity, we proposed a method for converting these correction matrices into projection operators or unitary matrices, ensuring it their applicability even under the Dirichlet boundary condition.

The remainder of this paper is structured as follows: Section 2 presents an overview of the extant studies on VQAs for solving the Poisson equation. Next, Section 3 reports the construction the objective function for the constantcoefficient second-order linear differential equations on the basis of VQAs. Thereafter, we demonstrated the decomposition method for the matrices included in the objective function under the periodic boundary condition. Next, we present the decomposition method for the correction matrices for cases involving the imposition of the Dirichlet boundary condition is imposed on section. Section 4 presents numerical examples using the proposed method implemented by Qiskit in this section. Finally, Section 5 concludes the paper and discussing the future directions.

2.

Variational Quantum Algorithm for the Poisson Equation

Here, we describe VQA for solving the Poisson equation based on Sato et al.’s study [18]. For the domain Ω = (L, R)^d, (L, R ∈ ℝ, L < R), the d-dimensional Poisson equation with a given function f : Ω → ℝ is given by the following equation: 1 $- Δ u (x) = f (x) x \in Ω .$

To simply explain the previous study, we consider here the Dirichlet boundary condition: $u (x) = 0 x \in \partial Ω .$

For simplicity, we restrict the study to the case where d = 1. Thus, by discretizing Eq. (1) using the finite difference method with the interval-partition width h = (R – L)/(N + 1), we derive the following linear equation: 2 $A u = f,$ where the matrix A ∈ ℝ^N×N is a symmetric positive-definite matrix, and u, f ∈ ℝ^N. Therefore, solving the differential equation Eq. (1) reduces to solving Eq. (2). Furthermore, as matrix A is positive definite and symmetric, solving the linear equation Eq. (2) is equivalent to determining the minimization solution ν_min of the following objective function E(ν): $E (v) : = \frac{1}{2} v^{†} A v - v^{†} f,$ where ν and f are assumed to be real vectors, and thus ν^†f = f^†ν.

In the following sections of this paper, we will consider N = 2ⁿ with respect to the number of qubits n ∈ ℤ⁺.

To minimize E_v(ν), the vector ν is assumed to be represented as ν = r |ψ(θ)⟩ and determine ν_min from this representation. Here, r ∈ ℝ is the norm of ν, and |ψ(θ)⟩ is a quantum state vector depending on the parameter θ. In this study, |ψ(θ)⟩ generated by unitary matrices U(θ) is conventionally referred to as an ansatz gate, which is constructed by stacking the layers of R_y gates and CNOT gates. Notably, the ansatz parameters θ is given by the following $θ = [^{θ_{1}^{1}, \dots, θ_{1}^{2 m}, θ_{2}^{1}, \dots, θ_{2}^{2 m}, \dots, θ_{n}^{1}, \dots, θ_{n}^{2 m}] ⊤},$ where the dimension of the parameter vector θ is linear with respect to the number of R_y gate m and number of qubits n.

On the basis of the foregoing, using r and θ, the objective function E_v (ν) is redefined as follows: 3 $E (v) = ℱ (r, θ) = \frac{1}{2} r^{2} ψ (θ) | A | ψ (θ) - 〈 r ψ (θ) | f 〉 .$

By minimizing ℱ(r, θ) with respect to r, we can reduce the problem to a minimization problem with respect to θ. $E (θ) = - \frac{1}{2} \frac{| ψ (θ) | f |^{2}}{〈 ψ (θ) | A | ψ (θ) 〉} .$

Furthermore, for any quantum state vector |φ⟩ or |ψ⟩, we define 4 $| ϕ, ψ 〉 : = \frac{1}{\sqrt{2}} (| 0 〉 | ϕ 〉 + | 1 〉 | ψ 〉) .$

Then, the following can be obtained using Eq. (4): 5 $E (θ) = - \frac{1}{2} \frac{{[〈 ψ (θ), f | X \otimes I^{\otimes n} | ψ (θ), f 〉]}^{2}}{〈 ψ (θ) | A | ψ (θ) 〉} .$

Furthermore, the quantum state |ϕ, ψ⟩ shown in Eq. (4) is constructed by the quantum circuit in Figure 1. Here, the unitary matrices U_ϕ and U_ψ satisfy |φ⟩ = U_φ|0⟩^⊗n and |ϕ⟩ = U_ϕ | 0⟩^⊗n, respectively. The construction of multiple controlled-U_ϕ and controlled-U_ψ requires O(n²) quantum gates and a circuit depth of O(n) [20].

The coefficient matrix A^(d) in the d-dimensional Poisson equation is expressed using matrix A as follows: 6 $A^{(d)} = A \otimes I^{\otimes n (d - 1)} + I^{\otimes n} \otimes A \otimes I^{\otimes n (d - 2)} + \dots + I^{\otimes n (d - 1)} \otimes A .$

Following Eq. (6), assuming a problem can be formulated in one dimension, it can also be formulated in higher dimensions. Therefore, in this paper, we will discuss the one-dimensional VQA in this paper.

3.

Variational Quantum Algorithms for Second-Order Linear Differential Equations

3.1.

Formulation

We consider the following d-dimensional second-order linear PDEs: 7 $- a \nabla \cdot (\nabla u (x)) + b \cdot \nabla u (x) + c u (x) = f (x), x \in Ω,$ where a, c ∈ ℝ, b ∈ ℝ^N, Ω = (L, R)^d and f : Ω → ℝ. When a ≠ 0 and c = 0, Eq. (7) represents the steady advection-diffusion equation and when a ≠ 0, c = 0 and b = 0, Eq. (7) represents the d-dimension Poisson equation. We discretize this equation (7) using the central difference method. In one dimension, let the mesh size be h = (R – L)/(N + 1). Then, we obtain the following difference equation: $- a \frac{u_{i + 1} - 2 u_{i} + u_{i - 1}}{h^{2}} + b \frac{u_{i + 1} - u_{i - 1}}{2 h} + c u_{i} = f_{i}, 1 \leq i \leq N .$

This difference equation is represented by the following linear equation: 8 $M u = f,$ where M ∈ ℝ^N×N and u, f ∈ ℝ^N. The components of the coefficient matrix M in Eq. (8) depend on the boundary conditions of PDEs. Additionally, let M_p and M_d denote M with the periodic boundary condition and Dirichlet boundary condition, u₀ = u_N and u₀ = u_N = 0 respectively. Here, the matrices M_p and M_d are given as follows: $\begin{array}{l} M_{p} : = [\begin{matrix} α & γ & 0 & \dots & 0 & 0 & β \\ β & α & γ & \dots & 0 & 0 & 0 \\ 0 & β & α & \dots & 0 & 0 & 0 \\ ⋮ & ⋱ & ⋮ \\ 0 & 0 & 0 & \dots & α & γ & 0 \\ 0 & 0 & 0 & \dots & β & α & γ \\ γ & 0 & 0 & \dots & 0 & β & α \end{matrix}], \\ M_{d} : = [\begin{matrix} α & γ & 0 & \dots & 0 & 0 & 0 \\ β & α & γ & \dots & 0 & 0 & 0 \\ 0 & β & α & \dots & 0 & 0 & 0 \\ ⋮ & ⋱ & ⋮ \\ 0 & 0 & 0 & \dots & α & γ & 0 \\ 0 & 0 & 0 & \dots & β & α & γ \\ 0 & 0 & 0 & \dots & 0 & β & α \end{matrix}], \end{array}$ where Ω = (0, 1), a mesh size h = 1/(N + 1), α = 2a/h² + c, β = – b/(2h) – a/h², and γ = b/(2h) – a/h². Let us solve the linear equation in Eq. (8). Assuming the matrix M is symmetric positive definite, solving Eq. (8) will be equivalent to minimizing an objective function $ℱ (v) : = \frac{1}{2} v^{†} M v - v^{†} f .$

However, as M is not always a symmetric positive-definite matrix, VQA for the Poisson equation cannot be applied to the optimization problem. Hence, we solve the following normal equation: 9 $M^{†} M u = M^{†} f$ instead of Eq. (8). Assuming M is nonsingular, the coefficient matrix M^†M will be a symmetric positive-definite matrix. Thus, we define the objective function E(ν), as follows: 10 $ℱ_{n o r m a l} (v) : = \frac{1}{2} v^{†} M^{†} M v - v^{†} M^{†} f$ to solve Eq. (9). Further, we express Eq. (10) as a function that depends on θ, similar to the previous section. Thus, we obtain Eq. (11) 11 $r_{opt} = \frac{〈 ψ (θ) | M^{†} | f 〉}{〈 ψ (θ) | M^{†} M | ψ (θ) 〉} .$

Additionally, we obtain the following objective function ε(θ) by substituting Eq. (11) for Eq. (3) 12 $ℰ (θ) = - \frac{1}{2} \frac{{[〈 ψ (θ) | M^{†} | f 〉]}^{2}}{〈 ψ (θ) | M^{†} M | ψ (θ) 〉} .$

Thus, Algorithm 1 expresses how to solve second-order linear PDEs using a VQA. Here, the formulation of the gradient of the objective function ▽ε (θ) is discussed in the Appendix A.

Algorithm 1

VQAs for second-order linear PDEs

Require: Objective function ε and gradient ▽ε;

an initial ansatz parameter θ₀ ; a tolerant parameter of gradient norm ϵ

Ensure: Optimized ansatz parameter θ_opt and an optimized quantum state |ψ(θ_opt)⟩

1: while ∥▽ε (θ)∥ ≥ ϵ do

2: Encode |ψ(θ)⟩ on a quantum computer

3: Evaluate ε(θ) and ▽ε (θ) on a quantum computer

4: Update θ using a optimization method employing ε(θ) and ▽ε (θ) on a classical computer.

5: end while

3.2.

Evaluation of the Objective Function

As M is not unitary, we cannot directly calculate $〈 ψ (θ) | M_{p}^{†} M_{p} | ψ (θ) 〉$ and $〈 ψ (θ) | M_{p}^{†} | f 〉$ in Eq. (3) on a quantum computer. Therefore, we first decompose $M_{p}^{†} M_{p}$ and $M_{p}^{†}$ into unitary matrices. Mp is decomposed into 13 $M_{p} = α I^{\otimes n} + β P + γ P^{†},$ where P is a shift operator defined by $P : = \sum_{i = 0}^{N - 1} | (i + 1) \mod N 〉〈 i | = [\begin{array}{l} 0 & 0 & 0 & \dots & 0 & 0 & 1 \\ 1 & 0 & 0 & \dots & 0 & 0 & 0 \\ 0 & 1 & 0 & \dots & 0 & 0 & 0 \\ 0 & 0 & 1 & \dots & 0 & 0 & 0 \\ ⋮ & \dots & ⋮ \\ 0 & 0 & 0 & \dots & 1 & 0 & 0 \\ 0 & 0 & 0 & \dots & 0 & 1 & 0 \end{array}] .$

Thus, matrix M^†M is represented by 14 $\begin{array}{l} M_{p}^{†} M_{p} & = (α I^{\otimes n} + β P^{†} + γ P) (α I^{\otimes n} + β P + γ P^{†}) \\ = (α^{2} + β^{2} + γ^{2}) I^{\otimes n} + α (β + γ) (P^{+} + P) + β γ (P P + P^{†} P^{†}) . \end{array}$

Employing Eq. (13). Let us denote $\begin{array}{l} | ψ_{+} (θ) 〉 : = P | ψ (θ) 〉, \\ | ψ_{-} (θ) 〉 : = P^{†} | ψ (θ) 〉 \end{array}$ with respect to the quantum state |ψ(θ)⟩. Further, we introduce |ψ₊ (θ)⟩ and |ψ – (θ)⟩ newly defined, into the following expectation terms of observable M^†M: $\begin{array}{l} 〈 ψ (θ) | M_{p}^{†} M_{p} | ψ (θ) 〉 & = (α^{2} + β^{2} + γ^{2}) \\ + α (β + γ) 〈 ψ (θ) | (P + P^{†}) | ψ (θ) 〉 \\ + β γ 〈 ψ (θ) | (P P + P^{†} P^{†}) | ψ (θ) 〉 \end{array}$ respectively. The first term ψ(θ)| (P + P^†) |ψ(θ)⟩ can yield $\begin{array}{l} 〈 ψ (θ) | (P + P^{†}) | ψ (θ) 〉 & = 〈 ψ (θ) | P | ψ (θ) 〉 + 〈 ψ (θ) | P^{†} | ψ (θ) 〉 \\ = 〈 ψ (θ) ∣ ψ_{+} (θ) 〉 + 〈 ψ_{+} (θ) ∣ ψ (θ) 〉 \\ = 2 〈 ψ_{+} (θ), ψ (θ) | X \otimes I^{\otimes n} | ψ_{+} (θ), ψ (θ) 〉 . \end{array}$

The second term ψ(θ)| (PP + P^†P^†) |ψ(θ) in Eq. (13) can be evaluated as $\begin{array}{l} 〈 ψ (θ) | (P P + P^{†} P^{†}) | ψ (θ) 〉 & = 〈 ψ (θ) | P P | ψ (θ) 〉 + 〈 ψ (θ) | P^{†} P^{†} | ψ (θ) 〉 \\ = 〈 ψ_{-} (θ) ∣ ψ_{+} (θ) 〉 + 〈 ψ_{+} (θ) ∣ ψ_{-} (θ) 〉 \\ = 〈 ψ_{-} (θ) ∣ ψ_{+} (θ) 〉 + 〈 ψ_{+} (θ) ∣ ψ_{-} (θ) 〉 \\ = 2 〈 ψ_{+} (θ), ψ_{-} (θ) | X \otimes I^{\otimes n} | ψ_{+} (θ), ψ_{-} (θ) 〉 . \end{array}$

From the foregoing, the denominator in Eq. (3) can be decomposed to $\begin{array}{l} 〈 ψ (θ) | M_{p}^{^{†}} M_{p} | ψ (θ) 〉 & = (α^{2} + β^{2} + γ^{2}) I^{\otimes n} \\ + 2 α (β + γ) 〈 ψ_{+} (θ), ψ (θ) | X \otimes I^{\otimes n} | ψ_{+} (θ), ψ (θ) 〉 \\ + 2 β γ 〈 ψ_{+} (θ), ψ_{-} (θ) | X \otimes I^{\otimes n} | ψ_{+} (θ), ψ_{-} (θ) 〉, \end{array}$ and we can calculate Eq. (14) on a quantum computer, where each term in the equation is represented as the expectation of the Pauli matrices. Thereafter, we evaluate the numerator in Eq. (3). The numerator can be decomposed into $\begin{array}{l} 〈 ψ (θ) | M_{p}^{†} | f 〉 & = α 〈 ψ (θ) ∣ f 〉 + β 〈 ψ (θ) | P^{†} | f 〉 + γ 〈 ψ (θ) | P | f 〉 \\ = α 〈 ψ (θ) ∣ f 〉 + β 〈 ψ_{+} (θ) ∣ f 〉 + γ 〈 ψ - (θ) ∣ f 〉 \\ = α 〈 ψ (θ), f | X \otimes I^{\otimes n} | ψ (θ), f 〉 \\ + β 〈 ψ_{+} (θ), f | X \otimes I^{\otimes n} | ψ_{+} (θ), f 〉 \\ + γ 〈 ψ - (θ), f | X \otimes I^{\otimes n} | ψ - (θ), f 〉 . \end{array}$

Note that if |ψ(θ)⟩, |ψ₊ (θ)⟩, |ψ_–(θ)⟩ and |f⟩ are real-valued vector, then the inner products ψ(θ)|f⟩, ⟨ψ₊(θ)|f⟩, ⟨ψ_–(θ)|f⟩ are commutative.

3.3.

Dirichlet Boundary Condition

The proposed method can be applied to the Dirichlet boundary condition. Notably, M_d can also be decomposed into 15 $M_{d} = α I^{\otimes n} + β P + γ P^{†} - M_{B C},$ where the matrix M_BC is given as follows: $M_{B C} = [\begin{matrix} 0 & 0 & 0 & \dots & 0 & 0 & β \\ 0 & 0 & 0 & \dots & 0 & 0 & 0 \\ 0 & 0 & 0 & \dots & 0 & 0 & 0 \\ ⋮ \\ 0 & 0 & 0 & \dots & 0 & 0 & 0 \\ 0 & 0 & 0 & \dots & 0 & 0 & 0 \\ γ & 0 & 0 & \dots & 0 & 0 & 0 \end{matrix}] = β | 0 〉〈 N - 1 | + γ | N - 1 〉〈 0 | .$

Thus, $M_{d}^{†} M_{d}$ can be decomposed into 16 $\begin{array}{l} M_{d}^{†} M_{d} & = (α I^{\otimes n} + β P^{†} + γ P - M_{B C}^{†}) (α I^{\otimes n} + β P + γ P^{†} - M_{B C}) \\ = M_{p}^{†} M_{p} - M_{p}^{†} M_{B C} - M_{B C}^{†} M_{p} + M_{B C}^{†} M_{B C} . \end{array}$

To calculate the expectation of Eq. (16) on a quantum computer, we decompose the matrix M_BC into symmetric and non-symmetric matrices, ν and W, respectively. 17 $M_{B C} = \frac{β + γ}{2} V + \frac{β - γ}{2} W .$

Here, ν and W are given as follows: $\begin{array}{l} V = [\begin{matrix} 0 & 0 & 0 & \dots & 0 & 0 & 1 \\ 0 & 0 & 0 & \dots & 0 & 0 & 0 \\ 0 & 0 & 0 & \dots & 0 & 0 & 0 \\ ⋮ \\ 0 & 0 & 0 & \dots & 0 & 0 & 0 \\ 0 & 0 & 0 & \dots & 0 & 0 & 0 \\ 1 & 0 & 0 & \dots & 0 & 0 & 0 \end{matrix}] = | 0 〉〈 N - 1 | + | N - 1 〉〈 0 |, \\ W = [\begin{matrix} 0 & 0 & 0 & \dots & 0 & 0 & 1 \\ 0 & 0 & 0 & \dots & 0 & 0 & 0 \\ 0 & 0 & 0 & \dots & 0 & 0 & 0 \\ ⋮ \\ 0 & 0 & 0 & \dots & 0 & 0 & 0 \\ 0 & 0 & 0 & \dots & 0 & 0 & 0 \\ - 1 & 0 & 0 & \dots & 0 & 0 & 0 \end{matrix}] = | N - 1 〉〈 0 | - | 0 〉〈 N - 1 | . \end{array}$

To compute the numerator and denominator of ε(θ) (Eq. (12)) on a quantum computer, we introduce two decomposition methods for V and W.

3.3.1.

Computing the Numerator of the Objective Function

Employing Eq. (16), the numerator of Eq. (12) is decomposed into the following terms: 18 $〈 ψ (θ) | M_{d}^{†} M_{d} | ψ (θ) 〉 = 〈 ψ (θ) | M_{p}^{†} M_{p} | ψ (θ) 〉 - 〈 ψ (θ) | M_{p}^{†} M_{B C} + M_{B C}^{†} M_{p} | ψ (θ) 〉 + 〈 ψ (θ) | M_{B C}^{†} M_{B C} | ψ (θ) 〉 .$

The first term in Eq. (18) can be computed using the proposed method with the periodic boundary condition. Therefore, the second term in Eq. (18) is decomposed into 19 $\begin{array}{l} 〈 ψ (θ) | M_{p}^{†} M_{B C} † M_{B C}^{+} M_{p} | ψ (θ) 〉 \\ = \frac{β + γ}{2} 〈 ψ (θ) | M_{p}^{†} V + V^{†} M_{p} | ψ (θ) 〉 + \frac{β - γ}{2} 〈 ψ (θ) | M_{p}^{†} W + W^{†} M_{p} | ψ (θ) 〉 \\ ⊨ α (β + γ) 〈 ψ (θ) | V | ψ (θ) 〉 \\ + β \frac{β + γ}{2} 〈 ψ (θ) | P^{†} V + V^{†} P | ψ (θ) 〉 + β \frac{β - γ}{2} 〈 ψ (θ) | P^{†} W + W^{†} P | ψ (θ) 〉 \\ + γ \frac{β + γ}{2} 〈 ψ (θ) | P^{†} V + V^{†} P | ψ (θ) 〉 + γ \frac{β - γ}{2} 〈 ψ (θ) | P^{†} W + W^{†} P | ψ (θ) 〉 . \end{array}$

We perform the decomposition using Eq. (17). Afterward, we extract the following outer products with a computational basis from the matrix included in each term in Eq. (19): $\begin{array}{l} P^{†} V + V^{†} P & = 2 | N - 1 〉〈 N - 1 | + | N - 2 〉〈 0 | + | 0 〉〈 N - 2 |, \\ P V + V^{†} P^{†} & = 2 | 0 〉〈 0 | + | 1 〉〈 N - 1 | + | N - 1 〉〈 1 |, \\ P^{†} W + W^{†} P & = | N - 2 〉〈 0 | + | 0 〉〈 N - 2 | - 2 | N - 1 〉〈 N - 1 |, \\ P W + W^{†} P^{†} & = | 1 〉〈 N - 1 | + | N - 1 〉〈 1 | - 2 | 0 〉〈 0 | . \end{array}$

To compute the expectations in a quantum computer, these outer products must be denoted as measurement operators. Further, by performing matrix decomposition using measurement operators, we localize the quantum circuits and reduce the number of measurements involved in expectation-value measurement.

Here, we define the measurement operators |GHZ⟩, |GHZ₊⟩, and |GHZ_–⟩, as 20 $| G H Z 〉 : = \frac{| 0 〉 + | N - 1 〉}{\sqrt{2}} = U_{G H Z} | 0 〉,$ 21 $| G H Z_{+} 〉 : = \frac{| 1 〉 + | N - 1 〉}{\sqrt{2}} = U_{G H Z_{+}} | 0 〉,$ 22 $| G H Z_{-} 〉 : = \frac{| 0 〉 + | N - 2 〉}{\sqrt{2}} = U_{G H Z_{-}} | 0 〉 .$

Here these unitary matrices, U_GHZ, U_GHZ₊, U_{GHZ_–} ∈ ℝ^N×N are equivalent to the quantum circuits shown in Figures 2, 3 and 4, respectively: these quantum circuits are constructed using 𝒪(n) fundamental quantum gates.

Furthermore, these outer products can be introduced as projection operators using Eqs. (20) and (21) and (22). $\begin{array}{l} | N - 1 〉〈 0 | + | 0 〉〈 N - 1 | & = 2 | G H Z 〉〈 G H Z | - | 0 〉〈 0 | - | N - 1 〉〈 N - 1 |, \\ | N - 1 〉〈 1 | + | 1 〉〈 N - 1 | & = 2 | G H Z_{+} 〉〈 G H Z_{+} | - | 1 〉〈 1 | - | N - 1 〉〈 N - 1 |, \\ | N - 2 〉〈 0 | + | 0 〉〈 N - 2 | & = 2 | G H Z_{-} 〉〈 G H Z_{-} | - | 0 〉〈 0 | - | N - 2 〉〈 N - 2 | . \end{array}$

From the above, Eq. (19) is given as follows, and the expectation can be obtained from the projection measurement for each term. $\begin{array}{l} 〈 ψ (θ) | M_{p}^{†} M_{B C} + M_{B C}^{†} M_{p} | ψ (θ) 〉 \\ = 2 α (β + γ) | 〈 ψ (θ) | U_{G H Z} | 0 〉 |^{2} + 2 β γ | 〈 ψ (θ) | U_{G H Z_{+}} | 0 〉 |^{2} \\ + 2 β γ | 〈 ψ (θ) | U_{G H Z_{-}} | 0 〉 |^{2} - {α (β + γ) + γ (β - 2 γ)} | 〈 ψ (θ) | 0 〉 |^{2} \\ - {α (β + γ) + β (γ - 2 β)} | 〈 ψ (θ) | N - 1 〉 |^{2} - β γ | 〈 ψ (θ) | 1 〉 |^{2} - β γ | 〈 ψ (θ) | N - 2 〉 |^{2} . \end{array}$

The third term of Eq. (18) is given as follows: $〈 ψ (θ) | M_{B C}^{†} M_{B C} | ψ (θ) 〉 = γ^{2} | 〈 ψ (θ) | 0 〉 |^{2} + β^{2} | 〈 ψ (θ) | N - 1 〉 |^{2} .$

From the above, Eq. (18) is decomposed into the following terms: $\begin{array}{l} 〈 ψ (θ) | M_{d}^{†} M_{d} | ψ (θ) 〉 = & (α^{2} + β^{2} + γ^{2}) + α (β + γ) 〈 ψ (θ) | (P + P^{†}) | ψ (θ) 〉 \\ + β γ 〈 ψ (θ) | (P P + P^{†} P^{†}) | ψ (θ) 〉 - 2 α (β + γ) | 〈 ψ (θ) | U_{G H Z} | 0 〉 |^{2} \\ - 2 β γ | 〈 ψ (θ) | U_{G H Z_{-}} | 0 〉 |^{2} - 2 β γ | 〈 ψ (θ) | U_{G H Z_{+}} | 0 〉 |^{2} \\ + {α (β + γ) + γ (β - γ)} | 〈 ψ (θ) ∣ 0 〉 |^{2} + β γ | 〈 ψ (θ) ∣ 1 〉 |^{2} \\ + β γ | 〈 ψ (θ) ∣ N - 2 〉 |^{2} + {α (β + γ) + β (γ - β)} | 〈 ψ (θ) ∣ N - 1 〉 |^{2} . \end{array}$

3.3.2.

Computing the Denominator of the Objective Function

Employing Eq. (15), we decomposed the denominator of Eq. (12) as follows: 23 $〈 ψ (θ) | M_{B C}^{†} | f 〉 = \frac{β + γ}{2} 〈 ψ (θ) | V | f 〉 + \frac{β - γ}{2} 〈 ψ (θ) | W | f 〉 .$

The given matrix factorization is limited to cases where V and W could be decomposed into measurement operators, represented as the expected values of V and W, respectively. Alternatively, we consider the decomposition method for calculation using a quantum computer. Furthermore, matrices V and W can be decomposed as follows: $\begin{array}{l} V & = \frac{V_{B C} + V_{B C}^{†}}{2}, \\ W & = \frac{- W_{B C} + W_{B C}^{†}}{2}, \end{array}$ where the unitary matrices V_BC and W_BC are defined $\begin{array}{l} V_{B C} : & = [\begin{matrix} 0 & 0 & 0 & \dots & 0 & 0 & \dots & 0 & 0 & 1 \\ 0 & 0 & 0 & \dots & 0 & 0 & \dots & 0 & 1 & 0 \\ 0 & 0 & 0 & \dots & 0 & 0 & \dots & 1 & 0 & 0 \\ ⋮ \\ 0 & 0 & 0 & \dots & 0 & 1 & \dots & 0 & 0 & 0 \\ 0 & 0 & 0 & \dots & - 1 & 0 & \dots & 0 & 0 & 0 \\ ⋮ \\ 0 & 0 & - 1 & \dots & 0 & 0 & \dots & 0 & 0 & 0 \\ 0 & - 1 & 0 & \dots & 0 & 0 & \dots & 0 & 0 & 0 \\ 1 & 0 & 0 & \dots & 0 & 0 & \dots & 0 & 0 & 0 \end{matrix}], \\ W_{B C} : & = [\begin{matrix} 0 & 0 & 0 & \dots & 0 & 0 & \dots & 0 & 0 & 1 \\ 0 & 0 & 0 & \dots & 0 & 0 & \dots & 0 & 1 & 0 \\ 0 & 0 & 0 & \dots & 0 & 0 & \dots & 1 & 0 & 0 \\ ⋮ \\ 0 & 0 & 0 & \dots & 0 & 1 & \dots & 0 & 0 & 0 \\ 0 & 0 & 0 & \dots & 1 & 0 & \dots & 0 & 0 & 0 \\ ⋮ \\ 0 & 0 & 1 & \dots & 0 & 0 & \dots & 0 & 0 & 0 \\ 0 & 1 & 0 & \dots & 0 & 0 & \dots & 0 & 0 & 0 \\ - 1 & 0 & 0 & \dots & 0 & 0 & \dots & 0 & 0 & 0 \end{matrix}] . \end{array}$

From the above, the expectation (23) is decomposed into the linear combination of the expectation of V_BC and W_BC. $\begin{array}{l} 〈 ψ (θ) | M_{B C}^{†} | f 〉 \\ = \frac{β + γ}{2} (〈 ψ (θ) | V_{B C} | f 〉 + 〈 ψ (θ) | V_{B C}^{†} | f 〉) + \frac{γ - β}{2} (〈 ψ (θ) | W_{B C} | f 〉 - 〈 ψ (θ) | W_{B C}^{†} | f 〉) \\ = \frac{β + γ}{4} (〈 ψ (θ), V_{B C} f | X \otimes I^{\otimes n} | ψ (θ), V_{B C} f 〉 + 〈 ψ (θ), V_{B C} f | X \otimes I^{\otimes n} | ψ (θ), V_{B C}^{†} f 〉) \\ + \frac{γ - β}{4} (〈 ψ (θ), W_{B C} f | X \otimes I^{\otimes n} | ψ (θ), W_{B C} f 〉 - 〈 ψ (θ), W_{B C}^{†} f | X \otimes I^{\otimes n} | ψ (θ), W_{B C}^{†} f 〉) . \end{array}$

Furthermore, the quantum circuits corresponding to V_BC and W_BC are shown in Figures 5 and 6, respectively. The quantum circuits of V_BC and W_BC comprise n – 1 control bits and one target multi-Toffoli gate. Since V_BC and W_BC include multi-controlled Toffoli gates, they require a circuit depth of O(n) and are composed of O(n) elementary quantum gates, with an ancilla qubit [21]. Furthermore, considering the circuit cost of Eq. (4), one can see that a total of O(n²) elementary quantum gates is required. Therefore, from the perspective of gate implementation costs [22], decomposing into projection operators is effective for obtaining the expected values of matrices $M_{B C}^{†}$ and $M_{B C}^{†}$ .

Classical shadows are a useful technique for efficiently estimating expectation values [23]. In particular, by using random Clifford circuits, it is possible to compute expectation values with respect to arbitrary matrices. However, since the upper bound of the variance for random Clifford measurements scales as O(tr(M²)), the variance depends on the values of α, β, and γ. That is, as these values increase, the variance also increases. In contrast, in our method, the error in the expectation value does not depend on the magnitude of the individual matrix components (α, β, and γ). Therefore, when tr(M) is large, our method can be considered advantageous

Additionally, expectation value estimation using the Extended Bell Measurement (XBM) is also a useful approach [24]. XBM is a technique for estimating expectation values of banded matrices and can compute expectation values involving tridiagonal matrices using O(n) measurement operators. The measurement operators used in XBM consist of at most n – 1 CNOT gates and one Hadamard gate. XBM is applicable to general problems including Eq. (7). However, simultaneous measurements using XBM yield not only the matrix elements required for computation but also additional and irrelevant elements, which may lead to a higher number of shots required to compute the expectation values. In contrast, our method, particularly with the periodic boundary condition, is tailored for second-order linear differential equations with constant coefficients, and therefore avoids computing irrelevant elements. We will investigate the relationship between the number of unnecessary elements and the required number of shots elsewhere in the future.

4.

Numerical Experiments

In this section, we present numerical-analysis experiments using the proposed method. In this experiment, a quantum computing simulation is performed using a classical computer. Further, Qiskit (version 0.38.0) is used as the quantum computation library [25], and Qiskit provided a simulator backend. In this experiment, we exclusively employed the Qiskit StatevectorSimulator. As StatevectorSimulator performed matrix-vector operations internally, it does not account for the sampling errors that may arise in the calculation of the expectation value.

The initial values of the variational parameters were randomly selected from the $0 \leq θ_{i}^{t} \leq 4 π$ range. Additionally, the BFGS method, a type of gradient-descent method, was employed as the parameter-update method.

Further, we constructed a specific ansatz gate layer (Figure 7) and the number of the ansatz-gate layer m was 3.

To compare the numerical solutions obtained with the proposed method, the solutions obtained using LU decomposition for the linear equation (Eq. (8)) were utilized.

The following PDEs was considered for the numerical calculations: 24 $- \frac{d^{2} u}{d x^{2}} + 0.2 \frac{d u}{d x} + 0.1 u = c_{f},$ where the boundary condition in Eq. (24) was the Dirichlet boundary condition. Furthermore, the right-hand vector is defined, as follows: $| f 〉 = 2^{n / 2} c_{f} (H^{\otimes n} | 0 〉^{\otimes n}),$ where the matrix H is the Hadamard gate and the coefficient c_f = 0.1.

Figure 8 shows that the solution obtained using the proposed method was plotted, exhibiting significant agreement with the analytical solution. Furthermore, Table 1 shows that regardless of the increase in the number of qubits, ϵ_c and ϵ_q exhibited nearly identical values, indicating that the solution obtained using the proposed method performed at a comparable level with that obtained via classical computation. Therefore, we concluded that the proposed method functioned appropriately as a method for solving second-order linear differential equations with constant coefficients.

Table 1.

Relative errors, ϵ_c and ϵ for n = 3, 4 and 5 qubit cases.

the number of qubits	the relative error ϵ_c	the relative error ϵ_q
3 qubits	7.8455 × 10⁻⁴	7.8478 × 10⁻⁴
4 qubits	2.2014 × 10⁻⁴	2.2027 × 10⁻⁴
5 qubits	5.8441 × 10⁻⁵	4.4997 × 10⁻⁴

5.

Conclusion

In this study, we proposed a VQA for solving constant-coefficient second-order linear differential equations with the periodic and the Dirichlet boundary condition, and our finding contributed the following two points to the extant literature:

We provided ε(θ) for the normal equations corresponding to the discretized linear equations and demonstrated a matrix decomposition method for reducing them to a computationally feasible form on a quantum computer.
For the case involving the imposition of the Dirichlet boundary condition, we proposed a method that focused on the difference between the matrices obtained with the periodic and Dirichlet boundary condition, following by decomposing this difference into a linear combination of unitary matrices and measurement operators.

To solve the linear equation, Mu = f, the condition number is given by $κ (M) = ‖ M ‖ \cdot ‖ M^{- 1} ‖ .$

When the condition number κ was large, the numerical solution tended to be unstable. Particularly, as κ of the normal equations M^†Mu = M^†f(κ(M^†M) was approximately the square of κ(M), directly solving the normal equations presented a numerically unstable problem. In classical computation, preconditioning is conventionally performed to improve κ. Therefore, QAs that improve κ of the coefficient matrices and reduce the computational cost of the objective function. However, studying the preconditioning matrices for efficient quantum computation is a challenging task, requiring entirely different approaches from the existing preconditioning methods [19]. In VQAs, It is important to avoid barren plateaus when computing the objective function is an important issue. Barren plateaus refer to the phenomenon where the gradient of the objective function vanishes exponentially with respect to the number of qubits. This issue arises due to the execution of deep quantum circuits or measurements involving global quantum operations [26]. Although a complete solution to the barren plateaus problem has not yet been discovered, strategies that mitigate its effects are preferable. For example, designing ansatz gates specialized for specific problems has been proposed as a potential approach [27].

Variational Quantum Algorithm for Solving Second-Order Linear Differential Equations

Full Article

Paradigm

My account