Variational Quantum Linear Solver via Block Encoding for the Poisson Equation

Abstract: We present a variational quantum linear solver (VQLS) for the Poisson equation built on an exact block encoding of the discrete Laplacian, and demonstrate its performance on physically motivated benchmarks. Unlike LCU-based VQLS where the number of distinct circuits required per cost-function evaluation is $\mathcal{O}(L<sup>2)$, where $L$ is the number of terms in the LCU decomposition of the discrete Laplacian operator, this approach requires only a single circuit for cost evaluation. We further empirically demonstrate that the choice of classical optimizer materially affects where the variational optimization ceases to make progress. The solver is benchmarked on three problems: a Poisson equation with sinusoidal forcing and a steady-state heat conduction problem with a localized Gaussian source, both with Dirichlet boundaries, and the pressure-Poisson equation of a two-dimensional lid-driven cavity flow, in which the solver is invoked once per time step under Neumann boundary conditions.

• The paper develops an exact four-ancilla block encoding of the finite-difference Laplacian that evaluates the VQLS cost with one circuit instead of the O(L²) circuits required by Pauli-LCU methods.
• The paper finds post-selection success probabilities scale as N^-0.57 for Dirichlet and N^-0.94 for Neumann conditions, while block-encoded circuits become less resource-intensive than LCU circuits beyond N=512.
• The paper shows CMA-ES outperforms COBYLA on larger systems and supports short-horizon cavity-flow simulation, but accuracy degradation, barren plateaus, noise, and full-solution readout remain open challenges.

Overview

This paper develops a variational quantum linear solver (VQLS) for the discretized Poisson equation in which the system operator is supplied by an exact block encoding of the finite-difference Laplacian, rather than by a linear combination of unitaries (LCU). The central methodological claim is that the global cost function can be evaluated from the measurement statistics of a single quantum circuit, whereas the standard LCU-based VQLS formulation requires $\mathcal{O}(L^2)$ Hadamard-test circuits per cost evaluation, where $L$ is the number of terms in the Pauli decomposition. The authors characterize the approach on three benchmarks—a Poisson problem with sinusoidal forcing, a steady-state heat conduction problem with a localized Gaussian source (both with Dirichlet boundaries), and the pressure-Poisson equation inside a projection-method lid-driven cavity simulation under Neumann conditions—and they empirically connect optimizer choice to the onset of barren-plateau stagnation.

Block encoding of the discrete Laplacian

The 2D Laplacian inherits a Kronecker-sum structure $$A = I_{N_y} \otimes T_{N_x} + T_{N_y} \otimes I_{N_x}$$ from the five-point stencil, and this separability is exploited directly. Because the spectral norm grows as $\mathcal{O}(1/h_{\min}^2)$, the operator is rescaled as $\widetilde{A} = (\Lambda/4)A$ with $\Lambda = (1/\Delta x^2 + 1/\Delta y^2)^{-1}$, guaranteeing $\|\widetilde{A}\|_2 \leq 1$. A one-qubit selector register prepared in the superposition $$\sqrt{\omega_1}|0\rangle + \sqrt{\omega_2}|1\rangle$$ routes between one-dimensional block-encoding sub-circuits for each axis; Dirichlet closure suppresses wrap-around couplings from cyclic shifts via a deletion ancilla, while Neumann closure additionally deletes one diagonal branch at each boundary node, reducing the boundary diagonal from $-2$ to $-1$.

The result is an exact $L$0-block encoding using four ancilla qubits, verified numerically to double-precision machine epsilon across all sizes from $L$1 to $L$2 for both boundary conditions. The construction follows a companion work (Boutot et al., 12 Mar 2026) and carries gate complexity $L$3, reducing to $L$4 in two dimensions. The Neumann operator is singular—its null space is spanned by the constant vector—and the paper fixes the gauge by projecting onto the zero-mean subspace, which leaves the pressure gradient and corrected velocity unchanged.

Cost evaluation from a single circuit

The normalized global cost function $L$5 measures the infidelity between $L$6 and $L$7. Using Bayes' theorem applied to post-selection statistics, the cost reduces to

$L$8

so both numerator and denominator come from one circuit's measurement outcomes. Solution magnitude is recovered coherently from the same post-selection probability, at shot cost $L$9 independent of $$A = I_{N_y} \otimes T_{N_x} + T_{N_y} \otimes I_{N_x}$$0.

The critical question is whether the single-circuit reduction survives the post-selection overhead. Empirically it does: fitting power laws to the success probability at optimized parameters gives $$A = I_{N_y} \otimes T_{N_x} + T_{N_y} \otimes I_{N_x}$$1 (Dirichlet) and $$A = I_{N_y} \otimes T_{N_x} + T_{N_y} \otimes I_{N_x}$$2 (Neumann), far slower than the $$A = I_{N_y} \otimes T_{N_x} + T_{N_y} \otimes I_{N_x}$$3 decay that the shrinking spectral gap might suggest. Against the Pauli-LCU formulation, which yields exactly $$A = I_{N_y} \otimes T_{N_x} + T_{N_y} \otimes I_{N_x}$$4 terms and hence $$A = I_{N_y} \otimes T_{N_x} + T_{N_y} \otimes I_{N_x}$$5 circuits per cost evaluation ($$A = I_{N_y} \otimes T_{N_x} + T_{N_y} \otimes I_{N_x}$$6, growing from 12 circuits at $$A = I_{N_y} \otimes T_{N_x} + T_{N_y} \otimes I_{N_x}$$7 to 9120 at $$A = I_{N_y} \otimes T_{N_x} + T_{N_y} \otimes I_{N_x}$$8), the block-encoded formulation needs exactly one circuit at every size. Per-circuit resources are more nuanced: LCU circuits are individually cheaper at small sizes, but the ordering reverses beyond $$A = I_{N_y} \otimes T_{N_x} + T_{N_y} \otimes I_{N_x}$$9; at $\mathcal{O}(1/h_{\min}^2)$0 the block-encoded circuit requires 16,034 two-qubit gates against 28,980 for the heaviest LCU circuit, at depths of 28,741 versus 45,744. The only axis on which the LCU wins is width, by three ancilla qubits. The authors note that the Pauli decomposition is close to the most favorable LCU available in per-circuit depth, so any compact non-Pauli decomposition would be expected to lose the crossover at an even smaller size.

Solution accuracy

Accuracy separates sharply between the two Dirichlet benchmarks, and the reason is structural. The sinusoidal source is an exact eigenvector of the discrete Dirichlet Laplacian, so its solution state coincides with $\mathcal{O}(1/h_{\min}^2)$1 and the ansatz need only prepare a Schmidt-rank-one state across the $\mathcal{O}(1/h_{\min}^2)$2–$\mathcal{O}(1/h_{\min}^2)$3 cut. The solver recovers it to arithmetic precision at every size—the median seed attains infidelity $\mathcal{O}(1/h_{\min}^2)$4 at final cost $\mathcal{O}(1/h_{\min}^2)$5 on the $\mathcal{O}(1/h_{\min}^2)$6 grid—but this verifies the pipeline end to end without probing scaling behavior.

Grid	$\mathcal{O}(1/h_{\min}^2)$7	Sine $\mathcal{O}(1/h_{\min}^2)$8	Heat $\mathcal{O}(1/h_{\min}^2)$9
$\widetilde{A} = (\Lambda/4)A$0	16	$\widetilde{A} = (\Lambda/4)A$1	$\widetilde{A} = (\Lambda/4)A$2
$\widetilde{A} = (\Lambda/4)A$3	64	$\widetilde{A} = (\Lambda/4)A$4	$\widetilde{A} = (\Lambda/4)A$5
$\widetilde{A} = (\Lambda/4)A$6	256	$\widetilde{A} = (\Lambda/4)A$7	$\widetilde{A} = (\Lambda/4)A$8
$\widetilde{A} = (\Lambda/4)A$9	1024	$\Lambda = (1/\Delta x^2 + 1/\Delta y^2)^{-1}$0	$\Lambda = (1/\Delta x^2 + 1/\Delta y^2)^{-1}$1

For the Gaussian heat source, whose solution has overlap $\Lambda = (1/\Delta x^2 + 1/\Delta y^2)^{-1}$2 with $\Lambda = (1/\Delta x^2 + 1/\Delta y^2)^{-1}$3 at $\Lambda = (1/\Delta x^2 + 1/\Delta y^2)^{-1}$4, accuracy holds at machine precision through $\Lambda = (1/\Delta x^2 + 1/\Delta y^2)^{-1}$5 and then degrades by roughly an order of magnitude per grid-size increase. The paper attributes this degradation to optimization rather than to the encoding.

The third benchmark embeds the solver in a projection-method Navier–Stokes loop for the lid-driven cavity at $\Lambda = (1/\Delta x^2 + 1/\Delta y^2)^{-1}$6 on a $\Lambda = (1/\Delta x^2 + 1/\Delta y^2)^{-1}$7 grid, invoking the Neumann VQLS once per time step with warm starts from the previous step's converged parameters. Over ten steps to $\Lambda = (1/\Delta x^2 + 1/\Delta y^2)^{-1}$8, both solvers place the primary vortex near $\Lambda = (1/\Delta x^2 + 1/\Delta y^2)^{-1}$9 and track each other in centreline velocity profiles, though the quantum pressure solve systematically damps velocity magnitude while preserving field structure. Since pressure errors propagate into subsequent initial conditions, the fact that trajectories remain together over ten steps indicates the attained accuracy suffices for time advancement—at least over this short horizon while the flow is still developing.

Barren plateaus and optimizer choice

Direct measurement confirms a barren plateau in the sense of exponentially decaying gradient variance under uniform sampling: fits give $\|\widetilde{A}\|_2 \leq 1$0 with $\|\widetilde{A}\|_2 \leq 1$1, $\|\widetilde{A}\|_2 \leq 1$2, $\|\widetilde{A}\|_2 \leq 1$3, $\|\widetilde{A}\|_2 \leq 1$4 for ansatz depths $\|\widetilde{A}\|_2 \leq 1$5 to $\|\widetilde{A}\|_2 \leq 1$6—statistically indistinguishable across depths, consistent with the plateau being a property of the global cost function rather than of ansatz expressivity [2209.xxxxx; cerezo2021cost]. Sampled instead along the CMA-ES trajectory, gradient variances sit one to three orders of magnitude above the uniform baseline at every size, with depth-dependent apparent decay rates. The authors are careful here: elevated variance along the trajectory is partly a consequence of successful optimization as much as a cause, so no causal mechanism is claimed—only that the regions CMA-ES occupies are atypical of the landscape, in the manner anticipated by the narrow-gorge picture [2107.xxxxx].

At equal budget ($\|\widetilde{A}\|_2 \leq 1$7 cost evaluations), COBYLA descends rapidly then flattens entirely at $\|\widetilde{A}\|_2 \leq 1$8, while CMA-ES continues reducing cost throughout. CMA-ES also offers embarrassingly parallel population evaluation, unlike sequential COBYLA. Importantly, the paper states plainly that CMA-ES does not remove the barren plateau—it postpones the system size at which stagnation sets in.

Limitations and open questions

Several caveats bear directly on the results. All experiments run on exact statevector simulation (Qiskit Aer); shot requirements implied by the success-probability scaling are quantified analytically but not measured, and noise resilience is entirely uncharacterized—both prerequisites for any hardware demonstration. The sinusoidal benchmark, being an exact eigenvector case, does not test scaling behavior, and the Gaussian-source degradation beyond $\|\widetilde{A}\|_2 \leq 1$9 remains attributed to optimization only by inference. The empirical advantage of CMA-ES lacks theoretical grounding, and whether covariance-adapting search continues to locate atypical high-gradient regions as $$\sqrt{\omega_1}|0\rangle + \sqrt{\omega_2}|1\rangle$$0 grows is the open question the work raises most directly. The hardware-efficient ansatz treats the Laplacian as a generic operator; ansätze exploiting Kronecker-sum structure are unexplored. Finally, the classical readout bottleneck—$$\sqrt{\omega_1}|0\rangle + \sqrt{\omega_2}|1\rangle$$1 measurements to extract the full solution vector—remains unresolved, and identifying engineering-relevant low-dimensional functionals admitting efficient observable estimation is left open.

Conclusion

This paper establishes block-encoded VQLS as a concrete alternative to LCU-based formulations for elliptic PDE systems, with a rigorous resource argument: one circuit per cost evaluation, slower-than-$$\sqrt{\omega_1}|0\rangle + \sqrt{\omega_2}|1\rangle$$2 post-selection decay, and per-circuit cost advantages beyond $$\sqrt{\omega_1}|0\rangle + \sqrt{\omega_2}|1\rangle$$3, at the price of three ancilla qubits. Accuracy results are honest about their limits—arithmetic precision on eigenvector-aligned problems, order-of-magnitude-per-grid degradation elsewhere, attributable to the barren plateau that the measurements confirm. The demonstration that the solver sustains a time-advancing cavity flow simulation, even briefly and at modest resolution, is the most application-relevant result. The encoding circuits are directly reusable within QSVT-based fault-tolerant solvers, giving the construction value beyond the NISQ variational setting.