Comparative Analysis on Two Quantum Algorithms for Solving the Heat Equation

As of now, an optimal quantum algorithm solving partial differential equations eludes us. There are several different methods, each with their own strengths and weaknesses. In past years comparisons of these existing methods have been made, but new work has emerged since then. Therefore, we conducted a survey on quantum methods developed post-2020, applying two such solvers to the heat equation in one spatial dimension. By analyzing their performance (including the cost of classical extraction), we explore their precision and runtime efficiency advancements between the two, identifying advantages and considerations.

💡 Research Summary

**
This paper conducts a comparative study of two quantum algorithms developed after 2020 for solving the one‑dimensional heat equation ∂u/∂t = α ∂²u/∂x². The motivation is to revisit the earlier benchmark by Liu, Montanaro, and Somma (LMS20), which reported at most a quantum speed‑up in terms of precision, but omitted the cost of converting the quantum output back to a classical representation. By selecting two recent methods—(i) the discrete‑adiabatic‑theorem‑based quantum‑walk algorithm of Costa, An, Sanders et al. (CAS + 22) and (ii) the Taylor‑expansion with Gauss‑Lobatto‑Chebyshev discretization followed by Quantum Amplitude Estimation (QAEA) of Oz, San, and Kara (OSK 23)—the authors evaluate runtime scaling, condition‑number dependence, and practical overheads.

Algorithm 1 (CAS + 22)
The heat equation is discretized in space using central differences and in time with forward differences, yielding a block‑triangular linear system A ũ = b. The matrix A consists of identity blocks and a tridiagonal Laplacian L, with spectral norm ‖A‖ = Θ(1) and ‖A⁻¹‖ = Θ(m), where m is the number of spatial grid points; thus the condition number κ = Θ(m). After normalizing A to unit norm, a block‑encoding U_A is constructed, exploiting the sparsity of L. The quantum linear system solver then proceeds via a sequence of quantum‑walk steps that implement a discretized adiabatic evolution from an initial Hamiltonian H₀ (whose ground state encodes |b⟩) to a final Hamiltonian H₁ (whose ground state encodes the solution |x⟩ = A⁻¹|b⟩/‖A⁻¹|b‖). The adiabatic theorem guarantees that, provided the spectral gap Δ(s) stays sufficiently large, the error after T steps satisfies ε ≤ ‖∂H/∂s‖ T/Δ². Choosing T = O(κ log 1/ε) yields a total runtime O(κ log 1/ε) = O(m log 1/ε). To suppress residual high‑frequency diffusion modes, Chebyshev polynomial filtering is applied before measurement. Finally, the solution amplitude is extracted via a small number of measurements (log 1/ε) and converted to a classical vector using standard sampling techniques.

Algorithm 2 (OSK 23)
Instead of a uniform grid, this method employs Gauss‑Lobatto‑Chebyshev points for spatial discretization, which produce a more uniform eigenvalue distribution of the discrete Laplacian and reduce discretization error. Time integration is performed by a Taylor expansion of the ODE system derived from the spatial discretization, with truncation error O(Δt^p). The resulting integral representation of the solution is fed into a Quantum Amplitude Estimation routine. QAEA achieves an ε‑accurate estimate with O(log 1/ε) quantum queries, a quadratic improvement over naïve Monte‑Carlo sampling. The initial state |b⟩ is prepared via amplitude encoding; for smooth initial conditions this can be done in polylog (m) time. The condition number of the underlying linear system remains κ = Θ(m), so the overall complexity is O(m log 1/ε) plus polylogarithmic overhead for state preparation and QAEA circuit depth. However, the circuit implementing the Chebyshev points, the Taylor‑step unitaries, and the QAEA oracle is considerably deeper than that of the quantum‑walk approach, leading to larger constant factors.

Inclusion of Quantum‑to‑Classical Conversion
Both algorithms require measurement and classical post‑processing to obtain a usable temperature profile. The authors model this conversion as O(polylog m) operations, assuming efficient sampling and error mitigation. In practice, measurement noise, the need for repeated runs to achieve the desired confidence, and error‑correction overhead can dominate the total runtime on near‑term devices.

Comparison and Findings

• Both methods achieve the theoretically optimal logarithmic dependence on the target precision ε, confirming the promise of quantum speed‑up for high‑precision regimes.
• The runtime scales linearly with the spatial grid size m due to the condition number κ = Θ(m), identical to the earlier LMS20 analysis.
• CAS + 22’s quantum‑walk/adiabatic framework offers a conceptually cleaner implementation with relatively shallow circuits, but it relies on maintaining a sizable spectral gap throughout the evolution, which may be challenging for more complex PDEs.
• OSK 23’s Chebyshev discretization reduces discretization error and leverages QAEA’s logarithmic query complexity, yet the required oracle construction and deeper circuits increase practical overhead.
• When the cost of preparing |b⟩ and converting |x⟩ back to classical data is accounted for, the overall advantage over classical finite‑difference methods narrows, especially on Noisy Intermediate‑Scale Quantum (NISQ) hardware where circuit depth is limited.

Conclusion and Outlook
The study confirms that recent quantum algorithms can outperform classical solvers in the asymptotic regime of very high precision, but practical implementation hurdles—state preparation, circuit depth, error correction, and measurement overhead—remain significant. Future work should explore pre‑conditioning techniques to reduce κ, more efficient block‑encoding schemes, and hardware‑aware designs that minimize depth while preserving the logarithmic ε‑scaling. Only with such advances will quantum PDE solvers become competitive for realistic engineering problems.