PaperlayerView exact revision on arXiv

Source-linked AI summary

Structure-preserving full- and low-rank exponential methods for the optimal control of Lindblad equations

Hao Chen, Alfio Borzi

arXiv:2610.09403v1math.NAcs.NA

TL;DR

Open-system optimal control requires repeated forward and adjoint Lindblad propagation that is accurate, structure-preserving, and efficient at large scale. The paper combines continuous Pontryagin optimality conditions and SQH with full- and low-rank exponential midpoint propagators and rigorous error estimates. Numerical experiments report lower computational time for SQH implementations than GRAPE in tested cases, with additional CPU-time reduction for the low-rank implementation.

  • Problem

    Repeated forward and adjoint propagation makes efficient, accurate, structure-preserving computation challenging for large Hilbert spaces, while optimal control must accommodate nonsmooth costs.

    Source: §1 · pp. 1–2

  • Method

    The framework derives continuous Pontryagin conditions, applies SQH, and discretizes forward and adjoint Lindblad equations with full- and low-rank exponential midpoint propagators.

    Source: §1 · p. 2

  • Results

    For tested discretizations and system dimensions, SQH-FREM and SQH-LREM required less computational time than GRAPE, and low-rank implementation further reduced CPU time versus full-rank realization.

    Source: §8.2 · p. 28

  • Takeaways & Limitations

    The experiments demonstrate the framework's effectiveness for tested smooth, nonsmooth, and nonconvex open-system control problems.

    Source: §8.2 · pp. 28–29

    The Pontryagin and SQH analysis is restricted to coherent Hamiltonian control, although the exponential midpoint discretization can accommodate incoherent controls.

    Source: §2 · p. 3

Abstract

from arXiv · show

Optimal control of open quantum systems governed by the Lindblad master equation requires the repeated solution of forward and adjoint evolution equations. For large Hilbert spaces, these computations demand numerical methods that are structure-preserving, accurate, and computationally efficient, while remaining compatible with nonsmooth optimization techniques. In this work, a unified framework for optimal control of Lindblad equations is presented. The continuous optimality system is derived from the Pontryagin maximum principle and solved by a sequential quadratic Hamiltonian (SQH) method. Its numerical realization is based on second-order exponential midpoint propagators for the forward and adjoint Lindblad equations with time-dependent Hamiltonians. Full-rank schemes preserving the Hermitian and positive-semidefinite structure of the density matrix are developed together with low-rank formulations that substantially reduce storage requirements and computational cost. Rigorous error estimates are established for the full- and low-rank forward and adjoint propagators. Numerical experiments confirm the predicted convergence rates, demonstrate the effectiveness of the low-rank approximations, and illustrate the performance of the proposed FREM-SQH and LREM-SQH algorithms for optimal control problems with smooth and nonsmooth control costs.

1. Introduction

Open-system optimal control requires repeated forward and adjoint Lindblad propagation, so the propagators must combine physical structure preservation, accuracy for time-dependent Hamiltonians, and efficiency at large scale. The paper addresses these needs with an optimize-then-discretize framework, structure-preserving full- and low-rank exponential midpoint methods, and rigorous error estimates.

Motivation

  • Repeated forward and backward-adjoint integration makes propagator quality central to optimization efficiency, while large systems also demand physical structure preservation and computational efficiency.For time-dependent Hamiltonians, high accuracy is another stated requirement.Source: §1 · pp. 1–2
  • Ten qutrits require dense 59 049×59 049 matrices occupying approximately 52 GiB each, illustrating the scaling pressure motivating low-rank methods.For n_q identical d-level qudits, the Hilbert-space dimension is m = d^n_q, density and adjoint matrices each have m^2 real degrees of freedom, and dominant operations scale as O(m^3).Source: §1 · p. 2

Related work

  • Existing methods are often limited to time-independent Hamiltonians or first-order accuracy for time-dependent cases, and rigorous low-rank convergence analysis remains limited.The introduction identifies unavailable second-order error estimates and low-rank adjoint convergence analysis as gaps.Source: §1 · p. 2

Contribution

  • The framework derives continuous optimality conditions by the Pontryagin maximum principle and uses an optimize-then-discretize strategy to retain matrix structure in state and adjoint discretizations.This contrasts with approaches that discretize first and then optimize.Source: §1 · p. 2
  • Full-rank exponential midpoint propagators preserve density-matrix Hermiticity and positive semidefiniteness, while low-rank versions reduce storage and computational cost.The work also establishes rigorous error estimates for the resulting forward and adjoint propagators.Source: §1 · p. 2
  • The SQH framework is extended to Lindblad dynamics and accommodates smooth, nonsmooth, and nonconvex control costs.The paper's numerical experiments assess the FREM-SQH and LREM-SQH algorithms and validate theoretical convergence results.Source: §1 · p. 3

2. Controlled Lindblad dynamics and optimal control

The section formulates coherent Hamiltonian control of Lindblad dynamics as an optimal control problem over density matrices. Its running cost need not be differentiable or convex, allowing the framework to encompass sparsity-promoting and nonconvex control costs.

Controlled dynamics

  • The Lindblad equation serves as the state dynamics, with the system represented by a density matrix that is Hermitian, positive semidefinite, and unit trace.The equation models Markovian open-system dynamics.Source: §2 · p. 3
  • The formulation uses coherent control through a Hermitian Hamiltonian H0 + uH1, and the presentation's single real-valued control extends directly to multiple coherent controls.The authors retain one control to simplify notation.Source: §2 · p. 3

Optimal control problem

  • The running cost is bounded below and lower semicontinuous but need not be differentiable or convex, encompassing smooth quadratic regularization and sparsity-promoting or nonconvex costs.Examples include L1-, nonconvex Lp-, and L0-type costs.Source: §2 · p. 4
  • The paper extends SQH methods for nonsmooth and sparsity-promoting quantum control to open systems governed by Lindblad dynamics.The terminal objectives include observable and quadratic tracking functionals; the latter can have a Hermitian but indefinite terminal derivative relevant to low-rank adjoint approximation.Source: §2 · p. 4
  • The continuous optimality conditions are derived before discretization, yielding forward and adjoint equations plus pointwise Hamiltonian minimization for later full- or low-rank exponential midpoint realization.This is the section's optimize-then-discretize ordering.Source: §2 · p. 4

3. Pontryagin optimality system

The section derives continuous first-order optimality conditions for the controlled Lindblad problem using the Pontryagin framework. The resulting system couples forward state evolution, backward adjoint evolution, and pointwise minimization of the Hamilton-Pontryagin function.

Adjoint equation

  • The adjoint of the controlled Lindblad generator includes Hamiltonian and dissipative contributions and evolves backward from the terminal time.The adjoint evolution preserves Hermiticity but generally does not preserve positive semidefiniteness.Source: §3 · p. 5

Hamilton-Pontryagin function

  • Because the control enters the Hamiltonian linearly, the Hamilton-Pontryagin function is affine in the control and its control-dependent part combines the running cost with a switching term.The switching quantity is real-valued for Hermitian Hamiltonian and state matrices.Source: §3 · pp. 5–6

Pontryagin optimality system

  • The optimality system combines the controlled state equation, an adjoint equation with terminal condition from the terminal objective, and the Hamiltonian stationarity condition.The optimality conditions are derived continuously before either equation is discretized.Source: §3 · pp. 5–7
  • The Pontryagin condition characterizes an optimal control through pointwise minimization of the Hamilton-Pontryagin function, including when the running cost is nonsmooth.The paper identifies this characterization as the basis for extending SQH to Lindblad dynamics.Source: §3 · pp. 6–7

4. Sequential quadratic Hamiltonian method

The SQH iteration repeatedly solves the state and adjoint equations and minimizes a stabilized Hamilton-Pontryagin function pointwise. An adaptive stabilization loop accepts sufficient-decrease steps and supports monotone objective decrease under the stated assumptions.

4. Sequential quadratic Hamiltonian method

  • The control update consists of independent finite-dimensional minimizations at each time once the current state and adjoint are computed.For common running costs, these minimizations may be solved analytically; otherwise, inexpensive local procedures can be used.Source: §4 · pp. 7–8
  • SQH uses a quadratic stabilization penalty to moderate control updates, while retaining pointwise minimization for smooth, nonsmooth, and nonconvex costs.Unrestricted pointwise minimization can produce large control changes and fail to decrease the objective; the stabilized formulation addresses this robustness issue.Source: §4 · pp. 7–8
  • Under the stated regularity assumptions, the SQH update satisfies a sufficient-decrease estimate, and accepted iterates have monotonically decreasing objective values that converge.The asymptotic regularity conclusion also follows under the accepted-step condition.Source: §4 · pp. 9–10

4.1. Adaptive choice of the stabilization parameter

  • The sufficient-decrease result guarantees that the adaptive inner loop accepts a trial after finitely many stabilization increases.Each rejection increases the stabilization parameter; after acceptance, the next initial trial value is reduced by a factor between zero and one.Source: §4 · p. 10, §4.1 · p. 11
  • The accepted stabilization parameters are bounded by the initial value and a bound involving the increase factor and the sufficient-decrease constants.The stated bound is ε_k < σ(2C_SQH + 2η) for a rejected trial followed by acceptance.Source: §4.1 · p. 11
  • The algorithm discretizes its continuous state and adjoint equations with second-order full- and low-rank exponential midpoint schemes, producing FREM-SQH and LREM-SQH.The SQH optimization procedure itself is independent of the particular time discretization.Source: §4 · p. 7, §4.1 · p. 12

5. Exponential midpoint propagation of the optimality system

The optimality system is discretized with exponential midpoint propagators for the forward and adjoint Lindblad equations. Full-rank updates preserve matrix structure, while low-rank variants work on factors to reduce computational and storage costs.

5.1. Full-rank midpoint propagators

  • The FREM schemes use quadrature approximations to derive exponential midpoint updates for both the forward and adjoint Lindblad equations.The forward update uses left-rectangle and midpoint quadrature; the adjoint update uses right-rectangle and midpoint quadrature.Source: §5.1 · pp. 13–14
  • Full-rank forward and adjoint updates preserve Hermiticity and positive semidefiniteness for every step size when their initial matrices are positive semidefinite.Each update is represented as a sum of congruence transformations.Source: §5.1 · p. 14
  • Because the unnormalized forward FREM scheme may not preserve unit trace, a normalized variant is proposed to restore unit trace.The normalization procedure is not adapted to the adjoint equation, which is usually not trace-preserving.Source: §5.1 · p. 14

5.2. Low-rank exponential integrators

  • Forward and adjoint low-rank solutions use different factorizations because density matrices are positive semidefinite, whereas adjoint matrices may be indefinite.The low-rank schemes approximate the solutions using factors whose ranks are much smaller than the Hilbert-space dimension.Source: §5.2 · pp. 14–15
  • Low-rank propagators replace full density matrices with factors, reducing storage from O(m^2) to O(mr) when r is much smaller than m.The matrix-exponential products are also applied to m × r factors rather than full density matrices.Source: §5.2 · pp. 14, 17

6. Error analysis of the forward problem

The forward error analysis establishes estimates for the full-rank FREM and low-rank LREM schemes. The LREM global error bound combines time-step, initial approximation, column-compression, and matrix-exponential approximation errors.

6.1. Error estimate of the FREM scheme

  • A single FREM step has trace within Cϱτ^3 of one, supported by a local error bound of Cϱτ^3 and positivity preservation.The local bound is in trace norm; the trace estimate is stated for n = 0, …, N−1.Source: §6.1 · p. 18
  • The analysis proves global error estimates for both unnormalized and normalized FREM schemes.The theorem statements identify the numerical solutions generated by schemes (49) and (54), respectively.Source: §6.1 · pp. 18–19

6.2. Error estimate of the LREM scheme

  • The LREM global error is bounded by c1 τ^2 + c2 δ + c3 ε1 + c4 ε2 under the stated error-tolerance assumptions.The bound applies for 0 ≤ n ≤ N; the tolerances satisfy ε1 = τϵ1 and ε2 = τϵ2.Source: §6.2 · p. 21

7. Error analysis of the adjoint problem

The adjoint error analysis establishes second-order convergence for FREM and bounds LREM error through discretization and approximation terms. For LREM, second-order accuracy requires suitably scaled terminal, compression, and exponential-approximation errors.

7.1. Error estimate of the FREM scheme

  • The FREM adjoint scheme has second-order convergence.Source: §7.1 · p. 23

7.2. Error estimate of the LREM scheme

  • For 0≤n≤N, ‖q(tₙ)−pₙ‖₁ ≤ c₁τ²+c₂δ+c₃ε₁+c₄ε₂; ε₁ and ε₂ are column-compression and matrix-exponential error tolerances, respectively.The theorem assumes ε₁=τϵ₁ and ε₂=τϵ₂ for some ϵ₁,ϵ₂>0.Source: §7.2 · p. 25
  • The complete LREM approximation is second-order when δ=O(τ²), ε1=O(τ³), and ε2=O(τ³).Source: §7.2 · p. 26

8. Numerical experiments

The experiments test the full- and low-rank propagators’ convergence, structure preservation, and computational performance, then apply them in SQH optimal control for smooth and nonsmooth costs.

8.1. Verification of the full- and low-rank propagators

  • Full-rank forward and adjoint propagators showed second-order convergence, and the forward scheme preserved positivity and unit trace.The error plots cover both equations; population and trace-error plots illustrate structure preservation.Source: §8.1 · p. 27, Figure 8.1 · p. 30, Figure 8.2 · p. 30
  • Low-rank schemes recovered second-order convergence when approximation and compression errors stayed below temporal error; when these errors dominated, they limited accuracy, while the forward scheme preserved positivity and trace.The experiments vary initial approximation, column-compression, and matrix-exponential tolerances.Source: §8.1 · p. 27, Figure 8.3 · p. 31, Figure 8.4 · p. 31
  • At comparable accuracy, the proposed low-rank propagators substantially reduced computational cost against QuTiP, particularly as Hilbert-space dimension increased.The comparison used QuTiP’s mesolve with the dop853 integrator on the vectorized equation.Source: §8.1 · p. 27, Figure 8.7 · p. 32, Figure 8.8 · p. 33

8.2. Numerical results for optimal control

  • For the smooth control problem, SQH-FREM and SQH-LREM matched GRAPE’s controls and population transfers, decreased the objective monotonically, and ran faster, with LREM further reducing CPU time.These comparisons apply to the tested discretizations and system dimensions.Source: §8.2 · p. 28, Figure 8.9 · p. 33, Figure 8.10 · p. 34
  • In two nonsmooth, nonconvex control regimes, SQH produced monotone objective decrease and controls that steered population toward the target state.GRAPE was omitted because it relies on differentiability of the running cost.Source: §8.2 · p. 29, Figure 8.11 · p. 34, Figure 8.12 · p. 34

9. Conclusions

The paper presents a unified Lindblad optimal-control framework combining a Pontryagin-derived optimality system solved by SQH with structure-preserving full- and low-rank exponential midpoint propagators. It establishes error estimates and reports convergence and computational advantages for smooth and nonsmooth problems, while identifying extensions for future work.

  • The framework combines a Pontryagin maximum-principle optimality system solved by SQH with structure-preserving full- and low-rank exponential midpoint propagators for forward and adjoint equations.Source: §9 · p. 35
  • Rigorous error estimates cover full- and low-rank propagators; experiments confirm predicted convergence and show computational advantages for smooth and nonsmooth control problems.Source: §9 · p. 35
  • Further work includes adaptive timestepping, improved low-rank compression, broader open-system classes, and parallel implementations for high-dimensional applications.Source: §9 · p. 35
Loading 2610.09403v1…