A central test case for quantum linear system algorithms (QLSA) is elliptic PDEs after a finite element discretization. Most existing analyses focus on preparing a normalized solution state. But an end-to-end quantum PDE solver must also extract physical quantities of interest, such as fluxes, currents, tractions, and energy. These outputs require quantum measurement, and their observable norms may grow like $h^{-χ} $ with mesh size $h $, creating a readout bottleneck even when a quantum preconditioner reduces the condition-number dependence on $h$.
We present a multilevel framework for this readout problem, motivated by the variance-reduction mechanism of multilevel Monte Carlo (MLMC), which is naturally compatible with a multi-level finite element discretization. Instead of estimating the full fine-grid observable directly, the method estimates a telescoping sum of interlevel corrections, so that the fine-coarse cancellation is exposed before quantum measurement. Our algorithm is based on Schur-complement factorization of the corrected Green's operator through a Ritz-complement map. For quantities of interest with readout order $χ\leq 2$, the multilevel estimator removes the polynomial $h$-dependent readout overhead. With amplitude estimation, the remaining statistical dependence is $ \widetilde{O}(1/\varepsilon)$, i.e., Heisenberg scaling in the inference precision up to logarithmic factors and with direct sampling, the complexity is reduced to standard Monte Carlo scaling $\widetilde{O}(1/ε^2)$.
A central test case for quantum linear system algorithms (QLSA) is elliptic PDEs after a finite element discretization. Most existing analyses focus on preparing a normalized solution state. But an end-to-end quantum PDE solver must also extract physical quantities of interest, such as fluxes, currents, tractions, and energy. These outputs require quantum measurement, and their observable norms may grow like h−χ with mesh size h, creating a readout bottleneck even when a quantum preconditioner reduces the condition-number dependence on h.
We present a multilevel framework for this readout problem, motivated by the variance-reduction mechanism of multilevel Monte Carlo (MLMC), which is naturally compatible with a multi-level finite element discretization. Instead of estimating the full fine-grid observable directly, the method estimates a telescoping sum of interlevel corrections, so that the fine-coarse cancellation is exposed before quantum measurement. Our algorithm is based on Schur-complement factorization of the corrected Green's operator through a Ritz-complement map. For quantities of interest with readout order χ≤2, the multilevel estimator removes the polynomial h-dependent readout overhead. With amplitude estimation, the remaining statistical dependence is O(1/ε), i.e., Heisenberg scaling in the inference precision up to logarithmic factors and with direct sampling, the complexity is reduced to standard Monte Carlo scaling O(1/ε2).
Can a quantum PDE algorithm avoid the polynomial cost of resolving a fine spatial mesh? For standard fixed-order discretizations, direct classical methods require work polynomial in h−1, or equivalently in the number of spatial degrees of freedom Nh=Θ(h−d). Direct quantum implementations of a parabolic semigroup still have coherent complexity O(T/h), and gradient-dependent observables such as heat flux and dissipation introduce additional mesh dependence. Decay of the solution norm will further suppress the postselection probability for preparing a normalized final state.
We develop a multilevel quantum algorithm that estimates linear and quadratic observables directly and places the fine--coarse cancellation inside the circuit before measurement. A contour-based LCU reconstructs each target-time correction from a coherent family of shifted resolvent differences. Rather than block encoding the fine and coarse inverses separately, we encode their difference through a shifted Ritz--Schur factorization, exposing its O(hℓ2) two-grid normalization. For Fourier hierarchies, the corresponding SELECT oracle consists of a quantum Fourier or sine transform, a spectral-band selector, and reversible diagonal arithmetic. We also give a non-Fourier realization based on energy-orthogonal dyadic midpoint details in one dimension, together with structured tensor-product extensions under fixed-rank coefficient and access assumptions.
For readouts with derivative order 0≤χ≤2, optimized amplitude estimation removes all polynomial dependence on the finest mesh size. Under the stated access assumptions, both linear and quadratic observables can be estimated with complexity O(1+(Tε)−1), with only polylogarithmic dependence on h−1.
We present a quantum preconditioning framework for solving linear systems arising from a finite difference discretization of the Poisson equation. It is based on the combination of the Schrödingerization technique \cite{JLY22b,JLYPRL24} and the BPX multilevel preconditioner in order to achieve near-optimal complexity. The Schrödingerization technique transforms linear partial and ordinary differential equations into Schrödinger-type systems with unitary evolution in one higher dimension, making them suitable for quantum simulation. A key contribution is a structure-aware construction of the block-encoding for the symmetrically preconditioned matrix AS=S⊤AS, where A is the stiffness matrix and S encodes the BPX preconditioner in factored form. By establishing a novel commuting identity, we avoid the unfavorable normalization scaling that would otherwise arise from naive multiplication of block-encodings. This yields an exact block-encoding of AS with normalization O(d2(L+1)), where d is the spatial dimension and L is the number of levels. Combined with the Schrödingerization-based Hamiltonian simulation, the overall quantum algorithm achieves a query complexity of O(poly(d)ε−1polylog(ε−1)) for estimating linear functionals of the solution to a given tolerance ε.
We present a structure-aware construction of quantum BPX preconditioners for Q1 finite element discretizations of constant-coefficient reaction--diffusion equations on the unit cube, including the Poisson equation. Building on existing quantum finite element and finite difference methods, we express the preconditioned operator through explicit one-dimensional maps, tensor products, and norm-preserving interlevel transfers. Symmetric preconditioning gives multilevel coefficients, and its favorable conditioning alone does not ensure efficient recovery of a state in the original finite element nodal basis. We instead use the local L2-orthonormal basis underlying the construction as the output representation of the same finite element function. Its coefficient norm equals the continuous L2 norm. Assuming direct access to the preconditioned right-hand side, we prepare this orthonormal coefficient state with a represented function error at most ε. The expected query complexity is O(dL3∥u∥L2∥f∥L2log(2+εL∥f∥L2)), where u is the exact solution, f is the right-hand side, d is the spatial dimension, and L=log2(1/h) for mesh width h. The bound holds for 0<ε≤∥u∥L2/2 and depends only polylogarithmically on mesh refinement for fixed continuous data. We also analyze linear-functional estimation and illustrate the construction and continuous error bounds with a one-dimensional UnitaryLab experiment.
In numerical approaches to solving differential equations on a lattice, a representation of the derivative operator that correctly matches the continuum behaviour of functions of momentum up to the band limit must be non-local. We present the construction of efficient linear-combination-of-unitaries (LCU)-based block-encodings for the first-order derivative and Laplacian operators in the non-local N=2n-dimensional SLAC representation. We use state-preparation techniques designed for smoothly decaying functions to prepare the dense LCU amplitudes with high success probability and low gate cost. Furthermore, we demonstrate how Shannon wavelet transforms can be applied to these block-encodings to obtain multiscale representations of the SLAC derivative operators. We then show how to apply a diagonal preconditioner that reduces the condition number of these matrices in the multiscale wavelet basis to a small constant. This enables the solution of partial differential equations (PDEs) with SLAC-discretised derivative operators on a finite lattice using the quantum linear solving algorithm (QLSA). For a d-dimensional PDE, after projection away from the nullspace, the resulting quantum linear-system algorithm has overall gate complexity O(dn3α(k)log(1/ε)), where α(k) is the subnormalisation factor of the order-k SLAC block-encoding and ε denotes the algorithmic approximation error.
Partial differential equations are a promising application area for fault-tolerant quantum algorithms, but the subject lies between two communities with different languages: numerical analysis and quantum computation. These lecture notes provide a numerically grounded introduction for readers entering from either field. Block encoding is the organizing principle: once a discretized differential operator is embedded in a unitary, primitives such as quantum singular value transformation, Hamiltonian simulation, linear combinations of unitaries, amplitude amplification, and measurement can be assembled into algorithms for elliptic, hyperbolic, and parabolic PDEs. Each chapter begins with a standard finite difference or finite element discretization and follows the full pipeline from the continuous PDE to quantum encoding, transformation, and extraction of a quantity of interest. Particular attention is paid to the factors governing end-to-end performance, including discretization error, state preparation, normalization, postselection, and measurement cost. A final chapter introduces nonlinear problems through Carleman and Koopman-von Neumann linearizations. The aim is not a comprehensive survey or a claim of universal quantum advantage, but a mathematically transparent entry point and a shared vocabulary for researchers in both communities.
Quantum methods for solving partial differential equations (PDEs) generally fall into a few general categories, each leveraging different quantum mechanisms to achieve potential speedups. One well-known approach is the quantum linear system solver, epitomized by the Harrow-Hassidim-Lloyd (HHL) algorithm, which uses phase estimation and eigenvalue inversion to solve large, sparse linear systems arising from discretized linear PDEs [1]. This method can yield exponential improvements over classical solvers under certain conditions (mostly when the matrix is well-conditioned and sparse). More recent developments, such as Quantum Singular Value Transformation (QSVT), provide a framework for matrix inversion in a more flexible way, often improving on the cost and error dependence of earlier solvers [2]. Time-dependent PDEs can be tackled via Hamiltonian simulations. The Trotter-Suzuki method decomposes the exponential of a Hamiltonian into products of exponentials of simpler sub-Hamiltonians, allowing a piecewise approximation of continuous dynamics [3]. Continuous-time quantum walks, which treat a PDE operator as a Hamiltonian that governs the walk, similarly simulate processes like diffusion or wave propagation on discretized grids or graphs [4]. Beyond direct simulation, variational quantum algorithms have emerged as a promising route in the noisy intermediate-scale quantum era. These approaches encode the PDE solution as a parameterized quantum state and optimize the parameters via classical feedback, thereby minimizing a cost function tied to the PDE residual or energy functional [5]. Quantum annealing or adiabatic methods, meanwhile, can solve optimization problems derived from PDE constraints by adiabatically evolving a simple initial Hamiltonian to the final problem Hamiltonian and staying in the system's ground state [3]. In our study, we implement a new parameterized quantum model, specifically designed to solve large linear systems on quantum hardware. This model is less deep and computationally faster than the proposed one in [5], reducing the complexity and speeding up convergence times. Such improvements are particularly valuable in cases where classical algorithms may fail due to the problem's size or complexity. Ongoing research is focusing on generalized quantum models that can be adapted to practical PDE problems for solving large electromagnetic problems where classical computers can fail.