PaperlayerView exact revision on arXiv

Source-linked AI summary

Well-conditioned iterative methods for large open quantum systems

Gaspard Beugnot, Paul Gregory, Rémi Robin, Antoine Tilloy

arXiv:2608.30860v1quant-phmath.NA

TL;DR

Large open quantum systems make exact linear-algebra approaches prohibitively expensive, motivating iterative methods with effective preconditioning. The paper efficiently inverts the no-jump evolution to construct a CPTP fixed-point method and precondition spectral and time-integration computations, reaching larger systems and strong CPU/GPU performance.

  • Problem

    Exact diagonalization and exponentiation scale as O(n^6), becoming prohibitive for systems with n ≳100, while iterative methods require effective preconditioning.

  • Method

    The paper inverts the no-jump part S via a continuous Lyapunov equation, then uses Φ = −KS^-1 for steady states and S^-1 to precondition spectral and implicit time-integration methods.

  • Results

    The methods manipulate n×n matrices without forming the n^2 × n^2 Liouvillian, scale as O(n^3) per iteration, and compute spectra roughly one order of magnitude larger than SciPy baselines.

  • Takeaways & Limitations

    The auxiliary CPTP map makes Arnoldi iterations effective for steady-state computation, while the shared preconditioner supports low-lying spectra and stiff implicit dynamics with gains amplified on GPUs.

  • Takeaways & Limitations

    The auxiliary map can have non-trivial peripheral eigenvalues, and optimizing the jump/no-jump split is left for future work.

Abstract

from arXiv · show

Markovian open quantum systems are well modeled by the Lindblad Master Equation (ME) $\frac{\mathrm{d}}{\mathrm{d} t} ρ_t = \mathcal{L} ρ_t$, where $\mathcal{L}$ is a linear (super-)operator and $ρ_t$ is the system state, a positive matrix. When designing or characterizing a quantum system, one is usually interested in the steady state $ρ_\infty$ (such that $\mathcal{L} ρ_\infty = 0$), the first few excited states, and trajectories $t\mapsto ρ_t$. In finite dimension, $ρ_t$ is an $n\times n$ matrix, $\mathcal{L}$ thus typically costs $n^4$ to store explicitly as a dense matrix, and $O(n^6)$ to diagonalize or invert exactly, making standard linear algebraic techniques expensive for large systems. However, $\mathcal{L}$ usually costs only $O(n^3)$ to apply. This makes iterative methods appealing, but they do not work without a good preconditioner. In this article, our main observation is that a part of the Lindblad equation, corresponding to the so-called no-jump evolution $\mathcal{S}$, can be inverted efficiently. Using this inverse map, we introduce an auxiliary completely positive trace-preserving (CPTP) map $Φ$ whose fixed point is directly related to $ρ_\infty$, all the other eigenvalues having smaller magnitude. The map $Φ$ is thus well suited to iterative methods, and $ρ_\infty$ can be found in a few Arnoldi iterations. Using the same inverse map $\mathcal{S}^{-1}$ as preconditioner, we compute the low-lying spectrum efficiently via shift-invert Arnoldi, and, as a proof of concept, build an implicit time integrator that is competitive on stiff systems in the low-precision regime. For the steady-state and low excited states problems, our methods scale like $O(n^3)$ per iteration and offer state-of-the-art performance on CPU and GPU.

1 Introduction

The introduction frames steady-state, spectral, and real-time questions as central problems in modeling open quantum systems, but standard linear algebra becomes prohibitive at large dimension. It motivates iterative methods using the Lindbladian’s cheaper application and an efficiently invertible no-jump component, while noting that the decomposition itself can be chosen in multiple ways.

  • Motivations and scope of the paper: Steady states describe long-time behavior, support experimental calibration, and enter expectation values for continuous matrix product states.Efficient steady-state computation is therefore relevant both to physical-system characterization and to larger tensor-network calculations.
  • Motivations and scope of the paper: The spectral gap and first few Lindbladian eigenvalues determine the speed of convergence to the steady state.
  • Motivations and scope of the paper: Real-time computation seeks an approximation to the master-equation solution for a specified initial state.
  • Motivations and scope of the paper: O(n^6) standard linear algebra on the n^2-dimensional operator space becomes prohibitive for n ≳100.This cost covers exact diagonalization and exponentiation, including the three main computational questions.
  • Motivations and scope of the paper: When the number of dissipators is O(1), applying the Lindbladian costs O(n^3), motivating Krylov-based iterative techniques.Direct iteration can nevertheless disappoint because the Lindbladian is often ill-conditioned, especially for discretizations of initially unbounded operators.
  • Methods overview: The no-jump generator can be inverted through a continuous Lyapunov equation in O(n^3), providing the paper’s preconditioning building block.The inverse is completely positive under the stated spectral assumption, and the no-jump/jump split determines which part is inverted exactly.
  • Methods overview: The Lindbladian decomposition is not unique: shifting jump operators by identity multiples changes the split while preserving an admissible jump map.The paper uses one fixed split and leaves optimization over this gauge freedom for future work.

2 Computing the steady state as the fixed point of CPTP maps

The paper reviews steady-state solvers and reformulates the problem using an auxiliary CPTP map whose fixed point yields the steady state. This approach avoids the poor scaling and conditioning of direct and unpreconditioned iterative methods.

  • The steady state can be computed by diagonalizing the Lindbladian, solving a linear system, or simulating long-time dynamics.
  • Iterative Arnoldi methods target low-magnitude eigenvalues inefficiently unless shift-invert is used, which requires an effective preconditioner.
  • Dense direct solvers construct an n^2 × n^2 Liouvillian and scale as O(n6).
  • The proposed map Φ = −KS−1 is CPTP, has a fixed point related to ρ∞, and can be applied at cost O(n3).When S is invertible, the steady state is obtained by normalizing −S−1ξ∞.
  • Arnoldi is effective because Φ has 1 as a simple eigenvalue of maximal modulus, unlike L whose target eigenvalue has smallest magnitude.Other peripheral eigenvalues can exist, so the Ritz value closest to 1 is retained rather than simply the dominant one.
  • Benchmarks show Arnoldi on Φ and preconditioned GMRES outperform competing solvers and gain further speedup on GPU.The comparison includes dense and cat-qubit benchmarks, with convergence assessed using ∥L(ρ)∥max < 10−8.

3 A new preconditioner to compute the resolvent

The paper uses the efficiently invertible no-jump resolvent as a right preconditioner for shifted Lindbladian systems. This yields contraction-based convergence guarantees and supports resolvent computation, steady-state solving, spectral calculations, and implicit time integration without forming the full Liouvillian.

  • 3 A new preconditioner to compute the resolvent: The no-jump resolvent R_S^λ is proposed as an inexpensive right preconditioner for λ−L and its deflated variant.It is obtained from the no-jump evolution and is defined for λ>0, with the λ=0 case recovering the no-jump generator under the stated invertibility assumption.
  • 3 A new preconditioner to compute the resolvent: The preconditioned error operator is a strict contraction for every λ>0, enabling a Neumann-series representation of the Lindbladian resolvent.The resulting resolvent belongs to the resolvent set of L, and the series is closely connected to the Krylov space explored by the iterative solver.
  • 3 A new preconditioner to compute the resolvent: Preconditioned GMRES converges at least geometrically with ratio ∥N∥<1, using a norm bound rather than spectral information and therefore remaining insensitive to non-normality.The Krylov space contains truncated Neumann iterates, while GMRES selects a residual-minimizing approximation.
  • 3 A new preconditioner to compute the resolvent: At zero shift, strict contraction is not established for the steady-state and shift-invert solvers, although GMRES converges in a few iterations in the reported experiments.Deflation makes the linear systems nonsingular, but the theory does not show that it restores strict contraction.
  • 3 A new preconditioner to compute the resolvent: At λmax/100 on the cat family, the method gains a factor ∼50, while the matrix-free solver reaches n=3000 on dense systems as competitors drop out near n≈100.The advantage increases as the shift approaches zero and is amplified on GPU; at the largest shift, unpreconditioned GMRES already converges quickly.

4 Applications

The paper applies the no-jump inverse as a preconditioner to steady-state solves, low-lying spectral computations, and implicit time integration. Deflation and shift-invert Arnoldi target the relevant Lindbladian modes without explicitly forming the n^2 × n^2 operator.

  • The no-jump resolvent preconditioner is applied to steady-state linear solves, low-lying-spectrum computation, and implicit time integration.These are the three applications enabled by the preconditioner.
  • Steady state as a preconditioned linear solve: Preconditioned GMRES solves the deflated steady-state system using the no-jump resolvent at λ = 0, with rank-one deflation folded in by Sherman–Morrison.The system is eLη ρ = η Id, equivalently (0 − eLη)ρ = −η Id.
  • New method based on a preconditioned shift-invert: The low-lying spectrum is computed by applying Arnoldi to (µ − eLη)^−1, which removes the zero mode, makes targeted eigenvalues dominant, and uses inexpensive preconditioned inner solves.The deflated operator preserves nonzero eigenpairs while shifting the steady-state eigenvalue to ηn.
  • New method based on a preconditioned shift-invert: The solver reaches systems roughly one order of magnitude larger than SciPy baselines because it never forms the n^2 × n^2 Liouvillian.The baselines form that matrix and run out of memory beyond a few hundred basis states.
  • Implicit time integration: The implicit Euler scheme takes a roughly constant number of steps as Kerr strength increases, producing a large speedup once the system is sufficiently stiff.The scheme is compared with the explicit Tsit5 integrator.

5 Conclusion and perspectives

The paper’s methods share exact inversion of the no-jump part, enabling steady-state, spectral, and stiff-dynamics solvers that operate on n × n matrices. The demonstrated applications include noise-bias recovery in a cat qubit and speedups for stiff time evolution, while split selection and higher-order time integration remain open directions.

  • The no-jump part S is inverted exactly through a continuous Lyapunov equation at cost O(n^3), despite possible ill-conditioning of L.This inversion is the common building block for all three solvers.
  • The CPTP map Φ = −KS^−1 has a simple maximal-modulus eigenvalue 1 whose fixed point yields the steady state.This makes the fixed-point problem suitable for Krylov methods.
  • Low-lying spectrum: The cat-qubit spectrum shows an exponentially suppressed bit-flip rate and phase-flip rates that grow linearly with cat size |α|^2.Residual single-photon loss produces one steady state and three slow modes.
  • The resulting solvers manipulate n × n matrices while supporting steady-state computation, shift-invert spectral calculation, and implicit stiff-dynamics integration.All three tasks use the same no-jump-based machinery.
  • Open directions: The present work does not optimize the gauge freedom in the jump/no-jump split, although such optimization is expected to improve quantitative performance.The authors leave this direction to future work.
  • Implicit time integration: The implicit Euler scheme gains speed on sufficiently stiff systems, but adaptive step sizing and higher-order methods are needed for practical competitiveness.The current integrator is presented as a proof of concept.

A Numerical implementation

The numerical implementation builds all three solvers from exact inversion of the no-jump operator via continuous Lyapunov equations. It packages this inversion into a no-jump resolvent and the induced right preconditioner.

  • All three main-text solvers rely on exact inversion of S, implemented by solving a continuous Lyapunov equation.The steady-state and low-lying-spectrum solvers are built on this shared operation.
  • The implementation introduces the no-jump resolvent R_S^λ and the right preconditioner it induces as the two reusable algorithms.These algorithms are leveraged by the three applications.
  • The shifted Lyapunov equation (λ − S)(X) = Y supplies the resolvent used by the preconditioner.The unshifted equation is S(X) = GX + XG† = Y.

A.1.1 Bartels–Stewart algorithm

The Bartels–Stewart implementation solves the Lyapunov equation by Schur decomposition and triangular back-substitution. Hermitian right-hand sides reduce cost, while sequential back-substitution limits GPU parallelism and can require refinement for stability.

  • Bartels–Stewart computes the Schur decomposition G = UTU† and transforms the Lyapunov equation into a triangular equation solved by back-substitution.The transformed variables are X̃ = U†XU and Ỹ = U†YU.
  • Hermitian Y produces Hermitian X, allowing the implementation to halve the computational cost.This exploits the symmetry of the Lyapunov solution.
  • The shift is applied to the triangular factor as T → T − λ/2 Id, realizing −R_S^λ.This produces the shifted resolvent needed by the solver.
  • Sequential back-substitution is difficult to parallelize on GPUs and may have stability issues, which mixed-precision refinement can mitigate.The limitation concerns the numerical implementation of the Lyapunov solve.

A.1.2 Eigendecomposition method

When G is diagonalizable, its eigendecomposition provides a basis that diagonalizes the superoperator S and enables element-wise solution of the shifted equation.

  • The method diagonalizes G as G = UΣU−1 and transforms X and Y using V = U−† for element-wise decoupling.The shifted equation is solved in the transformed basis.
  • The inverse map is recovered at λ = 0 as −S−1.
  • Because ℜ(σi) < 0 when S is invertible, denominators remain nonzero for λ ⩾ 0.
  • The diagonal basis enables a fast matrix-matrix implementation of the inverse.

A.1.3 Comparison and implementation

The eigendecomposition and Bartels–Stewart approaches both scale cubically, but their hardware behavior differs: GEMM favors GPUs, whereas Schur back substitution lacks an efficient GPU implementation.

  • Both the eigendecomposition and Bartels–Stewart methods scale as O(n^3).
  • Bartels–Stewart has the smaller leading constant.
  • The eigendecomposition route uses matrix–matrix products that are extremely optimized on GPUs.
  • Bartels–Stewart relies on Schur factorization and back substitution, for which the authors did not find an efficient GPU implementation.
  • Algorithm 1 precomputes the eigendecomposition G = UΣU−1 and stores U and V at O(n^3) cost.

A.2 Linear solver implementation

The implementation combines the efficiently inverted no-jump operator with recycled, right-preconditioned Krylov methods for steady states and low-lying Lindbladian modes.

  • Linear solver implementation: Every GMRES solve uses the no-jump resolvent as a right preconditioner without forming the n^2 × n^2 Liouvillian.
  • Linear solver implementation: The custom matrix-free JAX GMRES applies the Lindbladian directly and supports Krylov-subspace recycling across restarts.
  • Deflation: The deflated systems add a rank-one term and incorporate it exactly into the preconditioner using the Sherman–Morrison formula.The correction requires one extra base solve and a trace, and excludes a single η value far from those used.
  • Steady state: The steady state can be computed either as the fixed point of Φ or from a deflated linear system.
  • Solver assembly: The implementation uses recycled GMRES for steady-state solves and shift-invert Arnoldi for the low-lying spectrum, while Scipy’s default ARPACK is slower by a constant factor.
  • Steady state: Φ has dominant eigenvalue 1, and restarted Arnoldi selects the Ritz vector closest to 1 before reconstructing and testing the state.
  • Low-lying spectrum: Shift-invert Arnoldi targets the smallest-magnitude nonzero modes at μ = 0 by solving systems involving the deflated Lindbladian.
  • Low-lying spectrum: Residuals are measured on the true Lindbladian, using maxk ∥Lvk −λkvk∥max < ε for raw, traceless eigenvectors.

B.1 Description of the system

The benchmark system is a truncated memory–buffer model for cat-state stabilization, complemented by a generic dense random benchmark; truncations are chosen to resolve populated states and avoid artifacts.

  • Sparse benchmark: The sparse benchmark models cat-state stabilization in a memory–buffer architecture with memory and buffer creation and annihilation operators.
  • Sparse benchmark: The model captures two-photon exchange with a lossy buffer plus thermal excitation and relaxation, without a Kerr term.
  • Sparse benchmark: The resulting Lindbladian is sparse, has a small spectral gap that shrinks with cat size, and has a unique steady state.
  • Truncation: Infinite-dimensional operator algebra is approximated by truncating the memory and buffer modes to na and nb levels.
  • Truncation: The effective Hilbert-space dimension is n = na nb, with scaling choices that keep the target cat amplitude within the truncated memory space.The stated truncation choices are intended to resolve physically populated Fock states while avoiding truncation artifacts.
  • Dense benchmark: The dense benchmark draws generic unstructured matrices with independent standard-normal real and imaginary parts.
  • Dense benchmark: A single random stream is drawn in increasing n so the dense systems are reproducible, and all generators are dense.

B.2 Description of the competing methods

The evaluation compares the authors’ fixed-point Arnoldi and preconditioned GMRES methods with dense, sparse-direct, and unpreconditioned iterative baselines across steady-state, spectral, resolvent, and time-evolution tasks. The experiments use specified residual or fidelity criteria, while the implicit method is tested as a low-precision proof of concept on systems with tunable stiffness.

  • Steady state: Steady-state baselines include dense SVD, sparse LU with MUMPS, and CPU GMRES with or without incomplete-LU preconditioning.The dense SVD method scales as O(n^6), while the sparse-direct and iterative alternatives operate on trace-pinned systems.
  • Authors’ methods: The authors use fixed-point Arnoldi for the steady state and preconditioned GMRES for linear solves, both with an eigendecomposition Lyapunov backend.The steady-state evaluation includes the fixed-point Arnoldi method and the preconditioned GMRES linear solver.
  • Low-lying spectrum: For the low-lying spectrum, dense NumPy eigendecomposition and ARPACK shift-invert provide baselines against the authors’ shift-invert Arnoldi method.The baselines return four eigenvalues including the steady state; the authors’ method returns the three smallest nonzero modes.
  • Time evolution: On a quartic-stiffness sweep, implicit Euler uses one preconditioned GMRES cycle per step, beats Tsit5 above K ≈2.5, and reaches ∼34× speedup at K = 15.Its step count is essentially independent of K, but the advantage applies only in the low-precision regime; tighter tolerances make it uncompetitive.
  • Resolvent: Resolvent experiments compare dense LU, MUMPS, unpreconditioned GMRES, and GMRES with incomplete-LU preconditioning against recycled preconditioned GMRES.The shifted systems use λ values based on the largest-magnitude eigenvalue, with convergence assessed by a max-norm residual below 10^-8.
Loading 2608.30860v1…