January 01, 1970
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^{-\chi}\) 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 \(\chi\leq 2\), the multilevel estimator removes the polynomial \(h\)-dependent readout overhead. With amplitude estimation, the remaining statistical dependence is \(\widetilde{O}(\epsilon^{-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 \(\widetilde{O}(\epsilon^{-2})\).
Quantum computing technology is progressing rapidly across several hardware platforms. Recent milestones include noisy intermediate-scale demonstrations of quantum utility on large superconducting processors, below-threshold logical memory experiments, and manufacturable photonic quantum-computing modules [1]–[3]. These developments make it increasingly important to identify scientific-computing tasks for which quantum algorithms can provide end-to-end gains, rather than speedups for isolated substeps.
Partial differential equations (PDEs) are a natural target. They are central to scientific and engineering computation, with great potential to broaden the scope of quantum computing applications. In particular, elliptic equations arise in mechanics, conductivity, electrostatics and diffusion. A representative model problem is \[-\nabla\cdot(a(x)\nabla u(x))=f(x) \quad\text{in }\Omega, \qquad u|_{\partial\Omega}=0.\] After finite element discretization, such a problem typically leads to a large sparse linear system. Quantum linear-system algorithms (QLSAs) and block-encoding methods [4]–[10] are the natural algebraic tools for discretized elliptic PDEs. Subsequent linear system algorithms have improved precision dependence, block-encoding normalization, pseudoinverse implementation, and condition-number dependence for important subclasses [6], [11]–[14]. They have been used in quantum algorithms for Poisson and elliptic equations, finite element methods, and high-precision PDE solvers [13], [15]–[17].
For PDEs, the dependence of the condition number on the grid size \(h\) is essential because it introduces a direct dependence on the size of the linear system. Leveraging the factorization of the stiffness matrix from a Galerkin projection, one can reduce the condition number dependence from \(h^{-2}\) to \(h^{-1}\) [13]. Meanwhile, classical multigrid and subspace-correction methods remove the \(h^{-2}\) conditioning barrier by exploiting the nested hierarchy of finite element spaces [18]–[22]. In the quantum setting, under suitable access assumptions, the BPX preconditioner [18] can be realized by directly block encoding the preconditioned finite element operator [23], [24], or by mapping iteration schemes to unitary dynamics [25], [26].
Quantum preconditioners give an efficient route to the solve stage and linear-functional estimations. An end-to-end quantum solver, however, proceeds in three stages [16]: \[\boldsymbol{P}repare,\qquad \boldsymbol{S}imulate,\qquad \boldsymbol{I}nfer.\] The preparation stage loads the finite element right-hand side. For structured data, Grover–Rudolph-type binary trees and related state-preparation schemes can be efficient [27], [28] in that the circuit depth depend logarithmically on \(h\). These ideas have been incorporated in prior quantum PDE solvers as well [13], [16], [23]. The simulation stage applies a QLSA, a block-encoded inverse, or a related dynamical primitive. The inference stage extracts the desired physical output. Similar end-to-end accounting issues arise in other quantum scientific-computing pipelines [16], [29]. This paper focuses on the integration of inference and simulation stage for elliptic PDEs to achieve optimal complexity scaling.
The inference stage has two distinct subtleties. First, a QLSA prepares a normalized solution state, not the finite element vector itself. The norm is a separate scalar that must be estimated [30]. Second, physically meaningful PDE outputs are usually scalar quantities of interest (QoIs), not the full solution vector. Examples include averaged fields, boundary fluxes, tractions, electric currents, and energies. After discretization, these are linear or quadratic forms such as \[Q_h^{\rm lin}=\boldsymbol{g}_h^\dagger\boldsymbol{x}_h, \qquad Q_h^{\rm quad}=\boldsymbol{x}_h^\dagger M_h\boldsymbol{x}_h .\] More importantly, even if the elliptic solve has been preconditioned, the readout observable itself may have block-encoding scale \(O(h^{-\chi})\). e.g., first-derivative readouts such as fluxes or currents typically have \(\chi=1\), and energy quadratic forms may have \(\chi=2\). Thus a single-level estimator following a QLSA step still pays \(\widetilde{O}(\epsilon^{-1}h^{-\chi})\) with amplitude estimation and \(\widetilde{O}(\epsilon^{-2}h^{-2\chi})\) with standard sampling. This is a readout bottleneck rather than a linear-solver conditioning bottleneck.
We address this bottleneck by combining the multilevel Monte Carlo principle, which is naturally compatible with multilevel finite element. Classical MLMC uses the telescoping identity for QoIs at \(L\) levels, \[Q_L=Q_0+\sum_{\ell=1}^{L}(Q_\ell-Q_{\ell-1})\] and couples the fine and coarse quantities so that the level differences are much smaller than the full fine-grid output [31], [32]. Here it is assumed that \(Q_0\) is at a very coarse level and its computational cost can be neglected. In quantum observable estimation, the analogous gain is a reduction in the block-encoding normalization of the operator being measured. Estimating \(Q_\ell\) and \(Q_{\ell-1}\) separately and subtracting the results classically does not reduce the measurement cost. Rather, the cancellation must be represented inside the operator, i.e., inside the quantum circuit.
The key novel approach in this paper is the corrected Green’s operator \[D_\ell = B_\ell-P_\ell B_{\ell-1}P_\ell^\dagger, \qquad B_\ell=A_\ell^{-1},\] which is the fine-grid inverse response after subtracting the lifted coarse Galerkin response via a prolongation. A direct subtraction of two inverse block encodings would not expose the cancellation. Instead, we prove the Ritz-complement Schur factorization \[D_\ell = J_\ell\Sigma_\ell^{-1}J_\ell^\dagger .\] Here \(J_\ell\) is a Ritz-complement map and \(\Sigma_\ell\) is the associated two-level Schur complement. If \(J_\ell\) has block-encoding normalization \(O(h_\ell)\) and \(\Sigma_\ell\) is uniformly conditioned, then \(D_\ell\) has normalization \(O(h_\ell^2)\). This two-power factor can be obtained for \(P_1/Q_1\) \(L^2\)-response setup of the finite element methods. The simulation and inference stages are therefore integrated: the quantum circuit measures a corrected inverse-response operator whose scale has already been reduced by the multilevel coupling. The procedure no longer requires a separate QLSA step.
The comparison of the complexity is summarized in 1. The table is stated for linear elements with standard \(O(h_\ell^2)\) \(L^2\) error and for \(0\le\chi\le2\). It compares scalar QoI estimation costs; logarithmic factors are hidden in \(\widetilde{O}(\cdot)\). The second row represents positive-definite or factorized QLSA methods that reduce the elliptic solve scale from \(h^{-2}\) to \(h^{-1}\) when their structural assumptions apply [13]. The third row represents the BPX quantum preconditioner approach [23], where the solve scale is removed but the observable scale remains, while the interface with quadratic QoIs is still open[23]. The fourth row is the Schur-complement multilevel estimator developed here, which under the Ritz-complement oracle assumptions and for observables up to the energy level, has Heisenberg scaling in the statistical precision, up to logarithmic factors, while the direct-sampling version has the usual shot-noise scaling. In both cases, the polynomial mesh-dependent observable overhead is removed entirely. Montanaro and Pallister already adopted an end-to-end perspective by accounting for the cost of estimating a linear functional of the finite element solution [16], and the \(\epsilon^{-1}\) dependence of amplitude estimation is optimal in the usual query model, but did not consider the scaling in the observables.
| Method | Amplitude estimation | Standard sampling |
|---|---|---|
| Direct stiffness inverse estimator | \(\widetilde{O}(\epsilon^{-1}h^{-2-\chi})\) | \(\widetilde{O}(\epsilon^{-2}h^{-2-2\chi})\) |
| Positive-definite/factorized QLSA [13] | \(\widetilde{O}(\epsilon^{-1}h^{-1-\chi})\) | \(\widetilde{O}(\epsilon^{-2}h^{-1-2\chi})\) |
| BPX linear functional estimator [23] | \(\widetilde{O}(\epsilon^{-1}h^{-\chi})\) | \(\widetilde{O}(\epsilon^{-2}h^{-2\chi})\) |
| Schur-complement, \(0\le\chi\le2\) (This work) | \(\widetilde{O}(\epsilon^{-1})\) | \(\widetilde{O}(\epsilon^{-2})\) |
The rest of the paper is organized as follows. In 2 we briefly review the elliptic Galerkin discretization, the reduced scaling used for amplitude encoding, and representative linear and quadratic QoIs, to get the problem ready for quantum algorithms. 3 collects the block-encoding, QLSA, norm-estimation, and observable-estimation primitives needed for the overall complexity estimates. In 4 we outline the MLMC allocation principle, redirected for quantum observable estimation. The formulation of the Schur complement approach for QoI estimations is presented in 5 and 5.4 discusses finite element constructions supporting the Ritz-complement access oracle.
All finite-dimensional vectors are denoted by bold symbols, such as \(\boldsymbol{b},\boldsymbol{x},\boldsymbol{g}\). Matrices are denoted by capital letters. We use \(\dagger\) for the conjugate transpose of vectors and matrices, even when all quantities are real. The Euclidean vector norm is \[\left\lVert \boldsymbol{x}\right\rVert_2=(\boldsymbol{x}^\dagger \boldsymbol{x})^{1/2}.\] The induced matrix norm, the spectral norm, of a matrix is denoted by \(\left\lVert A\right\rVert_2\). If \(A\) is Hermitian positive definite, then \[\kappa(A)=\left\lVert A\right\rVert_2\left\lVert A^{-1}\right\rVert_2 .\] For a rectangular matrix \(A \in \mathbb{C}^{m,n}\), \(m>n,\) with full rank, its condition number is defined by \[\kappa(A)=\frac{\sigma_{\max}(A)}{\sigma_{\min}(A)},\] where \(\sigma_{\max}(A)\) and \(\sigma_{\min}(A)\) are the largest and smallest nonzero singular values of \(A\).
The symbol \(\ell\) is reserved exclusively for discretization levels. We write \(V_\ell\) for the finite-dimensional space on level \(\ell\), \(h_\ell\) for the corresponding mesh width, and \(A_\ell,\boldsymbol{b}_\ell,\boldsymbol{x}_\ell\) for the associated discrete operator, right-hand side, and coefficient vector, respectively. The finest level is denoted by \(L\).
For a domain \(\Omega\subset \mathbb{R}^d\), we use the standard Hilbert spaces \[L^2(\Omega) = \left\{ u:\Omega\to \mathbb{R}: \int_\Omega |u(x)|^2\,dx < \infty \right\},\] with inner product \[(u,v)_{L^2(\Omega)} = \int_\Omega {u(x)}v(x)\,dx,\] and \[H^1(\Omega) = \left\{ u\in L^2(\Omega): \partial_{x_i}u\in L^2(\Omega),\;i=1,\dots,d \right\},\] with norm \[\left\lVert u\right\rVert_{H^1(\Omega)} = \left( \left\lVert u\right\rVert_{L^2(\Omega)}^2 + \left\lVert \nabla u\right\rVert_{L^2(\Omega)}^2 \right)^{1/2}.\] When homogeneous Dirichlet boundary conditions are imposed, the corresponding energy space is \(H_0^1(\Omega)\). In most applications, \(d=2\) or 3.
Bra-ket notation is reserved for normalized quantum states. Thus \(\left|v\right\rangle\) denotes \[\left|v\right\rangle = \frac{\boldsymbol{v}}{\left\lVert \boldsymbol{v}\right\rVert_2},\] unless \(\boldsymbol{v}\) is already normalized. Physical coefficient vectors such as \(\boldsymbol{x}_\ell\) are not identified with kets unless normalized explicitly.
An \(n\)-qubit quantum register has Hilbert-space dimension \(N=2^n\), with computational basis states \(\{\left|j\right\rangle:0\le j\le N-1\}.\) When a discrete PDE nodal coefficient vector \(\boldsymbol{x} \in \mathbb{C}^N\) is encoded as a quantum state, we use the amplitude encoding \[\label{q-encode} \boldsymbol{x} \longmapsto \left|x\right\rangle = \frac{1}{\left\lVert \boldsymbol{x}\right\rVert_2} \sum_{j=0}^{N-1} x_j \left|j\right\rangle.\tag{1}\]
If the number of degrees of freedom is not exactly a power of two, the vector is padded with zeros to dimension \(2^n\). This padding is understood throughout and does not change the Euclidean norm of the coefficient vector.
We write \(O(\cdot)\) and \(\widetilde{O}(\cdot)\) with constants independent of the mesh size \(h\). These constants may depend on the fixed spatial dimension \(d\) of the PDE, shape-regularity constants, and the ellipticity coefficient \(a(x)\). The notation \(\widetilde{O}(\cdot)\) suppresses polylogarithmic factors in the relevant parameters.
Let \(\Omega\subset\mathbb{R}^d\) be a bounded polyhedral domain. We consider the scalar elliptic model problem \[-\nabla\cdot(a(x)\nabla u(x))=f(x),\qquad u|_{\partial\Omega}=0, \label{eq:pde}\tag{2}\] where \(a(x)\in\mathbb{R}^{d\times d}\) is symmetric and uniformly elliptic \(\forall \boldsymbol{\xi} \in \mathbb{R}^d\): \[a_{\min}\left|\boldsymbol{\xi}\right|^2\le \boldsymbol{\xi}^\dagger a(x)\boldsymbol{\xi}\le a_{\max}\left|\boldsymbol{\xi}\right|^2 .\] The weak formulation corresponds to finding \(u\in H_0^1(\Omega)\) such that \[a(u,v):=\int_\Omega \nabla v(x)^\dagger a(x)\nabla u(x)\,\mathrm dx = \int_\Omega f(x)v(x)\,\mathrm dx=:F(v), \label{eq:weak}\tag{3}\] for all \(v\in H_0^1(\Omega)\). The Galerkin discretization and the finite element error estimates used below are standard; see, for example, [33]. The positive-definite quadratic form in 3 also induces an energy norm, \(\left\lVert u\right\rVert_a= \sqrt{a(u,u)},\) which is quite useful in the analysis.
The multi-level structure of the PDE discretization is constructed from a sequence of gradually refined spaces [34]: \[\label{nested} V_0\subset V_1\subset\cdots\subset V_L\subset H_0^1(\Omega).\tag{4}\] We consider standard P1 or Q1 finite element spaces on a shape-regular dyadic mesh hierarchy with \[h_\ell=h_{\ell-1}/2,\qquad h=h_L.\] The finest mesh width \(h_L\) is chosen so that the spatial discretization error is below the prescribed accuracy \(\epsilon\). For example, for a second-order elliptic discretization with error \(O(h_L^2)\), achieving accuracy \(\epsilon\) requires \(h_L=O(\sqrt{\epsilon})\). For a geometrically refined hierarchy with \(h_\ell \simeq 2^{-\ell}h_0\), the number of levels therefore satisfies \(L = O\!\left(\log_2 \frac{1}{h_L}\right) = O\!\left(\log_2 \frac{1}{\epsilon}\right).\) Thus the number of levels usually plays an insignificant role in the overall complexity.
Let \(N_\ell=\dim(V_\ell)\), and let \(\Phi_\ell:=\{\varphi_i^{(\ell)}\}_{i=1}^{N_\ell}\) be the nodal basis of \(V_\ell\). For convenience, we consider the standard Lagrange nodal basis, normalized by \[\varphi_i^{(\ell)}(x_j^{(\ell)})=\delta_{ij}. \label{eq:nodal-basis-property}\tag{5}\] The Galerkin solution expressed in \(V_\ell\) in terms of the shape functions \(\Phi_\ell\) has coefficient vector \(\boldsymbol{x}_\ell^{\rm raw}\in\mathbb{R}^{N_\ell}\) satisfying \[A_\ell^{\rm raw}\boldsymbol{x}_\ell^{\rm raw}=\boldsymbol{b}_\ell^{\rm raw}, \qquad (A_\ell^{\rm raw})_{ij}=a(\varphi_j^{(\ell)},\varphi_i^{(\ell)}), \qquad (\boldsymbol{b}_\ell^{\rm raw})_i=F(\varphi_i^{(\ell)}). \label{eq:discrete-system-raw}\tag{6}\] All matrices are real at this stage. When used in quantum circuits, they are viewed as complex matrices without changing the notation.
We keep the finite element formulation in nodal coordinates. No mass-matrix square root is formed. However, raw nodal \(\ell^2\)-norms differ from physical \(L^2\)-norms by a known grid-volume factor. On a quasi-uniform mesh, a finite element function \(v_\ell(x) \in V_\ell\) with nodal values \(\boldsymbol{v}_\ell^{\rm raw}\) has the scaling of \[\|v_\ell\|_{L^2(\Omega)}\simeq h_\ell^{d/2}\|\boldsymbol{v}_\ell^{\rm raw}\|_2, \label{eq:nodal-L2-equivalence}\tag{7}\] with constants independent of \(h_\ell\). Thus a nodal vector representing an \(O(1)\) physical \(L^2\)-function has raw Euclidean norm \(O(h_\ell^{-d/2})\). Meanwhile, the load vector has a dual scaling: If \(f\in L^2(\Omega)\), then \[(\boldsymbol{b}_\ell^{\rm raw})_i=\int_\Omega f(x)\varphi_i^{(\ell)}(x)\,dx,\] has the typical cell-volume scale \(h_\ell^d\), and the dual norm estimate gives \[\|\boldsymbol{b}_\ell^{\rm raw}\|_2 \lesssim h_\ell^{d/2}\|f\|_{L^2(\Omega)},\] for standard smooth or \(L^2\)-bounded data, up to mesh-independent constants.
Quantum amplitude encoding in 1 automatically normalizes such a vector in \(\ell^2\) norm. However, these scalar factors \(h_\ell^{d/2}\) are known and should be tracked explicitly, as follows \[\boldsymbol{x}_\ell:=h_\ell^{d/2}\boldsymbol{x}_\ell^{\rm raw}, \qquad \boldsymbol{b}_\ell:=h_\ell^{-d/2}\boldsymbol{b}_\ell^{\rm raw}, \qquad A_\ell:=h_\ell^{-d}A_\ell^{\rm raw}. \label{eq:response-coordinate-def}\tag{8}\] Then we arrive at the linear system of equations that we consider for the quantum implementation, \[A_\ell\boldsymbol{x}_\ell=\boldsymbol{b}_\ell. \label{eq:response-discrete-system}\tag{9}\] The point of 8 is only to remove artificial dimension-dependent powers from raw nodal \(\ell^2\)-norms.
This convention preserves the finite element energy pairing, which is commonly used in FEM analysis [33]. If we consider two functions in the subspace \(V_\ell\) \[u_\ell(x)= \Phi_\ell(x)^T\boldsymbol{x}_\ell^{\rm raw}, \qquad v_\ell(x)= \Phi_\ell(x)^T\boldsymbol{y}_\ell^{\rm raw},\] and further map them into the reduced coordinates, \(\boldsymbol{x}_\ell=h_\ell^{d/2}\boldsymbol{x}_\ell^{\rm raw}, \, \boldsymbol{y}_\ell=h_\ell^{d/2}\boldsymbol{y}_\ell^{\rm raw},\) then \[a(u_\ell,v_\ell) = (\boldsymbol{y}_\ell^{\rm raw})^\dagger A_\ell^{\rm raw}\boldsymbol{x}_\ell^{\rm raw} = \boldsymbol{y}_\ell^\dagger A_\ell\boldsymbol{x}_\ell . \label{eq:energy-preserved-amplitude-coordinates}\tag{10}\] Likewise, the load pairing is preserved: \[f(v_\ell) = (\boldsymbol{y}_\ell^{\rm raw})^\dagger \boldsymbol{b}_\ell^{\rm raw} = \boldsymbol{y}_\ell^\dagger \boldsymbol{b}_\ell .\] Hence the rescaled system 9 represents exactly the same Galerkin problem, but in coordinates whose Euclidean norm has the physical \(L^2\)-scale.
In accordance with the reduced scaling of the variables and operators, the multilevel transfer matrices are scaled consistently as well. If \[P_\ell^{\rm raw}:\mathbb{R}^{N_{\ell-1}}\to \mathbb{R}^{N_\ell}\] is the usual nodal prolongation between raw nodal vectors, then in reduced coordinates we set \[P_\ell := h_\ell^{d/2}P_\ell^{\rm raw}h_{\ell-1}^{-d/2}. \label{eq:reduced-prolongation}\tag{11}\] Thus, if \(\boldsymbol{y}_{\ell-1}\) is a reduced coordinate vector on level \(\ell-1\), then \(P_\ell\boldsymbol{y}_{\ell-1}\) is the reduced coordinate vector of the same finite element function represented on level \(\ell\): \[h_\ell^{-d/2} \Phi_\ell(x)^\dagger P_\ell\boldsymbol{y}_{\ell-1} = h_{\ell-1}^{-d/2} \Phi_{\ell-1}(x)^\dagger \boldsymbol{y}_{\ell-1}. \label{eq:P-function-meaning}\tag{12}\] When we write \(P_\ell V_{\ell-1}\), we mean this embedded coarse finite element subspace inside \(V_\ell\), or equivalently the coordinate image \(\operatorname{range}(P_\ell)\subset \mathbb{R}^{N_\ell}\).
After this subsection, unless explicitly stated otherwise, \[A_\ell,\quad \boldsymbol{x}_\ell,\quad \boldsymbol{b}_\ell,\quad P_\ell\] denote the reduced-coordinate objects. With this convention, the familiar Galerkin identities keep the same form: \[A_{\ell-1}=P_\ell^\dagger A_\ell P_\ell. \label{eq:reduced-galerkin-identities}\tag{13}\]
Another useful observation is that in this reduced coordinate and under standard assumptions [33], we have, \[\|A_\ell\|_2=O(h_\ell^{-2}). \label{eq:A-response-norm-scale}\tag{14}\]
In addition to the prolongation \(P_\ell\), we choose the restriction as the adjoint transfer \(P_\ell^\dagger\) as standard in Galerkin multigrid [19]–[21], [34]. Due to the Galerkin projection on the nested finite element spaces, the load vector and stiffness matrices are related between two successive levels of discretization, \[A_{\ell-1}=P_\ell^\dagger A_\ell P_\ell, \qquad \boldsymbol{b}_{\ell-1}=P_\ell^\dagger\boldsymbol{b}_\ell. \label{eq:galerkin-consistency}\tag{15}\] In addition, we define the Green’s function at the same level as \[B_\ell:=A_\ell^{-1}. \label{eq:B-def}\tag{16}\] Certainly we have \(\boldsymbol{x}_\ell=B_\ell\boldsymbol{b}_\ell.\)
At the deepest level we frequently abbreviate \[A:=A_L, \qquad B:=B_L, \qquad \boldsymbol{b}:=\boldsymbol{b}_L, \qquad \boldsymbol{x}:=\boldsymbol{x}_L,\] which are the quantities one would work with in a single-level method.
In applications, PDE simulations are rarely used only to recover all nodal values of the solution. More often, one computes scalar quantities of interest (QoIs), such as averages, fluxes, currents, tractions, compliances, or energies. We focus on QoIs that are linear or quadratic in the finite element solution and those that admit compatible level-wise representations on the nested spaces.
A linear QoI has the form \[\mathcal{L}_\ell(\boldsymbol{x}_\ell)=\boldsymbol{g}_\ell^\dagger\boldsymbol{x}_\ell = \boldsymbol{g}_\ell^\dagger B_\ell\boldsymbol{b}_\ell, \label{eq:linear-qoi}\tag{17}\] where \(\boldsymbol{g}_\ell\) is the reduced-coordinate vector representing the readout functional. Typical examples include weighted averages, \[\int_\Omega \psi(x)u_\ell(x)\,\mathrm dx,\] as well as fluxes, tractions, or currents through an interface \(\Gamma\), \[\int_\Gamma a(x)\nabla u_\ell(x)\cdot n(x)\,\,\mathrm ds.\] For volume averages with bounded weights, the readout vector is usually \(O(1)\) in the reduced scaling. Derivative-based readouts are different. If \(G_\ell \boldsymbol{x}_\ell\) denotes a discrete gradient or flux vector associated with \(u_\ell\), then standard finite element inverse estimates give \[\left\lVert G_\ell\right\rVert=O(h_\ell^{-1})\] in the reduced \(L^2\)-compatible scaling. Thus, when such QoI is written in the form 17 , the transpose of the derivative acts on \(\boldsymbol{g}_\ell\), and the readout norm may carry an extra \(h_\ell^{-1}\).
A quadratic QoI has the form \[\mathcal{Q}_\ell(\boldsymbol{x}_\ell)=\boldsymbol{x}_\ell^\dagger M_\ell\boldsymbol{x}_\ell = \boldsymbol{b}_\ell^\dagger B_\ell^\dagger M_\ell B_\ell\boldsymbol{b}_\ell, \label{eq:quadratic-qoi}\tag{18}\] for some Hermitian \(M_\ell.\) Examples include total power or \(L^2\)-mass, \[\int_\Omega \left|u_\ell(x)\right|^2\,\mathrm dx,\] and elastic energy, \[\int_\Omega \nabla u_\ell(x)^\dagger a(x)\nabla u_\ell(x)\,\mathrm dx.\] For former class of QoIs, \(M_\ell\) is typically \(O(1)\) in reduced coordinates. For energy-type observables, however, \(M_\ell\) contains two derivatives. For instance, it has the schematic form \[M_\ell \simeq G_\ell^\dagger U_{a} G_\ell,\] where \(U_{a}\) is a diagonal block encoding of the coefficient \(a(x)\). Since \(\left\lVert G_\ell\right\rVert=O(h_\ell^{-1})\), such quadratic observables may have \[\left\lVert M_\ell\right\rVert=O(h_\ell^{-2}).\]
Because the spaces are nested and the discretization is Galerkin-compatible, readout operators on consecutive levels are related by prolongation and restriction. For linear QoIs, compatibility means that evaluating the coarse readout on a coarse vector gives the same value as evaluating the fine readout on its prolongation. Hence \[\boldsymbol{g}_{\ell-1}=P_\ell^\dagger\boldsymbol{g}_\ell. \label{eq:g-compatible}\tag{19}\] For quadratic QoIs, the corresponding compatibility condition is \[M_{\ell-1}=P_\ell^\dagger M_\ell P_\ell. \label{eq:M-compatible}\tag{20}\] These relations will be leveraged to design more efficient estimators.
Definition 1 (Effective readout order). For a family of linear readout vectors \(\boldsymbol{g}_\ell\), define its effective readout order \(\mathrm{ord} _{\rm lin}(\boldsymbol{g})=\chi_g\) if \(\chi_g\ge0\) is the smallest exponent such that, after known scalar input normalizations are included, \[\|\boldsymbol{g}_\ell\|_2\; = O(h_\ell^{-\chi_g}). \label{eq:g-effective-scale}\tag{21}\] For a family of quadratic observables \(M_\ell\), define its effective readout order \(\mathrm{ord} _{\rm quad}(M)=\chi_{_M}\) if \(\chi_{_M}\ge0\) is the smallest exponent such that \[\left\lVert M_\ell\right\rVert = O(h_\ell^{-\chi_{_M}}). \label{eq:M-effective-scale}\tag{22}\]
The algebra above is written for unnormalized finite element coefficient vectors. Even after we have eliminated the \(d\)-dependent scaling with \(h_\ell\), the states \(\boldsymbol{b}_\ell\) and \(\boldsymbol{g}_\ell\) may still contain other constants. They contribute to the actual values of the QoIs. The quantum circuit prepares normalized states, \[|b_\ell\rangle:=\frac{\boldsymbol{b}_\ell}{\|\boldsymbol{b}_\ell\|_2}, \qquad |g_\ell\rangle:=\frac{\boldsymbol{g}_\ell}{\|\boldsymbol{g}_\ell\|_2}, \label{eq:normalized-bg-states}\tag{23}\] whenever the corresponding vector is nonzero. Therefore the physical linear QoI is recovered from a normalized overlap as \[\mathcal{L}_\ell = \|\boldsymbol{g}_\ell\|_2\,\|\boldsymbol{b}_\ell\|_2 \langle g_\ell|B_\ell|b_\ell\rangle . \label{eq:linear-normalized-prefactor}\tag{24}\] Likewise, the physical quadratic QoI is \[\mathcal{Q}_\ell = \|\boldsymbol{b}_\ell\|_2^2 \langle b_\ell|B_\ell^\dagger M_\ell B_\ell|b_\ell\rangle . \label{eq:quadratic-normalized-prefactor}\tag{25}\] The scalar factors \(\|\boldsymbol{b}_\ell\|_2\) and \(\|\boldsymbol{g}_\ell\|_2\) are not obtained by sampling the solution state. They are known from the input preparation procedure or are classically computable from the data model. In observable estimation, these factors can be kept outside the expectation or absorbed into the effective observable scale 22 . Another important point is that while we can assume that \(\left\lVert \boldsymbol{b}_\ell\right\rVert\) is accessible, the norm of \(\boldsymbol{x}_\ell\) usually has to be estimated separately from a QLSA approach [30].
For PDE data, the load vector is not an arbitrary table. Its entries are finite element integrals, \[(\boldsymbol{b}_\ell^{\rm raw})_j=\int_\Omega f(x)\varphi_j^{(\ell)}(x)\,\mathrm dx. \label{eq:load-integrals}\tag{26}\] We assume a compressed finite element load model: for each level \(\ell\), a reversible classical subroutine provides the prefix weights and phases needed by the Grover–Rudolph/Kaye–Mosca binary-tree construction, so that the normalized state \(|b_\ell\rangle\) can be prepared with \(O(\log N_\ell)\) prefix-weight calls [27], [28]. A concrete route for such state preparation in finite element methods is described in [23] and earlier works [13], [16]. For analytic, piecewise polynomial, tensor-product, localized, or otherwise compressed right-hand sides \(f\), these prefix weights can often be evaluated in time polylogarithmic in \(N_\ell\).
In this model, the norm \(\left\lVert \boldsymbol{b}_\ell\right\rVert_2\) is not an additional quantum-estimation task: it is the root weight of the same binary tree used to prepare \(|b_\ell\rangle\). The same convention applies to readout vectors such as \(\boldsymbol{g}_\ell\). Therefore the scalar prefactors in 24 and 25 are treated as classically available whenever the corresponding state-preparation oracles are available.
This section summarizes the block-encoding tools needed for the quantum algorithms. The purpose is not to reprove standard results, but to expose the basic building blocks and to make clear where normalization factors enter the final measurement cost. These normalizations are not cosmetic: in observable estimation, amplitude estimation has cost linear in the observable block-encoding normalization, while direct sampling has variance controlled by its square.
Definition 2 (Block encoding). Let \(M\in\mathbb{C}^{m\times n}\). A unitary \(U_M\) is an \((\alpha_M,b,\epsilon_M)\)-block encoding of \(M\) if \[\left\lVert M-\alpha_M(\left\langle 0\right|^{\otimes b}\otimes I_m)U_M (\left|0\right\rangle^{\otimes b}\otimes I_n) \right\rVert_2 \le \epsilon_M .\] Here \(b\) is the number of ancilla qubits. For rectangular matrices, the input and output registers are understood to be embedded into a common unitary space by padding if needed. The scalar \(\alpha_M\) is the block-encoding normalization. In the exact case it must satisfy \(\alpha_M\ge\left\lVert M\right\rVert_2\). Following [23], one may also track the subnormalization \[\widetilde{\alpha}_M:=\alpha_M/\left\lVert M\right\rVert_2,\] which becomes important when block encodings are composed.
Lemma 1 (Sparse block encodings). Let \(M\in\mathbb{C}^{m\times n}\) be \(s_r\)-row-sparse and \(s_c\)-column-sparse. Suppose sparse-access oracles return the nonzero positions and entries, and suppose \[|M_{ij}|\le M_{\max}.\] Then \(M\) admits a block encoding with normalization \[\alpha_M=O\!\left(\sqrt{s_rs_c}\,M_{\max}\right)\] using \(O(1)\) sparse-oracle calls and polylogarithmic additional gates, up to entry-precision factors.
The normalization in 1 is consistent with the spectral norm bound \[\left\lVert M\right\rVert_2 \le \sqrt{\left\lVert M\right\rVert_1\left\lVert M\right\rVert_\infty} \le \sqrt{s_c s_r}\,M_{\max}.\]
Lemma 2 (Block-encoding algebra). Let \(U_A\) and \(U_B\) block encode compatible matrices \(A\) and \(B\) with normalizations \(\alpha_A\) and \(\alpha_B\). Then:
\(A^\dagger\) is block encoded with normalization \(\alpha_A\).
\(AB\) is block encoded with normalization \(\alpha_A\alpha_B\).
\(\mu_AA+\mu_BB\) is block encoded by LCU with normalization \[|\mu_A|\alpha_A+|\mu_B|\alpha_B.\]
The block diagonal matrix \(A\oplus B\) can be block encoded with normalization \[\max\{\alpha_A,\alpha_B\}.\]
The product and LCU constructions appear as Lemmas 53 and 52 of Gilyén–Su–Low–Wiebe [7]. Deiml and Peterseim give a form adapted to finite element block encodings [23].
The primitives in (i) and (ii) can be used to use the gradient matrix to build a block-encoding of the stiffness matrix. The normalization factors \(\alpha\) will directly impact the cost of estimating QoIs.
Lemma 3 (Inverse and pseudoinverse by QSVT). Let \(U_A\) be an \((\alpha_A,a,\epsilon_A)\)-block encoding of \(A\). Assume the nonzero singular values of \(A\) lie in \[[\sigma_{\min},\sigma_{\max}], \qquad \kappa_A:=\sigma_{\max}/\sigma_{\min}.\] Then QSVT gives a block encoding of \(A^+\) with normalization \[O(\sigma_{\min}^{-1})\] and queries to \(U_A\) \[O\!\left( \frac{\alpha_A}{\sigma_{\min}}\log(1/\epsilon) \right),\] up to standard logarithmic factors and polynomial-approximation conventions. In the common case \(\alpha_A=\Theta(\sigma_{\max})\), this becomes \(O\!\left(\kappa_A\log(1/\epsilon)\right).\)
This follows from polynomial approximation of \(1/x\) and singular value transformation; see [7].
This lemma is particularly useful for solving linear system of equations \(A \boldsymbol{x} = \boldsymbol{b}.\) On the other hand, once the solution \(\boldsymbol{x}\) is obtained from \(B=A^{-1}\) applied to \(\boldsymbol{b}\), further steps are needed to estimate quantities of interest [35].
Lemma 4 (Estimating linear and quadratic forms). Let \(B\) have a block encoding with normalization \(\alpha_B\), and let \(M\) have a Hermitian block encoding with normalization \(\alpha_M\). Suppose \(\left|b\right\rangle\) and, for linear forms, \(\left|g\right\rangle\), can be prepared.
The normalized overlap \(\left\langle g\right|B\left|b\right\rangle\) can be estimated with amplitude-estimation cost proportional to \[\alpha_B/\epsilon.\] The final value \(\boldsymbol{g}^\dagger B\boldsymbol{b}\) is then obtained by multiplying by \(\left\lVert \boldsymbol{g}\right\rVert_2\left\lVert \boldsymbol{b}\right\rVert_2.\)
The normalized quadratic form \(\left\langle b\right|B^\dagger M B\left|b\right\rangle\) can be estimated with amplitude-estimation cost proportional to \[\alpha_B^2\alpha_M/\epsilon.\] The final value \(\boldsymbol{b}^\dagger B^\dagger M B\boldsymbol{b}\) is then obtained by multiplying by \(\left\lVert \boldsymbol{b}\right\rVert_2^2.\)
Recall that multigrid methods are often constructed from a hierarchy of nested subspaces \[V_0\subset V_1\subset\cdots\subset V_L, \qquad h_\ell=h_{\ell-1}/2, \qquad h=h_L.\] Such a hierarchy naturally gives a sequence of level-dependent outputs \(Q_\ell\). The fundamental multilevel identity is the telescoping sum \[Q_L=Q_0+\sum_{\ell=1}^L \Delta_\ell, \qquad \Delta_\ell:=Q_\ell-Q_{\ell-1}. \label{eq:ml-telescoping}\tag{27}\] Classical MLMC combines 27 with a coupling of \(Q_\ell\) and \(Q_{\ell-1}\), so that \(\Delta_\ell\) has much smaller variance than \(Q_\ell\) itself [31], [32].
The same principle is useful in quantum observable estimation, but the relevant scale is not only a classical variance. Quantum readout has two standard regimes. With direct projective measurement, a level-\(\ell\) observable with effective scale \(\alpha_\ell\) gives a bounded random variable whose variance is at most \(O(\alpha_\ell^2)\), which will require 4. Thus estimating a level contribution to RMS accuracy \(\epsilon_\ell\) by direct sampling costs \[\widetilde{O}\!\left( \frac{T_\ell\alpha_\ell^2}{\epsilon_\ell^2} \right),\] where \(T_\ell\) is the cost of one level-\(\ell\) circuit call. Amplitude estimation [35] improves the precision scaling to \[\widetilde{O}\!\left( \frac{T_\ell\alpha_\ell}{\epsilon_\ell}. \right)\] In both regimes, the quantum analogue of MLMC variance reduction is to make the level operator itself smaller before measurement.
Estimating \(Q_\ell\) and \(Q_{\ell-1}\) separately and subtracting the two estimates classically does not lower the quantum measurement cost, because each circuit still sees the full single-level observable scale. The cancellation must be represented inside the operator or state entering the quantum circuit. The rest of the paper implements this idea for elliptic finite element quantities of interest.
We use the following allocation estimate repeatedly. For notational convenience, set \(\Delta_0:=Q_0\).
Lemma 5 (MLMC allocation for quantum observable estimation). Let \[Q_L=\sum_{\ell=0}^L\Delta_\ell\] be a multilevel decomposition. Let \(T_\ell\) denote the cost of one level-\(\ell\) circuit call, including state preparation, block encodings, and any levelwise normalization subroutines used in the estimator. Let \(\alpha_\ell\) denote the effective block-encoding scale of the level-\(\ell\) observable after known scalar prefactors have been included.
Suppose amplitude estimation estimates \(\Delta_\ell\) to additive accuracy \(\epsilon_\ell\), with \[\sum_{\ell=0}^L\epsilon_\ell\le \epsilon,\] and with level cost \[\mathcal{C}_\ell^{\rm AE}(\epsilon_\ell) = \widetilde{O}\!\left(\frac{c_\ell}{\epsilon_\ell}\right).\] Then the optimized total cost satisfies \[\mathcal{C}_{\rm ML}^{\rm AE} = \widetilde{O}\!\left( \frac{1}{\epsilon} \left[\sum_{\ell=0}^L c_\ell^{1/2}\right]^2 \right). \label{eq:mlmc-ae-allocation}\tag{28}\] In particular, for a block-encoding observable estimator, \(c_\ell=T_\ell\alpha_\ell,\) and therefore \[\mathcal{C}_{\rm ML}^{\rm AE} = \widetilde{O}\!\left( \frac{1}{\epsilon} \left[\sum_{\ell=0}^L (T_\ell\alpha_\ell)^{1/2}\right]^2 \right). \label{eq:mlmc-ae-alpha}\tag{29}\]
Suppose direct sampling uses independent level estimators. Let \(v_\ell\) be the variance of one level-\(\ell\) sample, and let \(T_\ell\) be the cost of producing one such sample. Then the optimized cost to make the total sampling variance at most \(\epsilon^2\) is \[\mathcal{C}_{\rm ML}^{\rm samp} = \widetilde{O}\!\left( \frac{1}{\epsilon^2} \left[\sum_{\ell=0}^L (T_\ell v_\ell)^{1/2}\right]^2 \right). \label{eq:mlmc-samp-allocation}\tag{30}\] If one level-\(\ell\) sample is bounded in magnitude by \(O(\alpha_\ell)\), then \[v_\ell\le C\alpha_\ell^2\] for a mesh-independent constant \(C\), and hence the worst-case sampling cost is \[\mathcal{C}_{\rm ML}^{\rm samp} = \widetilde{O}\!\left( \frac{1}{\epsilon^2} \left[\sum_{\ell=0}^L T_\ell^{1/2}\alpha_\ell\right]^2 \right). \label{eq:mlmc-samp-alpha}\tag{31}\]
Proof. For amplitude estimation, we minimize \[\sum_{\ell=0}^L \frac{c_\ell}{\epsilon_\ell}, \quad \text{subject to} \; \sum_{\ell=0}^L\epsilon_\ell\le \epsilon.\]
The Lagrange multiplier equations give \[\epsilon_\ell = \epsilon \frac{c_\ell^{1/2}}{\sum_{j=0}^L c_j^{1/2}}, \Longrightarrow \sum_{\ell=0}^L\frac{c_\ell}{\epsilon_\ell} = \frac{1}{\epsilon} \left(\sum_{\ell=0}^L c_\ell^{1/2}\right)^2.\]
For block-encoding observable estimation, the amplitude estimation approach gives \(c_\ell=T_\ell\alpha_\ell\), which yields 29 .
For direct sampling, let \(N_\ell\) be the number of independent level-\(\ell\) samples. The total sampling variance and total cost are \[\sum_{\ell=0}^L \frac{v_\ell}{N_\ell}, \qquad \sum_{\ell=0}^L N_\ell T_\ell.\] Minimizing the cost subject to \[\sum_{\ell=0}^L \frac{v_\ell}{N_\ell}\le \epsilon^2\] gives \[N_\ell\propto \sqrt{\frac{v_\ell}{T_\ell}} \Rightarrow \frac{1}{\epsilon^2} \left(\sum_{\ell=0}^L (T_\ell v_\ell)^{1/2}\right)^2.\]
If the sample is bounded by \(O(\alpha_\ell)\), then \(v_\ell\le C\alpha_\ell^2\), giving 31 . ◻
For piecewise linear finite elements and smooth solutions, the natural \(L^2\)-response correction scale is \(O(h_\ell^2)\). Therefore, throughout the main results, we focus on the case where the level construction exposes a two-power factor \(h_\ell^2\) in the correction. This keeps the presentation centered on the standard \(P_1/Q_1\) finite element setting. Higher-order elements or rougher data can be treated by replacing the exponent \(2\) by the corresponding order.
This section gives an algebraic realization of the coupled level difference for 27 . It is based on a Green’s function correction, which we will express in terms of a Schur and Ritz complements.
Recall from 2.2 that \(B_\ell= A_\ell^{-1}\) is the Green’s function. We now define a lifted Green’s function from level \(\ell-1\) and consider its discrepancy with \(B_\ell,\) \[\label{Cl-Dl} C_\ell =P_\ell B_{\ell-1}P_\ell^\dagger, \qquad D_\ell=B_\ell- C_\ell ,\tag{32}\]
Lemma 6 (Level differences as one overlap or one expectation). Assume 15 , 19 , and 20 .
The linear level difference satisfies \[\Delta_\ell^{\rm lin} :=\boldsymbol{g}_\ell^\dagger B_\ell\boldsymbol{b}_\ell -\boldsymbol{g}_{\ell-1}^\dagger B_{\ell-1}\boldsymbol{b}_{\ell-1} = \boldsymbol{g}_\ell^\dagger D_\ell\boldsymbol{b}_\ell. \label{eq:D-linear-difference}\tag{33}\] Thus the level increment is a single overlap involving the residual-correction response.
The quadratic level difference satisfies \[\begin{align} \Delta_\ell^{\rm quad} &:= \boldsymbol{b}_\ell^\dagger B_\ell^\dagger M_\ell B_\ell\boldsymbol{b}_\ell - \boldsymbol{b}_{\ell-1}^\dagger B_{\ell-1}^\dagger M_{\ell-1}B_{\ell-1}\boldsymbol{b}_{\ell-1} \notag\\ &= \boldsymbol{b}_\ell^\dagger \left( D_\ell^\dagger M_\ell D_\ell +D_\ell^\dagger M_\ell C_\ell + C_\ell ^\dagger M_\ell D_\ell \right) \boldsymbol{b}_\ell. \label{eq:D-quadratic-difference} \end{align}\tag{34}\]
Proof. The linear identity follows from \[\boldsymbol{g}_{\ell-1}^\dagger B_{\ell-1}\boldsymbol{b}_{\ell-1} = \boldsymbol{g}_\ell^\dagger P_\ell B_{\ell-1}P_\ell^\dagger\boldsymbol{b}_\ell = \boldsymbol{g}_\ell^\dagger C_\ell \boldsymbol{b}_\ell.\] Subtracting this from \(\boldsymbol{g}_\ell^\dagger B_\ell\boldsymbol{b}_\ell\) gives 33 . For the quadratic identity, compatibility gives \[\boldsymbol{b}_{\ell-1}^\dagger B_{\ell-1}^\dagger M_{\ell-1}B_{\ell-1}\boldsymbol{b}_{\ell-1} = \boldsymbol{b}_\ell^\dagger C_\ell ^\dagger M_\ell C_\ell \boldsymbol{b}_\ell.\] Since \(B_\ell= C_\ell +D_\ell\), expanding \(B_\ell^\dagger M_\ell B_\ell- C_\ell ^\dagger M_\ell C_\ell\) gives 34 . ◻
The lemma expressed the level differences as an overlap or expectation through \(D_\ell\). Therefore, an efficient quantum estimation procedure in 4 can be applied to enable the MLMC speedup in 5, provided that \(D_\ell\) can be efficiently implemented via block encoding. We explain a possible path in the next section.
To construct a quantum algorithm for the corrected Green’s operator \(D_\ell\), we use the Schur complement associated with the nested spaces \(V_{\ell-1}\subset V_\ell .\) Throughout this section, all coefficient vectors and matrices are written in the reduced coordinates of 2.1. Throughout this section, \(P_\ell\) acts on reduced coordinate vectors; expressions such as \(P_\ell V_{\ell-1}\) refer to the embedded finite element functions described in 12 . In addition the Euclidean norm of a coefficient vector has a corresponding physical \(L^2\)-scale up to mesh-independent constants, while the finite element energy pairing is preserved.
Recall that \(\Phi_\ell(x) = \bigl(\varphi_1^{(\ell)}(x),\ldots,\varphi_{N_\ell}^{(\ell)}(x)\bigr)^T\) are the nodal basis vector on level \(\ell\). If \(\boldsymbol{x}_\ell\) is a reduced-coordinate vector, then the corresponding finite element function is \[u_\ell(x) = h_\ell^{-d/2} \Phi_\ell(x)^\dagger \boldsymbol{x}_\ell . \label{eq:function-from-reduced-vector}\tag{35}\] In particular we have established in 2.1 that, if \(u_\ell,v_\ell\in V_\ell\) are represented by reduced-coordinate vectors \(\boldsymbol{x}_\ell,\boldsymbol{y}_\ell\), then \[a(u_\ell,v_\ell) = \boldsymbol{y}_\ell^\dagger A_\ell \boldsymbol{x}_\ell . \label{eq:energy-pairing-reduced}\tag{36}\]
Let \(m_\ell:=N_\ell-N_{\ell-1},\) we choose a matrix \(W_\ell\in\mathbb{R}^{N_\ell\times m_\ell}\) such that \[\mathbb{R}^{N_\ell} = \operatorname{range}(P_\ell)\oplus\operatorname{range}(W_\ell), \label{eq:complement-splitting}\tag{37}\] and \(T_\ell:=[P_\ell\; W_\ell]\in\mathbb{R}^{N_\ell\times N_\ell}\) is invertible. The columns of \(W_\ell\) define a complement subspace in \(V_\ell\). In function form, this complement is represented by \[\Psi_\ell^\perp(x)^\dagger := h_\ell^{-d/2} \Phi_\ell(x)^\dagger W_\ell , \label{eq:W-function-basis}\tag{38}\] so that a complement function is \[w_\ell(x)=\Psi_\ell^\perp(x)^\dagger \boldsymbol{z}, \qquad \boldsymbol{z}\in\mathbb{R}^{m_\ell}.\] Thus \[V_\ell = P_\ell V_{\ell-1}\oplus \mathcal{W}_\ell, \qquad \mathcal{W}_\ell := \{\Psi_\ell^\perp(\cdot)^\dagger \boldsymbol{z}:\boldsymbol{z}\in\mathbb{R}^{m_\ell}\}.\] Such subspace decompositions and projection-based corrections are standard in subspace correction theory [19]. Ritz projections and Galerkin orthogonality are standard tools in finite element analysis [33]. The Schur complement algebra below is the same block-elimination principle [20], [21].
In the coordinate system determined by \(T_\ell\), every reduced-coordinate vector \(\boldsymbol{x}_\ell\in\mathbb{R}^{N_\ell}\) is written as \(\boldsymbol{x}_\ell=P_\ell\boldsymbol{y}_{\ell-1}+W_\ell\boldsymbol{z},\) or in function form, \[u_\ell(x) = \Psi_{\ell-1}(x)^\dagger \boldsymbol{y}_{\ell-1} + \Psi_\ell^\perp(x)^\dagger \boldsymbol{z}.\] The stiffness matrix in these new coordinates is \[\hat{A}_\ell := T_\ell^\dagger A_\ell T_\ell = \begin{bmatrix} A_{\ell-1} & E_\ell\\ E_\ell^\dagger & H_\ell \end{bmatrix}, \label{eq:T-coordinate-stiffness}\tag{39}\] where \[E_\ell:=P_\ell^\dagger A_\ell W_\ell, \qquad H_\ell:=W_\ell^\dagger A_\ell W_\ell . \label{eq:E-H}\tag{40}\] The matrix \(E_\ell\) is the energy coupling between the embedded coarse space \(P_\ell V_{\ell-1}\) and the chosen complement \(\mathcal{W}_\ell\).
The Ritz projection onto the embedded coarse space is the map \[\Pi_{\ell-1}^{a} := P_\ell A_{\ell-1}^{-1}P_\ell^\dagger A_\ell . \label{eq:ritz-projection-coordinate}\tag{41}\] It is characterized by \[a(v_\ell-\Pi_{\ell-1}^{a}v_\ell,\;P_\ell q_{\ell-1})=0 \qquad \forall q_{\ell-1}\in V_{\ell-1}.\] The corresponding Ritz complement is \(I-\Pi_{\ell-1}^{a}\). Applying this complement to the chosen complement basis \(W_\ell\) gives the Ritz-complement coordinate matrix \[J_\ell := (I-\Pi_{\ell-1}^{a})W_\ell = W_\ell-P_\ell A_{\ell-1}^{-1}E_\ell . \label{eq:J-as-ritz-complement}\tag{42}\]
The associated Ritz-complement basis functions are \[\Psi_\ell^{J}(x)^\dagger := h_\ell^{-d/2}\Phi_\ell(x)^\dagger J_\ell . \label{eq:J-function-basis}\tag{43}\] Thus, if \(\boldsymbol{z}\in\mathbb{R}^{m_\ell}\), the finite element function represented by the reduced-coordinate vector \(J_\ell\boldsymbol{z}\) is \[u_\ell(x)=\Psi_\ell^J(x)^\dagger\boldsymbol{z} .\]
Finally, define the Schur complement \[\Sigma_\ell := H_\ell-E_\ell^\dagger A_{\ell-1}^{-1}E_\ell . \label{eq:Sigma}\tag{44}\] The interpretation of \(J_\ell\) and \(\Sigma_\ell\) is simple: \(J_\ell\) represents the chosen complement after removal of its coarse Ritz projection, while \(\Sigma_\ell\) is the energy matrix on this Ritz-complement coordinate space.
Lemma 7. The Ritz-complement \(J_\ell\) satisfies \[P_\ell^\dagger A_\ell J_\ell=0, \label{eq:J-ritz-orthogonality}\tag{45}\] and \[\Sigma_\ell=J_\ell^\dagger A_\ell J_\ell . \label{eq:Sigma-JAJ}\tag{46}\]
Consequently, if \(u_\ell(x)=\Psi_\ell^J(x)^\dagger\boldsymbol{z}\) is the finite element function represented by the reduced-coordinate vector \(J_\ell\boldsymbol{z}\), then \[\boldsymbol{z}^\dagger\Sigma_\ell\boldsymbol{z} =a(u_\ell,u_\ell),\] so \(\Sigma_\ell\) is the energy Gram matrix in that subspace.
Proof. Using \(E_\ell=P_\ell^\dagger A_\ell W_\ell\) and \(A_{\ell-1}=P_\ell^\dagger A_\ell P_\ell\), \[P_\ell^\dagger A_\ell J_\ell = P_\ell^\dagger A_\ell W_\ell - P_\ell^\dagger A_\ell P_\ell A_{\ell-1}^{-1}E_\ell = E_\ell-A_{\ell-1}A_{\ell-1}^{-1}E_\ell=0.\] For the energy identity, expand \[J_\ell^\dagger A_\ell J_\ell = \left(W_\ell-P_\ell A_{\ell-1}^{-1}E_\ell\right)^\dagger A_\ell \left(W_\ell-P_\ell A_{\ell-1}^{-1}E_\ell\right).\] Using again \(P_\ell^\dagger A_\ell W_\ell=E_\ell\) and \(P_\ell^\dagger A_\ell P_\ell=A_{\ell-1}\), the expansion gives \[J_\ell^\dagger A_\ell J_\ell = H_\ell-E_\ell^\dagger A_{\ell-1}^{-1}E_\ell = \Sigma_\ell .\] ◻
Lemma 8 (Schur-complement factorization ). Assume the splitting 37 and the Galerkin consistency. Then \[D_\ell = B_\ell-P_\ell B_{\ell-1}P_\ell^\dagger = J_\ell\Sigma_\ell^{-1}J_\ell^\dagger . \label{eq:D-schur}\tag{47}\]
Proof. Since \(T_\ell\) is invertible, \(B_\ell=A_\ell^{-1} = T_\ell \hat{A}_\ell^{-1}T_\ell^\dagger .\) The block inverse formula for 39 , with leading block \(A_{\ell-1}\), gives \[\hat{A}_\ell^{-1} - \begin{bmatrix} A_{\ell-1}^{-1} & 0\\ 0 & 0 \end{bmatrix} = \begin{bmatrix} -A_{\ell-1}^{-1}E_\ell\\ I \end{bmatrix} \Sigma_\ell^{-1} \begin{bmatrix} -E_\ell^\dagger A_{\ell-1}^{-1} & I \end{bmatrix}.\] Multiplying by \(T_\ell=[P_\ell\;W_\ell]\) on the left and \(T_\ell^\dagger\) on the right yields \[\left(W_\ell-P_\ell A_{\ell-1}^{-1}E_\ell\right) \Sigma_\ell^{-1} \left(W_\ell-P_\ell A_{\ell-1}^{-1}E_\ell\right)^\dagger = J_\ell\Sigma_\ell^{-1}J_\ell^\dagger .\] The subtracted block is \[T_\ell \begin{bmatrix} A_{\ell-1}^{-1} & 0\\ 0 & 0 \end{bmatrix} T_\ell^\dagger = P_\ell A_{\ell-1}^{-1}P_\ell^\dagger = P_\ell B_{\ell-1}P_\ell^\dagger.\] Therefore the left-hand side is precisely \(B_\ell-P_\ell B_{\ell-1}P_\ell^\dagger\). ◻
The factorization 47 reduces the implementation of the corrected Green’s operator \(D_\ell = J_\ell\Sigma_\ell^{-1}J_\ell^\dagger\) to two objects: the Ritz-complement coordinate matrix \(J_\ell\) and the Schur complement \(\Sigma_\ell\).
The scaling of \(J_\ell\) and \(\Sigma_\ell\) depends on how the complement coordinates are normalized. A raw complement \(W_\ell\) may be easy to construct, but its coordinates need not be suitable for block encoding. For the main estimates we use coordinates in which the Ritz-complement functions are energy stable.
Assumption 1 (Energy-stable Ritz-complement). There exist constants \(0<c_0\le C_0<\infty\), independent of \(\ell\), such that for every \(\boldsymbol{z}\in\mathbb{R}^{m_\ell}\), \[c_0\|\boldsymbol{z}\|_2^2 \le \boldsymbol{z}^\dagger\Sigma_\ell\boldsymbol{z}= a\!\left( \Psi_\ell^J(\cdot)^\dagger\boldsymbol{z}, \Psi_\ell^J(\cdot)^\dagger\boldsymbol{z} \right) \le C_0\|\boldsymbol{z}\|_2^2 . \label{eq:stable-ritz-complement}\tag{48}\]
This is a normalization condition on the chosen complement coordinates after the coarse Ritz component has been removed. In particular, it implies that the Schur complement \(\Sigma_\ell\) is uniformly well conditioned. The condition is automatic if the Ritz-complement basis is chosen energy-orthonormally, and it is verified in the structured examples in the next section.
Lemma 9. Under 1 and the standard \(H^2\)-regularity. If \(u_\ell(x) = \Psi_\ell^J(x)^\dagger\boldsymbol{z}, \, \boldsymbol{z}\in\mathbb{R}^{m_\ell},\) then \[\|u_\ell\|_{L^2(\Omega)} \le C h_\ell \|\boldsymbol{z}\|_2 . \label{eq:J-L2-scale}\tag{49}\] Equivalently, in the reduced coordinates of 2.1, \[\|J_\ell\|_2=O(h_\ell). \label{eq:J-reduced-norm}\tag{50}\]
Proof. By 7, the function \(u_\ell\) is Ritz-orthogonal to the embedded coarse space: \[a(u_\ell,q_{\ell-1})=0 \qquad \forall q_{\ell-1}\in P_\ell V_{\ell-1}.\] Let \(\phi\in H_0^1(\Omega)\) solve the dual problem \[a(v,\phi)=(v,u_\ell)_{L^2(\Omega)} \qquad \forall v\in H_0^1(\Omega).\] Then \(\|u_\ell\|_{L^2(\Omega)}^2 = a(u_\ell,\phi).\) Let \(I_{\ell-1}\phi\in P_\ell V_{\ell-1}\) be a coarse interpolant. Ritz orthogonality gives \(a(u_\ell,I_{\ell-1}\phi)=0,\) and therefore \[\|u_\ell\|_{L^2}^2 = a(u_\ell,\phi-I_{\ell-1}\phi) \le \|u_\ell\|_a\, \|\phi-I_{\ell-1}\phi\|_a .\] The \(P_1/Q_1\) interpolation estimate and elliptic regularity imply \[\|\phi-I_{\ell-1}\phi\|_a \le C h_\ell \|\phi\|_{H^2(\Omega)} \le C h_\ell \|u_\ell\|_{L^2(\Omega)} .\] Canceling one factor of \(\|u_\ell\|_{L^2}\), we obtain \[\|u_\ell\|_{L^2} \le C h_\ell\|u_\ell\|_a .\] Finally, by 1, \[\|u_\ell\|_a^2 = a(u_\ell,u_\ell) = \boldsymbol{z}^\dagger\Sigma_\ell\boldsymbol{z} \le C_0\|\boldsymbol{z}\|_2^2 .\] This proves 49 . Since reduced-coordinate Euclidean norm is equivalent to the physical \(L^2\)-scale, 50 follows. The argument is the same Aubin–Nitsche duality mechanism used in the standard finite element \(L^2\)-error estimate. ◻
The analytic estimate 50 is necessary for the desired quantum normalization, but it is not by itself an implementation. A block encoding must realize this scale at the circuit level.
Assumption 2 (Ritz-complement oracle). For each level \(\ell\), the reduced-coordinate Ritz-complement matrix \(J_\ell\) admits a block encoding with normalization \[\alpha_{J_\ell}=O(h_\ell) \label{eq:J-oracle-normalization}\tag{51}\] and gate/query cost \(T_{J_\ell}\) that is polylogarithmic in \(h_\ell^{-1}\) and in the inverse target precision.
A direct LCU implementation of \(J_\ell = W_\ell-P_\ell A_{\ell-1}^{-1}E_\ell\) would generally normalize by the sizes of the two summands and would therefore have normalization \(O(1)\), not \(O(h_\ell)\). Thus 2 is a structural implementation assumption: the \(O(h_\ell)\) scale must be exposed by the construction of the complement, for example through spectral band filters, local energy-orthogonal coordinates, or localized Ritz correctors. We give concrete examples in which this assumption can be verified.
Theorem 1 (Block encoding of \(D_\ell\)). Assume [ass:stable-ritz-complement,ass:J-oracle], together with the standard block encoding normalization \(\alpha_{A_\ell}=O(h_\ell^{-2})\) for \(A_\ell\) with gate/query cost \(T_{A_\ell}\). Then \(\Sigma_\ell=J_\ell^\dagger A_\ell J_\ell\) has an \(O(1)\)-normalized block encoding with cost \[\widetilde{O}(T_{J_\ell}+T_{A_\ell}).\] Moreover, \(\Sigma_\ell^{-1}\) has an \(O(1)\)-normalized block encoding with the same cost up to polylogarithmic precision factors. Consequently, \(D_\ell=J_\ell\Sigma_\ell^{-1}J_\ell^\dagger\) has a block encoding with normalization \(\alpha_{D_\ell}=O(h_\ell^2)\) and cost \(\widetilde{O}(T_{J_\ell}+T_{A_\ell}).\)
Proof. By 1, \(c_0 I\preceq \Sigma_\ell\preceq C_0 I,\) so \(\kappa(\Sigma_\ell)\le C_0/c_0=O(1).\)
For the block encoding, the product rule from 2 gives \[\alpha_{\Sigma_\ell} \le \alpha_{J_\ell}^2\alpha_{A_\ell} = O(h_\ell^2)\,O(h_\ell^{-2}) = O(1).\] Since \(\Sigma_\ell\) is uniformly conditioned, QSVT inversion from 3 gives an \(O(1)\)-normalized block encoding of \(\Sigma_\ell^{-1}\) with only polylogarithmic precision overhead. Consequently, we arrive at a block encoding of \(D_\ell = J_\ell\Sigma_\ell^{-1}J_\ell^\dagger\) with normalization \[\alpha_{D_\ell} \le \alpha_{J_\ell}^2\,\alpha_{\Sigma^{-1}_\ell} = O(h_\ell^2).\] ◻
Remark 1 (First-order access to the same Schur complement). Instead of block encoding \(A_\ell\) directly with normalization \(O(h_\ell^{-2})\), one may use the elliptic factorization \[A_\ell=G_\ell^\dagger G_\ell, \qquad \alpha_{G_\ell}=O(h_\ell^{-1}),\] using the gradient of the shape functions. Then \(\Sigma_\ell = J_\ell^\dagger A_\ell J_\ell = (G_\ell J_\ell)^\dagger(G_\ell J_\ell).\) and \(\alpha_{G_\ell J_\ell } \le \alpha_{G_\ell}\alpha_{J_\ell} = O(h_\ell^{-1})O(h_\ell) = O(1).\) This is consistent with factorized positive-definite QLSA constructions and makes the first-order nature of the Schur complement explicit.
The only component left in 34 is the operator \(C_\ell\). We show that it can be constructed recursively once we have an efficient protocol for \(D_\ell\).
Lemma 10. Let \(P_{j\to \ell}:=P_\ell P_{\ell-1}\cdots P_{j+1}\) denote the reduced-coordinate prolongation from level \(j\) to level \(\ell\). Assume the coarsest inverse \(B_0\) has an \(O(1)\)-normalized block encoding. If, for each \(1\le j<\ell\), the corrected Green’s operator \(D_j\) has a block encoding with normalization \(\alpha_{D_j}=O(h_j^2),\) then the lifted coarse response \(C_\ell:=P_\ell B_{\ell-1}P_\ell^\dagger\) has an \(O(1)\)-normalized block encoding with polylogarithmic overhead in the number of levels and precision.
Proof. By recursively applying \(B_j=P_jB_{j-1}P_j^\dagger+D_j,\) we obtain \[C_\ell = P_{0\to \ell}B_0P_{0\to\ell}^\dagger + \sum_{j=1}^{\ell-1} P_{j\to \ell}D_jP_{j\to\ell}^\dagger .\] Each transfer \(P_{j\to\ell}\) has \(O(1)\) normalization by stability, so the normalization of the \(j\)-th summand is \(O(h_j^2)\). Since \(\sum_{j=1}^{\ell-1}h_j^2=O(1),\) the LCU normalization of the whole sum is \(O(1)\). The number of terms is \(O(L)\), and the controlled transfer and level-selection overheads are polylogarithmic in the problem size under the usual multilevel access model. ◻
The role of \(W_\ell\) is to choose coordinates for the degrees of freedom added when passing from \(V_{\ell-1}\) to \(V_\ell\). This is a classical issue in finite element multilevel methods. Given nested spaces 4 one often seeks levelwise decompositions \(V_\ell=V_{\ell-1}\oplus W_\ell\) in order to obtain stable multilevel representations, sparse transforms, wavelet-like bases, or well-conditioned preconditioners. There are many ways to build such complements, but the stability in the energy norm is key. A substantial FEM literature therefore modifies raw new-level functions by subtracting suitable coarse components, leading to prewavelets, approximate wavelets, stabilized hierarchical bases, and multilevel Riesz bases [36]–[41].
On the other hand, our use of \(W_\ell\) is local. We do not use it to build a global multilevel preconditioner. It only needs to be constructed for a pair \(V_{\ell-1}\subset V_\ell\) thanks to the decomposition 27 and our QoI-driven task. In our Schur complement approach, the choice of \(W_\ell\) determines the complement coordinate matrices \(J_\ell\) and \(\Sigma_\ell\) in [eq:J-as-ritz-complement,eq:Sigma-JAJ].
The properties needed for the corrected Green’s operator \(D_\ell\) in [Cl-Dl,eq:D-schur] for the QoI estimations in 6 are \[c_0I\preceq \Sigma_\ell\preceq C_0I, \qquad \|J_\ell\|_2=O(h_\ell).\] Classical prewavelet theory gives useful precedents for this type of stability. For example, Floater et al. analyze piecewise linear prewavelets on bounded triangulations and prove that the Schur complement arising in their two-scale matrix is uniformly conditioned [42]. Their Schur complement is not identical to our elliptic energy Schur complement \(\Sigma_\ell\), but it is a close finite element analogue: a carefully chosen two-scale complement can produce a stable Schur complement.
The simplest sufficient condition for our purposes is to choose \(W_\ell\) so that it is energy-orthogonal to the embedded coarse space: \[P_\ell^\dagger A_\ell W_\ell=0.\] Then \[J_\ell=W_\ell, \qquad \Sigma_\ell=W_\ell^\dagger A_\ell W_\ell.\]
For a conforming \(P_1\) finite element space obtained by uniform refinement of a shape-regular polygonal triangulation, a raw complement basis \(W_\ell^{\rm raw}\) can be obtained from the fine nodal basis functions associated with the new vertices created during refinement. Except in special cases such as the one-dimensional midpoint basis with constant coefficient, this raw complement generally has nonzero energy coupling with the coarse space. It can, however, be locally orthogonalized. For this purpose, let \(\widetilde{\psi}_{\omega,j}\) from \(W_\ell^{\rm raw}\) denote a raw fine function associated with a new degree of freedom in a uniformly bounded patch \(\omega\), and let \(p_{\omega,k}\) denote the embedded coarse basis functions whose supports intersect \(\omega\). Here a patch \(\omega\) means a uniformly bounded union of coarse elements surrounding the support of a new fine-level degree of freedom. In two dimensions, for example, the new degree of freedom can be easily created by inserting a nodal point at an edge midpoint, in which case one may take \(\omega\) to be the union of the two coarse triangles sharing that edge, enlarged by one ring of neighboring coarse triangles if needed to include all coarse basis functions whose supports interact through the energy form. The number of fine and coarse degrees of freedom in such a patch is bounded independently of \(h_\ell\) on a shape-regular refinement hierarchy.
Specifically, one can seek modified local functions \[\psi_{\omega,j}^{a} = \widetilde{\psi}_{\omega,j} - \sum_k c_{kj}p_{\omega,k}\] with coefficients chosen so that \[a_\omega(\psi_{\omega,j}^{a},p_{\omega,k}) := \int_\omega \nabla p_{\omega,k}(x)^\dagger a(x)\nabla \psi_{\omega,j}^{a}(x)\,dx =0, \quad \text{for all included coarse function }p_{\omega,k}. \label{eq:local-energy-orthogonalization}\tag{52}\] When \(a(x)\) is constant or elementwise constant on the patch, these integrals are computed from a fixed-size local stiffness matrix. The coefficients \(c_{ij}\) are therefore obtained from a local coarse stiffness problem of size independent of \(h_\ell\). After the coarse component has been removed from this energy orthogonalization, a further local energy diagonalization or Gram–Schmidt step can be applied within each patch to make the local energy Gram matrices uniformly conditioned. Since this diagonalization is performed inside the local subspace already satisfying 52 , it preserves energy orthogonality to the included coarse functions. Thus the resulting local functions remain orthogonal to the coarse functions accounted for in the patch: \(P_\ell^\dagger A_\ell W_\ell^a=0.\) Uniform conditioning of the assembled global matrix \((W_\ell^a)^\dagger A_\ell W_\ell^a\) follows from the finite-overlap stability argument below.
To pass from local patches to the global matrix \(W_\ell^a\), we use the standard finite-overlap argument. Each corrected local function has support in a patch of diameter \(O(h_\ell)\), and each coarse element belongs to only \(O(1)\) such patches. Consequently, the global energy Gram matrix is a sparse sum of uniformly bounded local Gram matrices. If the assembled patch functions form a stable local complement, then \[c_0 I\preceq (W_\ell^a)^\dagger A_\ell W_\ell^a\preceq C_0 I,\] with constants independent of \(h_\ell\).
The required scaling then follows from standard finite element estimates. The corrected local functions have support of diameter \(O(h_\ell)\) and are energy-normalized, so the local Poincaré inequality gives \[\|W_\ell^a\boldsymbol{z}\|_{L^2(\Omega)} \le C h_\ell \|W_\ell^a\boldsymbol{z}\|_a \simeq C h_\ell\|\boldsymbol{z}\|_2.\] In reduced coordinates, this becomes \[\|J_\ell\|_2=\|W_\ell^a\|_2=O(h_\ell).\] Moreover, \[c_0I\preceq \Sigma_\ell\preceq C_0I\] with constants independent of \(h_\ell\). In a nonoverlapping macrocell construction, \(\Sigma_\ell\) can be made block diagonal with constant-size blocks.
This construction is closely related to finite element prewavelets and stabilized hierarchical bases. Vassilevski–Wang modify hierarchical basis functions by subtracting computable approximate \(L^2\)-projections onto coarser finite element spaces, leading to spectrally equivalent multilevel preconditioners [38], [39]. Our \(W_\ell^a\) is the elliptic-energy analogue of these coarse-component subtraction ideas. It is also closely related to the localized orthogonal decomposition viewpoint of [43], where correctors are computed on vertex patches and the localization error decays exponentially in the number of coarse element layers; patch diameters of order \(h\log(h^{-1})\) preserve the coarse convergence rate for heterogeneous elliptic problems.
Thus exact energy orthogonality, \(P_\ell^\dagger A_\ell W_\ell=0\), should be viewed as a clean sufficient condition rather than a necessary one. A natural future direction is to adapt the stable complement constructions from prewavelets, approximate wavelets, multilevel Riesz bases, and localized orthogonal decomposition so that they yield block-encodable Ritz-complement coordinates \(J_\ell\), or alternative corrected-Green’s factorizations, with the same \(O(h_\ell^2)\) response scale but with more flexibility.
We now summarize the cost of the Schur-complement route by applying the MLMC allocation result from 4 and the observable-estimation scaling from 4. The main level-dependent quantities are the effective observable scale and the gate/query cost of one level estimator. By 1, under the energy-stable Ritz-complement coordinate condition and the Ritz-complement oracle assumption, the corrected Green’s operator \(D_\ell=J_\ell\Sigma_\ell^{-1}J_\ell^\dagger\) has a block encoding with normalization \(\alpha_{D_\ell}=O(h_\ell^2).\)
The MLMC identity 27 reduces the QoI estimation problem to contributions from neighboring pairs \((V_{\ell-1},V_\ell)\). The coarsest contribution \(Q_0\) is treated as a fixed coarse problem and is included in the sums below with \(O(1)\) cost. For \(\ell\ge1\), let \(T_\ell\) denote the gate/query cost of one level-\((\ell-1,\ell)\) circuit call for the Schur-complement estimator, including the block encodings of \(J_\ell\), \(A_\ell\), \(\Sigma_\ell^{-1}\), the observable, and, for quadratic QoIs, the lifted coarse response.
In our setting, the stiffness matrix has \(O(1)\) row and column sparsity, and in reduced coordinates has normalization \(\alpha_{A_\ell}=O(h_\ell^{-2}).\) The transfer \(P_\ell\) has \(O(1)\) sparsity and \(O(1)\) norm. In the efficient Ritz-complement constructions of 5.4, \(J_\ell\) is local or localized on patches of radius \(m\), with \[\alpha_{J_\ell}=O(h_\ell), \qquad T_{J_\ell}=\operatorname{poly}(m)\operatorname{polylog}(h_\ell^{-1},\epsilon^{-1}).\] For exact local constructions, \(m=O(1)\); for localized correctors, one typically takes \(m=O(\log(L/\epsilon))\). Since \(\Sigma_\ell=J_\ell^\dagger A_\ell J_\ell\) is \(O(1)\)-normalized and uniformly conditioned, 3 gives a block encoding of \(\Sigma_\ell^{-1}\) with only polylogarithmic precision overhead. Thus the corrected Green’s operator \(D_\ell\) has level cost \(T_\ell=\widetilde{O}(1)\) with respect to powers of \(h_\ell^{-1}\), in the efficient constructions considered here. For quadratic QoIs, the lifted coarse response \(C_\ell=P_\ell B_{\ell-1}P_\ell^\dagger\) also appears in the level-difference identity. We include its implementation cost in \(T_\ell\) estimated in 10.
For a linear QoI, let \(\chi=\chi_{\rm lin}\) be the effective readout order from the level-difference identity 33 . Combining the readout scale \(O(h_\ell^{-\chi})\) with \(\alpha_{D_\ell}=O(h_\ell^2)\), the level observable scale is \[\alpha_\ell^{\rm lin} = O(h_\ell^{2-\chi}). \label{eq:linear-level-scale-schur}\tag{53}\]
For a quadratic QoI, let \(\chi=\chi_{\rm quad}\) be the effective readout order of \(M_\ell\). By 34 , the quadratic level difference contains the terms \(D_\ell^\dagger M_\ell D_\ell, \, D_\ell^\dagger M_\ell C_\ell, \text{and}\, C_\ell^\dagger M_\ell D_\ell .\) The first term has scale \(O(h_\ell^{4-\chi})\), while the cross terms have scale \(O(h_\ell^{2-\chi})\). Hence the effective quadratic level scale is \[\alpha_\ell^{\rm quad} = O(h_\ell^{2-\chi}). \label{eq:quadratic-level-scale-schur}\tag{54}\] The coarse response affects the circuit cost \(T_\ell\), but it does not change the fine-level correction scale.
Therefore, for both linear and quadratic QoIs, one has \(\alpha_\ell=O(h_\ell^{2-\chi}),\) with \(0\le \chi\le 2.\) Using the amplitude-estimation allocation in 5, with \(c_\ell=T_\ell\alpha_\ell\), gives \[\mathcal{C}_{\rm Schur}^{\rm AE} = \widetilde{O}\!\left( \frac{1}{\epsilon} \left[ \sum_{\ell=0}^{L} \left(T_\ell h_\ell^{2-\chi}\right)^{1/2} \right]^2 \right). \label{eq:schur-AE-complexity}\tag{55}\] Since \(0\le \chi\le2\) and \(h_\ell=h_{\ell-1}/2\), the geometric sum is bounded independently of \(h=h_L\), with only logarithmic factors at the endpoint \(\chi=2\). Hence \[\mathcal{C}_{\rm Schur}^{\rm AE} = \widetilde{O}(\epsilon^{-1}), \qquad 0\le\chi\le2. \label{eq:schur-AE-complexity-simplified}\tag{56}\]
For direct sampling, one level-\(\ell\) sample has variance bounded by \(O(\alpha_\ell^2)\). Applying the sampling allocation in 5 gives \[\mathcal{C}_{\rm Schur}^{\rm samp} = \widetilde{O}\!\left( \frac{1}{\epsilon^2} \left[ \sum_{\ell=0}^{L} T_\ell^{1/2}h_\ell^{2-\chi} \right]^2 \right). \label{eq:schur-sampling-complexity}\tag{57}\] Again, the sum is bounded up to logarithmic factors at \(\chi=2\). Therefore \[\mathcal{C}_{\rm Schur}^{\rm samp} = \widetilde{O}(\epsilon^{-2}). \label{eq:schur-sampling-complexity-simplified}\tag{58}\] This gives the Schur-complement multilevel complexity summarized in 1.
Classically, MLMC couples solvers at two consecutive levels [32]. It is therefore natural to ask whether a quantum algorithm can also prepare, at level \(\ell\), a state containing both the coarse response and the fine–coarse correction. This section speculates such a possible route. The algebra is simple and the readout scaling is favorable. The open point is whether the resulting scaled two-level system can be solved by a QLSA while preserving the physical scaled state needed for the QoI readout.
We start from the Galerkin-compatible coarse load and response \(\boldsymbol{b}_{\ell-1}:=P_\ell^\dagger\boldsymbol{b}_\ell, \qquad \boldsymbol{x}_{\ell-1}:=B_{\ell-1}\boldsymbol{b}_{\ell-1},\) and define the interlevel correction \[\boldsymbol{e}_\ell := \boldsymbol{x}_\ell-P_\ell\boldsymbol{x}_{\ell-1}. \label{eq:two-level-error}\tag{59}\] Then \[A_\ell\boldsymbol{e}_\ell+A_\ell P_\ell\boldsymbol{x}_{\ell-1}=\boldsymbol{b}_\ell, \qquad A_{\ell-1}\boldsymbol{x}_{\ell-1}=\boldsymbol{b}_{\ell-1}. \label{eq:two-level-correction-equations}\tag{60}\]
For standard \(P_1/Q_1\) elements, the standard finite element estimate gives \[\left\lVert \boldsymbol{x}_\ell-P_\ell\boldsymbol{x}_{\ell-1}\right\rVert_2=O(h_\ell^2)\] under the usual regularity assumptions. Motivated by this scale, define the extended unknown \[\boldsymbol{z}_\ell := \begin{bmatrix} h_\ell^{-2}\boldsymbol{e}_\ell\\[1mm] \boldsymbol{x}_{\ell-1} \end{bmatrix}. \label{eq:scaled-extended-unknown}\tag{61}\] Then \(\boldsymbol{z}_\ell\) solves the Hermitian positive definite system \[S_\ell \boldsymbol{z}_\ell=\boldsymbol{r}_\ell, \label{eq:scaled-spd-extended-system}\tag{62}\] with \[S_\ell := \begin{bmatrix} h_\ell^{4}A_\ell & h_\ell^{2}A_\ell P_\ell\\[1mm] h_\ell^{2}P_\ell^\dagger A_\ell & 2A_{\ell-1} \end{bmatrix}, \qquad \boldsymbol{r}_\ell := \begin{bmatrix} h_\ell^{2}\boldsymbol{b}_\ell\\[1mm] 2\boldsymbol{b}_{\ell-1} \end{bmatrix}. \label{eq:S-r-scaled}\tag{63}\]
The advantage of the scaling is visible in the readout. For a linear QoI, \[\Delta_\ell^{\rm lin} = \boldsymbol{g}_\ell^\dagger\boldsymbol{e}_\ell = (\boldsymbol{g}_\ell^{\rm ext})^\dagger\boldsymbol{z}_\ell, \qquad \boldsymbol{g}_\ell^{\rm ext} := \begin{bmatrix} h_\ell^2\boldsymbol{g}_\ell\\[1mm] 0 \end{bmatrix}. \label{eq:linear-ext-readout}\tag{64}\] For a quadratic QoI, \[\Delta_\ell^{\rm quad} = \boldsymbol{z}_\ell^\dagger M_\ell^{\rm ext}\boldsymbol{z}_\ell, \quad M_\ell^{\rm ext} := \begin{bmatrix} h_\ell^4 M_\ell & h_\ell^2M_\ell P_\ell\\[1mm] h_\ell^2P_\ell^\dagger M_\ell & 0 \end{bmatrix}. \label{eq:M-ext-scaled}\tag{65}\] The \(h_\ell^4M_\ell\) block gives \(\boldsymbol{e}_\ell^\dagger M_\ell\boldsymbol{e}_\ell\), while the off-diagonal blocks give the two cross terms.
If the effective quadratic observable scale satisfies \(\alpha_{M_\ell}^{\rm eff}=O(h_\ell^{-\chi}), \qquad \chi\le 2,\) and \(\left\lVert P_\ell\right\rVert_2=O(1)\), then the extended observable has scale \[\alpha_{M_\ell^{\rm ext}}^{\rm eff} = O(h_\ell^{2-\chi}). \label{eq:M-ext-scale}\tag{66}\] Thus, if a QLSA could prepare the normalized extended state \(\left|z_\ell\right\rangle = \frac{\boldsymbol{z}_\ell}{Z_\ell}, \qquad Z_\ell:=\left\lVert \boldsymbol{z}_\ell\right\rVert_2,\) efficiently, then the physical quadratic level difference would be recovered as \(\Delta_\ell^{\rm quad} = Z_\ell^2 \left\langle z_\ell\right|M_\ell^{\rm ext}\left|z_\ell\right\rangle.\) The norm \(Z_\ell\) must be estimated separately. Since \(Z_\ell^2 = h_\ell^{-4}\left\lVert \boldsymbol{e}_\ell\right\rVert_2^2+\left\lVert \boldsymbol{x}_{\ell-1}\right\rVert_2^2 = O(1),\) the norm estimation adds precision dependence but no additional power of \(h_\ell\). By 4, the scale 66 would then give the same level readout improvement as in the Schur-complement route.
The difficulty, however, is the linear-system step. A naive block encoding of \(S_\ell\) is not well scaled: the blocks carry different powers of \(h_\ell\), and the scaled correction component may introduce an unfavorable condition number if treated as an ordinary linear system. A preconditioner for 62 must therefore do more than condition the solve; it must also allow preparation of the physical scaled state \(\left|z_\ell\right\rangle\), or an equivalent readout representation in which the small factor \(h_\ell^2\) in 64 –65 remains visible to the measurement procedure. This is not automatic for a BPX-frame QLSA, we therefore regard the scaled two-level QLSA formulation as a promising open route rather than as a completed algorithm in this paper.
This paper studied the end-to-end complexity of quantum elliptic PDE solvers. The central point of this paper is that the simulation stage and the inference stage should be co-designed. By coupling neighboring finite element levels before measurement, the quantum circuit can see the scale of the interlevel correction rather than the scale of the full fine-grid observable.
We developed a Schur-complement approach in which the corrected Green’s operator avoids solution-state preparation and expresses the level increment through a two-level Ritz-complement and Schur-complement. Under the Ritz-complement oracle assumptions, this exposes the \(O(h_\ell^2)\) scale of the finite element correction inside the block-encoded observable. As a result, for observables up to the energy level, the polynomial \(h\)-dependent readout overhead is removed up to logarithmic factors. With amplitude estimation, the remaining dependence on the statistical precision is \(\widetilde{O}(\epsilon^{-1})\), i.e., Heisenberg scaling for the inference stage; with direct sampling, the dependence is the usual \(\widetilde{O}(\epsilon^{-2})\). This is the sense in which the approach moves from a linear-system subroutine toward a more complete prepare–simulate–infer (PSI) pipeline.
From a circuit perspective, the multilevel estimator reduces how often expensive fine-level circuits must be executed. Amplitude estimation repeats controlled state-preparation and block-encoding circuits a number of times proportional to the scale of the observable being estimated; standard sampling repeats the measurement circuit according to the corresponding variance. By inserting the fine–coarse cancellation before measurement, the level-\(\ell\) circuit measures a small correction rather than a full fine-grid observable. This reduces the total number of fine-level oracle calls and measurements. When the Ritz-complement oracle is local or two-level, it may also reduce the depth and structural complexity of the circuit used at that level. Similar to [44] our current observable driven approach also reduces the overall complexity.
The work also points to a broader connection with finite element multilevel theory. The Ritz-complement construction used here is only one way of exposing a stable interlevel correction. The finite element literature contains many powerful constructions of preconditioners, hierarchical complements, prewavelets, approximate wavelets,and localized orthogonal decompositions. These methods were developed for classical stability and preconditioning, but they may also provide new block-encoding primitives for quantum observable estimation. Understanding which of these classical constructions lead to small-normalization quantum oracles is an important direction for future work.
The same readout principle also appears naturally in spectral discretizations. Recent quantum spectral-filtering approaches exploit Fourier-space filters to solve structured second-order PDEs by block encoding the inverse response in the diagonal spectral basis [45]. In such settings, dyadic frequency bands provide an especially transparent version of our construction: energy orthogonality is automatic, and the inverse restricted to a band with frequency scale \(h_\ell^{-1}\) has size \(O(h_\ell^2)\). Thus the multilevel readout idea can be viewed as inserting dyadic Green’s filters directly into the observable estimator, rather than using spectral filtering only to prepare a normalized solution state.
Several other extensions are natural. Nonsymmetric or indefinite PDEs may require saddle-point analogues of the present Schur-complement formulation. Time-dependent problems, including diffusion and wave equations are another important direction. In all of these settings, the main lesson is the same: quantum PDE algorithms should be analyzed end-to-end, with the quantity of interest and the measurement procedure built into the algorithm from the get go.
This work was supported by NSF Grant DMS-2411120.