Variational Quantum Linear Solver (VQLS)

• VQLS is a hybrid quantum-classical framework that probabilistically prepares a quantum state representing the solution A⁻¹|b⟩ to linear systems.
• It employs variational state preparation with hardware-efficient ansätze and diverse cost functions to optimize convergence and mitigate noise.
• Advanced decomposition techniques and distributed measurement strategies enable VQLS to efficiently address structured PDEs, quantum transport, and other applications.

Variational Quantum Linear Solver (VQLS) is a hybrid quantum-classical algorithm for solving linear systems on near-term quantum computers. In its standard form, the task is to prepare a parameterized quantum state $|x(\theta)\rangle = V(\theta)|0^n\rangle$ such that $A|x(\theta)\rangle \propto |b\rangle$, where the coefficient matrix is represented as a linear combination of unitaries $A=\sum_{l=1}^L c_l A_l$ and the right-hand side is prepared as $|b\rangle = U|0^n\rangle$ (Bravo-Prieto et al., 2019). The method replaces deep phase-estimation-style subroutines with variational state preparation, repeated quantum estimation of a cost function, and classical parameter optimization. Subsequent work has applied and extended this framework to tridiagonal Poisson systems, heat conduction, potential and Stokes flows, nonlinear dynamics, PDE-constrained optimization, modulo-2 systems, quantum machine learning, quantum transport, and distributed solvers (Balducci et al., 2024, 2207.14630, Liu et al., 2023, Hafshejani et al., 2024, Aboumrad et al., 2023, Yi et al., 2023, Yang et al., 6 Sep 2025, Lu et al., 15 Apr 2026).

1. Formal definition and problem setting

The core problem addressed by VQLS is the quantum linear system problem. Given a full-rank matrix $A$, condition number $\kappa$, and an efficiently preparable state $|b\rangle$, the goal is to produce a quantum state proportional to $A^{-1}|b\rangle$. In the canonical formulation, $A$ is written as a linear combination of unitaries, $A=\sum_{l=1}^L c_l A_l$, and a variational circuit $A|x(\theta)\rangle \propto |b\rangle$0 prepares $A|x(\theta)\rangle \propto |b\rangle$1. Defining $A|x(\theta)\rangle \propto |b\rangle$2, training seeks a parameter set for which $A|x(\theta)\rangle \propto |b\rangle$3 is proportional to $A|x(\theta)\rangle \propto |b\rangle$4 (Bravo-Prieto et al., 2019).

This formulation naturally fits amplitude encoding. For a system of size $A|x(\theta)\rangle \propto |b\rangle$5, both the unknown vector and the right-hand side are represented on $A|x(\theta)\rangle \propto |b\rangle$6 qubits. A perfect solution satisfies $A|x(\theta)\rangle \propto |b\rangle$7, with $A|x(\theta)\rangle \propto |b\rangle$8, so the quantum state captures the classical solution up to overall scale (Turati et al., 2024). In PDE-oriented applications, this encoding is used after discretization. The 1D Poisson problem yields a tridiagonal matrix with $A|x(\theta)\rangle \propto |b\rangle$9 and $A=\sum_{l=1}^L c_l A_l$0 on $A=\sum_{l=1}^L c_l A_l$1 grid points (Balducci et al., 2024), while the heat conduction equation leads to 1D tridiagonal and 2D block-Toeplitz systems under finite-difference discretization (2207.14630).

The original framework assumes Hermitian structure or an embedding into a Hermitian block, but later work broadens this scope. For Carleman-linearized nonlinear dynamics, one may use regularized normal equations $A=\sum_{l=1}^L c_l A_l$2 or an augmented-system dilation $A=\sum_{l=1}^L c_l A_l$3 (Liu et al., 14 May 2026). For quantum transport, where the NEGF coefficient matrix is complex and non-Hermitian, alternative cost constructions are introduced directly for $A=\sum_{l=1}^L c_l A_l$4 (Yang et al., 6 Sep 2025). This establishes VQLS not as a single fixed routine but as a family of variational linear-solver constructions sharing a common objective.

2. Cost functions and solution certification

The standard VQLS cost is the normalized global overlap cost

$A=\sum_{l=1}^L c_l A_l$5

which vanishes if and only if $A=\sum_{l=1}^L c_l A_l$6. A local alternative replaces the global projector with qubit-wise projectors and is designed to improve trainability for larger systems (Bravo-Prieto et al., 2019). An equivalent expression used in later empirical studies is

$A=\sum_{l=1}^L c_l A_l$7

with $A=\sum_{l=1}^L c_l A_l$8 and $A=\sum_{l=1}^L c_l A_l$9 (Turati et al., 2024).

A distinctive feature of the original VQLS analysis is its operational termination condition. The global and local costs satisfy

$|b\rangle = U|0^n\rangle$0

where $|b\rangle = U|0^n\rangle$1 is the trace-distance error of the prepared solution state and $|b\rangle = U|0^n\rangle$2 is the condition number of $|b\rangle = U|0^n\rangle$3 (Bravo-Prieto et al., 2019). This gives VQLS a stopping rule tied to solution precision rather than only optimizer stagnation.

Later work introduces alternative objectives for settings not well served by the original overlap form. For complex, non-symmetric transport matrices, the normalized-residual cost is

$|b\rangle = U|0^n\rangle$4

with an analytically optimized real scale $|b\rangle = U|0^n\rangle$5, and a hybrid cost

$|b\rangle = U|0^n\rangle$6

is used to improve convergence (Yang et al., 6 Sep 2025). For compressible-flow and Carleman-linearized applications, the squared-residual form $|b\rangle = U|0^n\rangle$7 and the Hamiltonian form $|b\rangle = U|0^n\rangle$8 are both employed (Yao et al., 28 Aug 2025, Gnanasekaran et al., 2024). A plausible implication is that the choice of cost function is now treated as a problem-dependent design variable rather than a fixed component of the algorithm.

3. Ansatz families and classical optimization

Most VQLS implementations use hardware-efficient ansätze. A common architecture alternates single-qubit $|b\rangle = U|0^n\rangle$9 rotations with entangling layers built from CNOT or CZ gates (Pellow-Jarman et al., 2021, Turati et al., 2024). In the 1D Poisson pilot study, the ansatz is deliberately minimal: $A$0 with no entangling gates in the ansatz itself, making the variational layer extremely shallow (Balducci et al., 2024). Other applications adopt more expressive templates, including Strongly Entangling Layers with $A$1 rotations for the Lorenz system (Hafshejani et al., 2024) and real-valued multi-block ansätze for PDE-constrained optimization (Gnanasekaran et al., 2024, Surana et al., 2024).

Ansatz adaptation has become a significant theme. The dynamic-ansatz VQLS begins with one layer and appends new layers when progress stalls, using the convergence threshold $A$2 and a switching hyperparameter $A$3 (Patil et al., 2021). Its Total Resource Cost, defined as cumulative layer count across iterations, is used to quantify quantum-resource usage. Benchmarks show that the dynamic ansatz can outperform a fixed-depth ansatz especially for ill-conditioned matrices, larger qubit counts, and in the presence of noise (Patil et al., 2021). In another direction, a multi-ansatz tree represents the solution as

$A$4

so that expressivity is increased by superposing several shallow branches rather than deepening a single circuit; increasing the number of branches and using domain decomposition improved convergence and stability in compressible-flow simulations (Yao et al., 28 Aug 2025).

Classical optimization is likewise problem-dependent. COBYLA appears in several application studies (Balducci et al., 2024, Yi et al., 2023), while ADAM, Adagrad, gradient descent, SPSA, BFGS, AMSGrad, Powell, Nelder-Mead, and conjugate gradient have all been benchmarked (2207.14630, Gnanasekaran et al., 2024, Pellow-Jarman et al., 2021). A direct optimizer comparison for VQLS reports that SPSA appears to be the best performing method under realistic noise, with AMSGrad and BFGS also relatively robust, whereas COBYLA, Nelder-Mead and Conjugate-Gradient methods appear to be the most heavily affected by noise (Pellow-Jarman et al., 2021). More recently, PVLS introduces graph-neural-network-based parameter prediction and reports up to a $A$5 speedup in optimization with $A$6 time saving while maintaining comparable solution accuracy (Yang, 4 Dec 2025).

4. Cost evaluation, LCU structure, and measurement overhead

The main practical bottleneck in VQLS is cost evaluation. In a Pauli-LCU decomposition $A$7, expectation values such as $A$8, $A$9, and their local-cost analogues are estimated through Hadamard-test or overlap-test circuits. For a naive Pauli-LCU decomposition one needs $\kappa$0 distinct controlled-unitary circuits per cost evaluation, and for a generic Hermitian matrix $\kappa$1 can be as large as $\kappa$2 in the Pauli basis; algebraic simplifications reduce this to roughly $\kappa$3 Hadamard tests, which is still $\kappa$4 (Turati et al., 2024). Circuit depth also matters: after transpilation, denominator circuits scale approximately linearly with $\kappa$5, while numerator circuits are $\kappa$6–$\kappa$7 deeper than denominator circuits in the reported experiments (Turati et al., 2024).

Because of this bottleneck, decomposition strategy has become central to VQLS research. For the 1D Poisson tridiagonal matrix, a pure Pauli decomposition uses exactly $\kappa$8 terms, whereas a mixed SWAP/center-switch decomposition uses $\kappa$9 terms. In the $|b\rangle$0 example, the maximum circuit-depth ratio is reported as $|b\rangle$1 for Pauli versus multi-qubit decomposition, so the gain in term count is offset by deeper Hadamard tests (Balducci et al., 2024). For structured sparse matrices arising from PDE discretizations, an alternate non-unitary $|b\rangle$2-basis can reduce the number of tensor-product terms from $|b\rangle$3 in the worst-case Pauli basis to $|b\rangle$4, using unitary completion to make the corresponding circuits measurable (Gnanasekaran et al., 2024). In a heat-equation example with $|b\rangle$5, the reported counts are approximately $|b\rangle$6 Pauli strings versus $|b\rangle$7 $|b\rangle$8-terms (Gnanasekaran et al., 2024).

Recent work pushes this compression further. A fast Walsh-Hadamard transform based Pauli decomposition with $|b\rangle$9 coefficient thresholding reduces the number of retained terms from $A^{-1}|b\rangle$0 to $A^{-1}|b\rangle$1 for structured sparse matrices with $A^{-1}|b\rangle$2 qubits, and for a $A^{-1}|b\rangle$3-qubit tridiagonal Toeplitz system the per-iteration circuit count drops from $A^{-1}|b\rangle$4 million to $A^{-1}|b\rangle$5, a $A^{-1}|b\rangle$6 reduction, while preserving over $A^{-1}|b\rangle$7 solution fidelity (Lu et al., 15 Apr 2026). For Carleman-linearized nonlinear dynamics, symmetry-grouped Hadamard tests reduce the number of distinct measurements by approximately $A^{-1}|b\rangle$8 by exploiting $A^{-1}|b\rangle$9 in the symmetric cost sum (Liu et al., 14 May 2026). These developments indicate that VQLS performance depends at least as much on matrix decomposition and measurement design as on the variational ansatz itself.

5. Applications and domain-specific variants

VQLS has been used extensively for discretized differential equations. In heat conduction, statevector simulations on Pennylane’s backend reached final cost $A$0 and $A$1 fidelity for a $A$2 1D system on $A$3 qubits, while 2D cases on $A$4 and $A$5 interior grids achieved $A$6 in approximately $A$7 evaluations (2207.14630). For potential and Stokes flows, VQLS was combined with finite differences and generalized differential quadrature, with numerical solutions reported on up to $A$8 qubits and heuristic scaling $A$9 (Liu et al., 2023). In the 1D Poisson study, the authors report the first simulated and real-hardware results obtained by solving tridiagonal linear systems with VQLS; on ibmq_athens, the $A=\sum_{l=1}^L c_l A_l$0 case yielded fidelities per run all $A=\sum_{l=1}^L c_l A_l$1, and the $A=\sum_{l=1}^L c_l A_l$2 case produced final fidelities approximately $A=\sum_{l=1}^L c_l A_l$3 and $A=\sum_{l=1}^L c_l A_l$4 in two repeats despite hardware noise (Balducci et al., 2024).

Nonlinear dynamics has motivated several further adaptations. A hybrid quantum solver for the Lorenz system rewrites each forward-Euler step as an $A=\sum_{l=1}^L c_l A_l$5 linear system and solves it on $A=\sum_{l=1}^L c_l A_l$6 qubits using VQLS; the relative trajectory error remains $A=\sum_{l=1}^L c_l A_l$7 over $A=\sum_{l=1}^L c_l A_l$8 steps, and at least $A=\sum_{l=1}^L c_l A_l$9 Strongly Entangling Layers are needed for good convergence (Hafshejani et al., 2024). Variational quantum frameworks for linear and nonlinear PDE-constrained optimization use VQLS as the inner linear-system solver inside a bi-level loop, with Bayesian optimization or other black-box optimizers in the outer loop (Surana et al., 2024, Gnanasekaran et al., 2024). For compressible flows, a multi-ansatz VQLS with domain decomposition solves linear systems from implicit one-dimensional Navier-Stokes updates and captures shock, rarefaction, and contact discontinuities in shock-tube tests (Yao et al., 28 Aug 2025).

The framework has also been specialized beyond conventional real or Hermitian linear algebra. Mod2VQLS solves $A|x(\theta)\rangle \propto |b\rangle$00 by using a Clifford matrix-multiply circuit $A|x(\theta)\rangle \propto |b\rangle$01 and a cost $A|x(\theta)\rangle \propto |b\rangle$02; low-dimensional experiments indicate performance on a par with the block Wiedemann algorithm in the small-$A|x(\theta)\rangle \propto |b\rangle$03 regime (Aboumrad et al., 2023). A VQLS-enhanced QSVM embeds the least-squares SVM linear system into a $A|x(\theta)\rangle \propto |b\rangle$04-qubit VQLS routine; SVD plus Pauli decomposition reduced the number of Pauli terms from $A|x(\theta)\rangle \propto |b\rangle$05 to $A|x(\theta)\rangle \propto |b\rangle$06 and the Hadamard-test loops from $A|x(\theta)\rangle \propto |b\rangle$07 to $A|x(\theta)\rangle \propto |b\rangle$08, with runtime dropping from approximately $A|x(\theta)\rangle \propto |b\rangle$09 minutes to approximately $A|x(\theta)\rangle \propto |b\rangle$10 minutes (Yi et al., 2023). For quantum transport in nanoscale semiconductor devices, VQLS is extended to complex, non-Hermitian NEGF systems; an $A|x(\theta)\rangle \propto |b\rangle$11-qubit implementation solves a $A|x(\theta)\rangle \propto |b\rangle$12 block-diagonal multi-energy system and reconstructs LDOS, transmission $A|x(\theta)\rangle \propto |b\rangle$13, and current in excellent agreement with classical NEGF (Yang et al., 6 Sep 2025).

6. Scalability, distributed formulations, and unresolved issues

A recurrent misconception is that logarithmic qubit scaling alone implies practical quantum advantage. The broad empirical study of VQLS on larger and less structured instances identifies three main obstacles: ansatz expressibility, measurement overhead, and circuit overhead. In the most general Pauli-LCU setting, $A|x(\theta)\rangle \propto |b\rangle$14 can scale as $A|x(\theta)\rangle \propto |b\rangle$15, making cost evaluation $A|x(\theta)\rangle \propto |b\rangle$16 in circuit calls; for random Pauli systems at $A|x(\theta)\rangle \propto |b\rangle$17, denominator circuits require up to $A|x(\theta)\rangle \propto |b\rangle$18 gates, and for Darcy flow matrices the hardware-efficient ansatz with up to five layers reaches only $A|x(\theta)\rangle \propto |b\rangle$19 at best (Turati et al., 2024). The same study concludes that further research is necessary to fully leverage the algorithm’s capabilities in addressing real-world problems (Turati et al., 2024).

Recent work therefore targets distributed and measurement-efficient execution rather than only better ansätze. One distributed VQLS framework, implemented with NVIDIA CUDA-Q, asynchronously distributes the $A|x(\theta)\rangle \propto |b\rangle$20 overlap circuits across GPUs and reports near-ideal strong scaling up to $A|x(\theta)\rangle \propto |b\rangle$21 GPUs, $A|x(\theta)\rangle \propto |b\rangle$22 weak scaling efficiency at $A|x(\theta)\rangle \propto |b\rangle$23 GPUs processing $A|x(\theta)\rangle \propto |b\rangle$24 circuits per iteration for a $A|x(\theta)\rangle \propto |b\rangle$25-qubit system, and a $A|x(\theta)\rangle \propto |b\rangle$26 speedup for the studied configurations (Lu et al., 15 Apr 2026). A separate distributed formulation partitions a large matrix into square block submatrices, assigns each block to a NISQ computer, and combines local VQLS objectives with DIGing-style gradient tracking and Adam-type updates; numerical simulations show that four $A|x(\theta)\rangle \propto |b\rangle$27-qubit devices can cooperatively solve a $A|x(\theta)\rangle \propto |b\rangle$28-qubit problem (Shen et al., 1 Apr 2026). Measurement-efficient VQLS for Carleman-linearized Duffing dynamics reports near-unity fidelity and vanishing relative residuals across block-banded test cases by combining topology-agnostic ansätze, Hermitianization choices, local cost functions, and symmetry-grouped Hadamard tests (Liu et al., 14 May 2026).

The resulting picture is technically specific rather than uniform. On structured sparse systems, especially those arising from banded PDE discretizations, VQLS can achieve high fidelities with shallow circuits and carefully engineered decompositions (2207.14630, Balducci et al., 2024, Lu et al., 15 Apr 2026). On general or weakly structured systems, measurement cost, state preparation, and ansatz expressibility remain dominant constraints (Turati et al., 2024). This suggests that VQLS is best understood not as a single solver with a fixed scaling law, but as a variational linear-algebra framework whose practicality depends on matrix structure, cost design, decomposition strategy, optimizer robustness, and the architecture used to distribute or compress the required measurements.