From 135f1654cfee92534451fac0174f49239a3185c1 Mon Sep 17 00:00:00 2001 From: Malte Krug Date: Wed, 4 Feb 2026 02:51:08 +0100 Subject: [PATCH 01/13] Add ExponentialUtilities GRAPE support + tests --- Project.toml | 6 + docs/Project.toml | 1 + docs/make.jl | 1 + docs/src/methods.md | 47 ++++ ...antumPropagatorsExponentialUtilitiesExt.jl | 201 ++++++++++++++++++ ...PropagatorsQuantumGradientGeneratorsExt.jl | 11 + src/newton.jl | 14 +- test/Project.toml | 1 + test/runtests.jl | 10 + test/test_grape_exponentialutilities.jl | 30 +++ test/test_grape_exputils_liouvillian.jl | 111 ++++++++++ test/test_prop_interfaces.jl | 186 ++++++++++------ test/test_propagate.jl | 3 + 13 files changed, 546 insertions(+), 76 deletions(-) create mode 100644 ext/QuantumPropagatorsExponentialUtilitiesExt.jl create mode 100644 ext/QuantumPropagatorsQuantumGradientGeneratorsExt.jl create mode 100644 test/test_grape_exponentialutilities.jl create mode 100644 test/test_grape_exputils_liouvillian.jl diff --git a/Project.toml b/Project.toml index d3c0501c..cb7666ff 100644 --- a/Project.toml +++ b/Project.toml @@ -16,20 +16,26 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" [weakdeps] +ExponentialUtilities = "d4d017d3-3776-5f7e-afef-a10c40355c18" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" +QuantumGradientGenerators = "a563f35e-61db-434d-8c01-8b9e3ccdfd85" RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" [extensions] +QuantumPropagatorsExponentialUtilitiesExt = "ExponentialUtilities" QuantumPropagatorsODEExt = "OrdinaryDiffEq" +QuantumPropagatorsQuantumGradientGeneratorsExt = "QuantumGradientGenerators" QuantumPropagatorsRecursiveArrayToolsExt = "RecursiveArrayTools" QuantumPropagatorsStaticArraysExt = "StaticArrays" [compat] ArrayInterface = "7.0" +ExponentialUtilities = "1.11" OffsetArrays = "1" OrdinaryDiffEq = "6.59" ProgressMeter = "1" +QuantumGradientGenerators = "0.1" Random = "1" RecursiveArrayTools = "3.12" SpecialFunctions = "1, 2" diff --git a/docs/Project.toml b/docs/Project.toml index c62945d7..905dd99b 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -6,6 +6,7 @@ DocInventories = "43dc2714-ed3b-44b5-b226-857eda1aa7de" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244" DocumenterInterLinks = "d12716ef-a0f6-4df4-a9f1-a5a34e75c656" +ExponentialUtilities = "d4d017d3-3776-5f7e-afef-a10c40355c18" LiveServer = "16fef848-5104-11e9-1b77-fb7a48bbb589" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" diff --git a/docs/make.jl b/docs/make.jl index b0533025..9830c8cb 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -3,6 +3,7 @@ using Pkg using Documenter using QuantumPropagators using QuantumPropagators: AbstractPropagator, set_t!, set_state! +import ExponentialUtilities # ensure ExponentialUtilities extension is loaded import OrdinaryDiffEq # ensure ODE extension is loaded using DocumenterCitations using DocumenterInterLinks diff --git a/docs/src/methods.md b/docs/src/methods.md index 0e02b1e7..f5cd1fbc 100644 --- a/docs/src/methods.md +++ b/docs/src/methods.md @@ -15,6 +15,7 @@ As discussed in the [Overview](@ref overview_approaches), time propagation can b We consider this especially in the piecewise-constant case (`pwc=true` in [`propagate`](@ref)/[`init_prop`](@ref)), which is required for the traditional optimization methods [GRAPE](https://juliaquantumcontrol.github.io/GRAPE.jl/stable/) and [Krotov](https://juliaquantumcontrol.github.io/Krotov.jl/stable/). In these propagations, the time-dependent generator ``\op{H}(t)`` is [evaluated](@ref QuantumPropagators.Controls.evaluate) to a constant operator ``\op{H}`` on each interval of the time grid. The analytical solution to the Schrödinger or Liouville equation is well known, and propagation step simply has to evaluate the application of the time evolution operator ``\op{U} = \exp[-i \op{H} dt]`` to the state ``|Ψ⟩``. The following methods are built in to `QuantumPropagators`: * [`ExpProp`](@ref method_expprop) – constructs ``\op{U}`` explicitly and then applies it to ``|Ψ⟩`` + * [`ExponentialUtilities`](@ref method_exponentialutilities) – applies ``\exp(\op{H} dt) |Ψ⟩`` using Krylov methods without explicitly forming ``\op{U}`` * [`Cheby`](@ref method_cheby) — expansion of ``\op{U} |Ψ⟩`` into Chebychev polynomials, valid if ``\op{H}`` has real eigenvalues * [`Newton`](@ref method_newton) – expansion of ``\op{U} |Ψ⟩`` into Newton polynomials, valid if ``\op{H}`` has complex eigenvalues (non-Hermitian Hamiltonian, Liouvillian) @@ -65,6 +66,52 @@ init_prop(state, generator, tlist, method::Val{:ExpProp}; kwargs...) * Comparing against another propagator +## [`ExponentialUtilities`](@id method_exponentialutilities) + +The method requires that the [ExponentialUtilities.jl](https://docs.sciml.ai/ExponentialUtilities/stable/) +package is loaded + +``` +using ExponentialUtilities +``` + +and then passed as `method=ExponentialUtilities` (or `method=:expv`) to +[`propagate`](@ref) or [`init_prop`](@ref): + +```@docs +init_prop(state, generator, tlist, method::Val{:ExponentialUtilities}; kwargs...) +``` + +This method evaluates ``\exp(-i \op{H} dt) |Ψ⟩`` via a Krylov expv algorithm +without explicitly forming the matrix exponential. It is therefore often a +good fit for larger systems or matrix-free operators where direct matrix +exponentiation is too costly. + +**Advantages** + +* Avoids explicit construction of ``\op{U}`` +* Works with matrix-free operators that support `mul!` +* Good scaling for large sparse systems + +**Disadvantages** + +* Requires ExponentialUtilities.jl +* Performance depends on Krylov parameters and operator structure + +**When to use** + +* Large, sparse, or matrix-free generators +* Piecewise-constant GRAPE-style workflows that need expv + +**GRAPE note** + +* When using `method=:expv` with GRAPE, set `gradient_method=:taylor`. + The default `:gradgen` path uses `GradVector`/`GradgenOperator`, which is + not currently supported by the ExponentialUtilities Krylov backend. +* For non-Hermitian generators (e.g., Liouvillians), pass + `prop_expv_kwargs=(; ishermitian=false)` to avoid Hermitian-only paths. + + ## [`Cheby`](@id method_cheby) The method should be loaded with diff --git a/ext/QuantumPropagatorsExponentialUtilitiesExt.jl b/ext/QuantumPropagatorsExponentialUtilitiesExt.jl new file mode 100644 index 00000000..da9c59dc --- /dev/null +++ b/ext/QuantumPropagatorsExponentialUtilitiesExt.jl @@ -0,0 +1,201 @@ +module QuantumPropagatorsExponentialUtilitiesExt + +using LinearAlgebra +using TimerOutputs: @timeit_debug, TimerOutput +using ExponentialUtilities + +using QuantumPropagators: + QuantumPropagators, + PWCPropagator, + _pwc_get_max_genop, + _pwc_get_genop, + _pwc_set_genop!, + _pwc_set_t!, + _pwc_advance_time!, + _pwc_process_parameters +using QuantumPropagators.Controls: get_controls +using QuantumPropagators.Interfaces: supports_inplace +using QuantumPropagators.Generators: ScaledOperator +import QuantumPropagators: init_prop, prop_step!, set_t! + +struct _SizedOp{A} + A::A +end + +Base.size(op::_SizedOp) = size(op.A) +Base.size(op::_SizedOp, dim::Integer) = size(op.A)[dim] +Base.eltype(op::_SizedOp) = eltype(op.A) +LinearAlgebra.ishermitian(op::_SizedOp) = LinearAlgebra.ishermitian(op.A) +LinearAlgebra.mul!(y, op::_SizedOp, x) = mul!(y, op.A, x) + +_ensure_size_dim(A) = + hasmethod(size, Tuple{typeof(A), Int}) ? A : + (hasmethod(size, Tuple{typeof(A)}) ? _SizedOp(A) : A) + + +"""Propagator for Krylov expv propagation via ExponentialUtilities (`method=ExponentialUtilities`). + +This is a [`PWCPropagator`](@ref). +""" +mutable struct ExpvPropagator{GT,OT,ST} <: PWCPropagator + const generator::GT + state::ST + t::Float64 # time at which current `state` is defined + n::Int64 # index of next interval to propagate + const tlist::Vector{Float64} + parameters::AbstractDict + controls + genop::OT + backward::Bool + inplace::Bool + convert_state::Type + convert_operator::Type + expv_kwargs::NamedTuple + const timing_data::TimerOutput +end + + +set_t!(propagator::ExpvPropagator, t) = _pwc_set_t!(propagator, t) + + +_expv_convert_state(state) = typeof(state) +_expv_convert_operator(::Any) = Any +_expv_convert_operator(::QuantumPropagators.Generator) = Matrix{ComplexF64} +_expv_convert_operator(::QuantumPropagators.Operator) = Matrix{ComplexF64} +_expv_convert_operator(::QuantumPropagators.ScaledOperator) = Matrix{ComplexF64} + + +""" +```julia +using ExponentialUtilities + +expv_propagator = init_prop( + state, + generator, + tlist; + method=ExponentialUtilities, + inplace=QuantumPropagators.Interfaces.supports_inplace(state), + backward=false, + verbose=false, + parameters=nothing, + expv_kwargs=(;), + convert_state=_expv_convert_state(state), + convert_operator=_expv_convert_operator(generator), + _... +) +``` + +initializes an [`ExpvPropagator`](@ref). + +# Method-specific keyword arguments + +* `expv_kwargs`: NamedTuple of keyword arguments forwarded to + `ExponentialUtilities.expv`. +* `convert_state`: Type to which to temporarily convert the state before + calling `expv`. +* `convert_operator`: Type to which to convert the operator before calling + `expv`. +""" +function init_prop( + state, + generator, + tlist, + method::Val{:ExponentialUtilities}; + inplace = supports_inplace(state), + backward = false, + verbose = false, + parameters = nothing, + expv_kwargs = (;), + convert_state = _expv_convert_state(state), + convert_operator = _expv_convert_operator(generator), + _... +) + tlist = convert(Vector{Float64}, tlist) + controls = get_controls(generator) + G = _pwc_get_max_genop(generator, controls, tlist) + + parameters = _pwc_process_parameters(parameters, controls, tlist) + timing_data = TimerOutput() + + n = 1 + t = tlist[1] + if backward + n = length(tlist) - 1 + t = float(tlist[n + 1]) + end + GT = typeof(generator) + OT = typeof(G) + ST = typeof(state) + return ExpvPropagator{GT,OT,ST}( + generator, + inplace ? copy(state) : state, + t, + n, + tlist, + parameters, + controls, + G, + backward, + inplace, + convert_state, + convert_operator, + expv_kwargs, + timing_data, + ) +end + + +# Aliases +init_prop(state, generator, tlist, method::Val{:expv}; kwargs...) = + init_prop(state, generator, tlist, Val(:ExponentialUtilities); kwargs...) + + +function prop_step!(propagator::ExpvPropagator) + @timeit_debug propagator.timing_data "prop_step!" begin + if nameof(typeof(propagator.state)) == :GradVector && + nameof(parentmodule(typeof(propagator.state))) == :QuantumGradientGenerators + throw(ArgumentError( + "ExponentialUtilities propagation does not support GRAPE `gradient_method=:gradgen`. " * + "Use `gradient_method=:taylor` instead." + )) + end + H = propagator.genop + n = propagator.n + tlist = getfield(propagator, :tlist) + (0 < n < length(tlist)) || return nothing + dt = tlist[n + 1] - tlist[n] + if propagator.backward + dt = -dt + end + dt_expv = complex(dt) + + Ψ = convert(propagator.convert_state, propagator.state) + if propagator.inplace + if supports_inplace(propagator.genop) + _pwc_set_genop!(propagator, n) + H = convert(propagator.convert_operator, propagator.genop) + else + H = convert(propagator.convert_operator, _pwc_get_genop(propagator, n)) + end + H = _ensure_size_dim(H) + H = ScaledOperator(-1im, H) + @timeit_debug propagator.timing_data "expv" begin + Ψ = ExponentialUtilities.expv(dt_expv, H, Ψ; propagator.expv_kwargs...) + end + copyto!(propagator.state, convert(typeof(propagator.state), Ψ)) + else + H = convert(propagator.convert_operator, _pwc_get_genop(propagator, n)) + H = _ensure_size_dim(H) + H = ScaledOperator(-1im, H) + @timeit_debug propagator.timing_data "expv" begin + Ψ = ExponentialUtilities.expv(dt_expv, H, Ψ; propagator.expv_kwargs...) + end + setfield!(propagator, :state, convert(typeof(propagator.state), Ψ)) + end + + _pwc_advance_time!(propagator) + return propagator.state + end +end + +end diff --git a/ext/QuantumPropagatorsQuantumGradientGeneratorsExt.jl b/ext/QuantumPropagatorsQuantumGradientGeneratorsExt.jl new file mode 100644 index 00000000..484d5b32 --- /dev/null +++ b/ext/QuantumPropagatorsQuantumGradientGeneratorsExt.jl @@ -0,0 +1,11 @@ +module QuantumPropagatorsQuantumGradientGeneratorsExt + +using QuantumGradientGenerators: GradgenOperator, GradVector + +Base.size(G::GradgenOperator, dim::Integer) = size(G.G, dim) +Base.similar(::GradVector, ::Type{T}, dims::Tuple{Int,Int}) where {T} = + Matrix{T}(undef, dims...) +Base.similar(::GradVector, ::Type{T}, dims::Tuple{Int}) where {T} = + Vector{T}(undef, dims[1]) + +end diff --git a/src/newton.jl b/src/newton.jl index 65f755b3..b5cafd97 100644 --- a/src/newton.jl +++ b/src/newton.jl @@ -34,7 +34,7 @@ mutable struct NewtonWrk{T} restarts::Int64 m_max::Int64 timing_data::TimerOutput - function NewtonWrk(v0::T; m_max::Int64 = 10, _timing_data = TimerOutput()) where {T} + function NewtonWrk(v0::T; m_max::Int64=10, _timing_data=TimerOutput()) where {T} if m_max <= 2 error("Newton propagation requires m_max > 2") end @@ -60,8 +60,8 @@ mutable struct NewtonWrk{T} end -lbound(array::OffsetArray, dim = 1) = first(axes(array)[dim]) -ubound(array::OffsetArray, dim = 1) = last(axes(array)[dim]) +lbound(array::OffsetArray, dim=1) = first(axes(array)[dim]) +ubound(array::OffsetArray, dim=1) = last(axes(array)[dim]) function leja_radius(z) @@ -281,9 +281,9 @@ function newton!(Ψ, H, dt, wrk; kwargs...) wrk.v, H, _dt; - extended = true, - norm_min = norm_min, - _timing_data = wrk.timing_data + extended=true, + norm_min=norm_min, + _timing_data=wrk.timing_data ) end if m == 1 && s == 0 @@ -294,7 +294,7 @@ function newton!(Ψ, H, dt, wrk; kwargs...) break end @timeit_debug wrk.timing_data "diagonalize_hessenberg_matrix" begin - ritz = diagonalize_hessenberg_matrix(Hess, m, accumulate = true) + ritz = diagonalize_hessenberg_matrix(Hess, m, accumulate=true) end # In the first iteration, the radius will be determined diff --git a/test/Project.toml b/test/Project.toml index 3df51a9c..fa8c5de6 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -9,6 +9,7 @@ DocInventories = "43dc2714-ed3b-44b5-b226-857eda1aa7de" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244" DocumenterInterLinks = "d12716ef-a0f6-4df4-a9f1-a5a34e75c656" +ExponentialUtilities = "d4d017d3-3776-5f7e-afef-a10c40355c18" GRAPE = "6b52fcaf-80fe-489a-93e9-9f92080510be" GRAPELinesearchAnalysis = "290eba36-e2d8-4488-81b6-f66cc44f2186" IOCapture = "b5f81e59-6552-4d32-b1f0-c071b021bf89" diff --git a/test/runtests.jl b/test/runtests.jl index cc3ac140..f2cbd52b 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -93,6 +93,16 @@ using SafeTestsets include("test_propagate_sequence.jl") end + println("\n* GRAPE with ExponentialUtilities (test_grape_exponentialutilities.jl):") + @time @safetestset "GRAPE ExponentialUtilities" begin + include("test_grape_exponentialutilities.jl") + end + + println("\n* GRAPE ExponentialUtilities Liouvillian (test_grape_exputils_liouvillian.jl):") + @time @safetestset "GRAPE ExponentialUtilities Liouvillian" begin + include("test_grape_exputils_liouvillian.jl") + end + println("\n* Time-dependent observables (test_timedependent_observables.jl):") @time @safetestset "Time-dependent observables" begin include("test_timedependent_observables.jl") diff --git a/test/test_grape_exponentialutilities.jl b/test/test_grape_exponentialutilities.jl new file mode 100644 index 00000000..e9801e0d --- /dev/null +++ b/test/test_grape_exponentialutilities.jl @@ -0,0 +1,30 @@ +using Test +using GRAPE +using QuantumControl +using QuantumControl.Functionals: J_T_sm +using QuantumControlTestUtils.DummyOptimization: dummy_control_problem +using StableRNGs: StableRNG +using ExponentialUtilities + +@testset "GRAPE with ExponentialUtilities" begin + rng = StableRNG(20260204) + problem = dummy_control_problem( + N=4, + n_trajectories=1, + n_controls=1, + n_steps=10, + dt=0.1, + density=1.0, + hermitian=true, + complex_operators=true, + prop_method=ExponentialUtilities, + J_T=J_T_sm, + iter_stop=1, + rng=rng + ) + + res = QuantumControl.optimize(problem; method=GRAPE, iter_stop=1, verbose=false) + + @test res.iter >= res.iter_start + @test isfinite(res.J_T) +end diff --git a/test/test_grape_exputils_liouvillian.jl b/test/test_grape_exputils_liouvillian.jl new file mode 100644 index 00000000..a3229caa --- /dev/null +++ b/test/test_grape_exputils_liouvillian.jl @@ -0,0 +1,111 @@ +using QuantumControl +using GRAPE +using LinearAlgebra +using QuantumPropagators.Controls: get_controls +using ExponentialUtilities +using Test + +# Scale +const GHz = 1.0 +const MHz = 0.001 * GHz +const ns = 1.0 + +function qubit_lvn(eps) + # Basis vectors + sigma_x = [0 1; 1 0] + sigma_z = [1 0; 0 -1] + + H0 = 1 / 2 * sigma_z + + return liouvillian((H0, (sigma_x, eps[1])); convention=:TDSE) +end + +@testset "GRAPE ExponentialUtilities Liouvillian" begin + # Pulse Parameters (Same as unitary) + tf = 20ns + Nsteps = 100 + tgrid = collect(range(0, tf; length=Nsteps + 1)) + + # Target gate on the qubit subspace (embed 2×2 into 3×3) + U_tgt = ComplexF64[0 1; 1 0] + + # Basis for process verification (3×3 density operators) + basis = Matrix{ComplexF64}[] + push!(basis, [1 0; 0 0]) + push!(basis, [0 0; 0 1]) + push!(basis, 1 / 2 * [1 1; 1 1]) + push!(basis, 1 / 2 * [1 im; -im 1]) + basis_tgt = [U_tgt * b * U_tgt' for b in basis] + + eps = [ + 0.2 * (1 + 0.05 * rand()) * + QuantumControl.Shapes.flattop.(tgrid, T=tf, t_rise=0.3, func=:blackman) + for _ in 1:2 + ] + + lvn = qubit_lvn(eps) + + trajectories = [ + Trajectory( + initial_state=reshape(basis[i], :), + generator=lvn, + target_state=reshape(basis_tgt[i], :) + ) for i in eachindex(basis) + ] + + problem = ControlProblem( + trajectories=trajectories, + tlist=tgrid, + iter_stop=10, + prop_method=:expv, + gradient_method=:taylor, + prop_expv_kwargs=(; ishermitian=false), + J_T=QuantumControl.Functionals.J_T_re + ) + + println("Starting optimization with QuantumPropagators: ExponentialUtilities...") + result = QuantumControl.optimize(problem; method=GRAPE, print_iters=true) + display(result) + + @test isfinite(result.J_T) +end + +@testset "GRAPE ExponentialUtilities Gradgen Error" begin + tf = 1.0 + Nsteps = 10 + tgrid = collect(range(0, tf; length=Nsteps + 1)) + + basis = Matrix{ComplexF64}[] + push!(basis, [1 0; 0 0]) + push!(basis, [0 0; 0 1]) + push!(basis, 1 / 2 * [1 1; 1 1]) + push!(basis, 1 / 2 * [1 im; -im 1]) + + U_tgt = ComplexF64[0 1; 1 0] + basis_tgt = [U_tgt * b * U_tgt' for b in basis] + + eps = [0.2 * QuantumControl.Shapes.flattop.(tgrid, T=tf, t_rise=0.3) for _ in 1:2] + lvn = qubit_lvn(eps) + + trajectories = [ + Trajectory( + initial_state=reshape(basis[i], :), + generator=lvn, + target_state=reshape(basis_tgt[i], :), + ) for i in eachindex(basis) + ] + + problem = ControlProblem( + trajectories=trajectories, + tlist=tgrid, + iter_stop=1, + prop_method=:expv, + prop_expv_kwargs=(; ishermitian=false), + grad_prop_method=:expv, + gradient_method=:gradgen, + J_T=QuantumControl.Functionals.J_T_re + ) + + res = QuantumControl.optimize(problem; method=GRAPE, print_iters=false) + @test occursin("gradient_method=:gradgen", res.message) +end diff --git a/test/test_prop_interfaces.jl b/test/test_prop_interfaces.jl index a6f9a733..68f0e858 100644 --- a/test/test_prop_interfaces.jl +++ b/test/test_prop_interfaces.jl @@ -4,6 +4,7 @@ using QuantumControlTestUtils.RandomObjects: random_state_vector, random_dynamic using QuantumPropagators: QuantumPropagators, init_prop using QuantumPropagators.Interfaces: check_propagator using StableRNGs: StableRNG +using ExponentialUtilities using OrdinaryDiffEq: OrdinaryDiffEq using QuantumPropagators.Shapes: flattop @@ -11,7 +12,7 @@ using QuantumPropagators.Shapes: flattop @testset "Cheby Propagator Interface" begin N = 10 - tlist = collect(range(0, 10, length = 101)) + tlist = collect(range(0, 10, length=101)) rng = StableRNG(677918056) Ψ = random_state_vector(N; rng) Ĥ = random_dynamic_generator(N, tlist; rng) @@ -20,10 +21,10 @@ using QuantumPropagators.Shapes: flattop Ψ, Ĥ, tlist; - method = :cheby, - backward = false, - inplace = true, - verbose = false + method=:cheby, + backward=false, + inplace=true, + verbose=false ) @test check_propagator(propagator) @@ -32,10 +33,10 @@ using QuantumPropagators.Shapes: flattop Ψ, Ĥ, tlist; - method = :cheby, - backward = false, - inplace = false, - verbose = false + method=:cheby, + backward=false, + inplace=false, + verbose=false ) @test check_propagator(propagator) @@ -44,10 +45,10 @@ using QuantumPropagators.Shapes: flattop Ψ, Ĥ, tlist; - method = :cheby, - backward = true, - inplace = true, - verbose = false + method=:cheby, + backward=true, + inplace=true, + verbose=false ) @test check_propagator(propagator) @@ -58,7 +59,7 @@ end @testset "Newton Propagator Interface" begin N = 10 - tlist = collect(range(0, 10, length = 101)) + tlist = collect(range(0, 10, length=101)) rng = StableRNG(677918057) Ψ = random_state_vector(N; rng) Ĥ = random_dynamic_generator(N, tlist; rng) @@ -67,10 +68,10 @@ end Ψ, Ĥ, tlist; - method = :newton, - backward = false, - inplace = true, - verbose = false + method=:newton, + backward=false, + inplace=true, + verbose=false ) @test check_propagator(propagator) @@ -80,10 +81,10 @@ end Ψ, Ĥ, tlist; - method = :newton, - backward = false, - inplace = false, - verbose = false + method=:newton, + backward=false, + inplace=false, + verbose=false ) end @@ -91,10 +92,10 @@ end Ψ, Ĥ, tlist; - method = :newton, - backward = true, - inplace = true, - verbose = false + method=:newton, + backward=true, + inplace=true, + verbose=false ) @test check_propagator(propagator) @@ -106,7 +107,54 @@ end N = 10 rng = StableRNG(677918057) - tlist = collect(range(0, 10, length = 101)) + tlist = collect(range(0, 10, length=101)) + Ψ = random_state_vector(N; rng) + Ĥ = random_dynamic_generator(N, tlist; rng) + + propagator = init_prop( + Ψ, + Ĥ, + tlist; + method=:expprop, + backward=false, + inplace=true, + verbose=false + ) + + @test check_propagator(propagator) + + propagator = init_prop( + Ψ, + Ĥ, + tlist; + method=:expprop, + backward=false, + inplace=false, + verbose=false + ) + + @test check_propagator(propagator) + + propagator = init_prop( + Ψ, + Ĥ, + tlist; + method=:expprop, + backward=true, + inplace=true, + verbose=false + ) + + @test check_propagator(propagator) + +end + + +@testset "ExponentialUtilities Propagator Interface" begin + + N = 10 + rng = StableRNG(677918057) + tlist = collect(range(0, 10, length=101)) Ψ = random_state_vector(N; rng) Ĥ = random_dynamic_generator(N, tlist; rng) @@ -114,10 +162,10 @@ end Ψ, Ĥ, tlist; - method = :expprop, - backward = false, - inplace = true, - verbose = false + method=ExponentialUtilities, + backward=false, + inplace=true, + verbose=false ) @test check_propagator(propagator) @@ -126,10 +174,10 @@ end Ψ, Ĥ, tlist; - method = :expprop, - backward = false, - inplace = false, - verbose = false + method=:expv, + backward=false, + inplace=false, + verbose=false ) @test check_propagator(propagator) @@ -138,10 +186,10 @@ end Ψ, Ĥ, tlist; - method = :expprop, - backward = true, - inplace = true, - verbose = false + method=ExponentialUtilities, + backward=true, + inplace=true, + verbose=false ) @test check_propagator(propagator) @@ -153,20 +201,20 @@ end N = 10 T = 10.0 - tlist = collect(range(0, T, length = 101)) + tlist = collect(range(0, T, length=101)) rng = StableRNG(677918057) Ψ = random_state_vector(N; rng) - amplitudes = [t -> flattop(t; T, t_rise = 0.3 * T),] + amplitudes = [t -> flattop(t; T, t_rise=0.3 * T),] Ĥ = random_dynamic_generator(N, tlist; amplitudes, rng) propagator = init_prop( Ψ, Ĥ, tlist; - method = OrdinaryDiffEq, - backward = false, - inplace = true, - verbose = false + method=OrdinaryDiffEq, + backward=false, + inplace=true, + verbose=false ) @test check_propagator(propagator) @@ -174,10 +222,10 @@ end Ψ, Ĥ, tlist; - method = OrdinaryDiffEq, - backward = false, - inplace = false, - verbose = false + method=OrdinaryDiffEq, + backward=false, + inplace=false, + verbose=false ) @test check_propagator(propagator) @@ -185,10 +233,10 @@ end Ψ, Ĥ, tlist; - method = OrdinaryDiffEq, - backward = true, - inplace = true, - verbose = false + method=OrdinaryDiffEq, + backward=true, + inplace=true, + verbose=false ) @test check_propagator(propagator) @@ -199,7 +247,7 @@ end N = 10 T = 10.0 - tlist = collect(range(0, T, length = 101)) + tlist = collect(range(0, T, length=101)) rng = StableRNG(677918057) Ψ = random_state_vector(N; rng) Ĥ = random_dynamic_generator(N, tlist; rng) @@ -208,11 +256,11 @@ end Ψ, Ĥ, tlist; - method = OrdinaryDiffEq, - pwc = true, - backward = false, - inplace = true, - verbose = false + method=OrdinaryDiffEq, + pwc=true, + backward=false, + inplace=true, + verbose=false ) @test check_propagator(propagator) @@ -220,11 +268,11 @@ end Ψ, Ĥ, tlist; - method = OrdinaryDiffEq, - pwc = true, - backward = false, - inplace = false, - verbose = false + method=OrdinaryDiffEq, + pwc=true, + backward=false, + inplace=false, + verbose=false ) @test check_propagator(propagator) @@ -232,11 +280,11 @@ end Ψ, Ĥ, tlist; - method = OrdinaryDiffEq, - pwc = true, - backward = true, - inplace = true, - verbose = false + method=OrdinaryDiffEq, + pwc=true, + backward=true, + inplace=true, + verbose=false ) @test check_propagator(propagator) diff --git a/test/test_propagate.jl b/test/test_propagate.jl index cb855112..b4b42dff 100644 --- a/test/test_propagate.jl +++ b/test/test_propagate.jl @@ -2,6 +2,7 @@ using Test using QuantumPropagators using QuantumPropagators: Cheby using QuantumPropagators.Storage +using ExponentialUtilities using LinearAlgebra using UnicodePlots using StaticArrays: @SMatrix, SVector, @SVector @@ -27,12 +28,14 @@ using StaticArrays: @SMatrix, SVector, @SVector @test eltype(storage) == ComplexF64 Ψ_out = propagate(Ψ0, generator, tlist; method = :expprop, storage = storage) + Ψ_out_expv = propagate(Ψ0, generator, tlist; method = ExponentialUtilities) Ψ_expected = ComplexF64[-1/√2, -1im/√2] # note the phases pop0 = abs.(storage[1, :]) .^ 2 SHOWPLOT && println(lineplot(tlist ./ π, pop0, ylim = [0, 1], title = "fw prop")) @test norm(Ψ_out - Ψ_expected) < 1e-12 + @test norm(Ψ_out_expv - Ψ_expected) < 1e-9 @test pop0[end] ≈ 0.5 # Propagating backward in time should exactly reverse the dynamics (since From 1e3e9843886a6b5f141242851b51b0f797156621 Mon Sep 17 00:00:00 2001 From: Malte Krug Date: Wed, 4 Feb 2026 03:42:29 +0100 Subject: [PATCH 02/13] docs: clarify GRAPE usage for ExponentialUtilities --- docs/src/methods.md | 50 ++++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 47 insertions(+), 3 deletions(-) diff --git a/docs/src/methods.md b/docs/src/methods.md index f5cd1fbc..a553049a 100644 --- a/docs/src/methods.md +++ b/docs/src/methods.md @@ -15,7 +15,7 @@ As discussed in the [Overview](@ref overview_approaches), time propagation can b We consider this especially in the piecewise-constant case (`pwc=true` in [`propagate`](@ref)/[`init_prop`](@ref)), which is required for the traditional optimization methods [GRAPE](https://juliaquantumcontrol.github.io/GRAPE.jl/stable/) and [Krotov](https://juliaquantumcontrol.github.io/Krotov.jl/stable/). In these propagations, the time-dependent generator ``\op{H}(t)`` is [evaluated](@ref QuantumPropagators.Controls.evaluate) to a constant operator ``\op{H}`` on each interval of the time grid. The analytical solution to the Schrödinger or Liouville equation is well known, and propagation step simply has to evaluate the application of the time evolution operator ``\op{U} = \exp[-i \op{H} dt]`` to the state ``|Ψ⟩``. The following methods are built in to `QuantumPropagators`: * [`ExpProp`](@ref method_expprop) – constructs ``\op{U}`` explicitly and then applies it to ``|Ψ⟩`` - * [`ExponentialUtilities`](@ref method_exponentialutilities) – applies ``\exp(\op{H} dt) |Ψ⟩`` using Krylov methods without explicitly forming ``\op{U}`` + * [`ExponentialUtilities`](@ref method_exponentialutilities) – applies ``\op{U} |Ψ⟩`` using Krylov methods without explicitly forming ``\op{U}`` * [`Cheby`](@ref method_cheby) — expansion of ``\op{U} |Ψ⟩`` into Chebychev polynomials, valid if ``\op{H}`` has real eigenvalues * [`Newton`](@ref method_newton) – expansion of ``\op{U} |Ψ⟩`` into Newton polynomials, valid if ``\op{H}`` has complex eigenvalues (non-Hermitian Hamiltonian, Liouvillian) @@ -78,8 +78,14 @@ using ExponentialUtilities and then passed as `method=ExponentialUtilities` (or `method=:expv`) to [`propagate`](@ref) or [`init_prop`](@ref): -```@docs -init_prop(state, generator, tlist, method::Val{:ExponentialUtilities}; kwargs...) +```julia +init_prop( + state, + generator, + tlist, + method::Val{:ExponentialUtilities}; + kwargs... +) ``` This method evaluates ``\exp(-i \op{H} dt) |Ψ⟩`` via a Krylov expv algorithm @@ -87,6 +93,44 @@ without explicitly forming the matrix exponential. It is therefore often a good fit for larger systems or matrix-free operators where direct matrix exponentiation is too costly. +Example initialization: + +```julia +using ExponentialUtilities + +expv_propagator = init_prop( + state, + generator, + tlist; + method=ExponentialUtilities, # or :expv + inplace=QuantumPropagators.Interfaces.supports_inplace(state), + backward=false, + verbose=false, + parameters=nothing, + expv_kwargs=(; ishermitian=false), # set for non-Hermitian generators + convert_state=typeof(state), + convert_operator=Matrix{ComplexF64}, +) +``` + +**Method-specific keyword arguments** + +* `expv_kwargs`: NamedTuple of keyword arguments forwarded to + `ExponentialUtilities.expv`. Use this to set `ishermitian=false` for + non-Hermitian generators (e.g., Liouvillians). +* `convert_state`: Type to which to temporarily convert the state before + calling `expv`. +* `convert_operator`: Type to which to convert the operator before calling + `expv`. + +**GRAPE note** + +* When using `method=:expv` with GRAPE, set `gradient_method=:taylor`. + The default `:gradgen` path uses `GradVector`/`GradgenOperator`, which is + not currently supported by the ExponentialUtilities Krylov backend. +* For non-Hermitian generators (e.g., Liouvillians), pass + `prop_expv_kwargs=(; ishermitian=false)` to avoid Hermitian-only paths. + **Advantages** * Avoids explicit construction of ``\op{U}`` From 2ac8b7d6bc984bf5cc1a871651cb3edf95c7e911 Mon Sep 17 00:00:00 2001 From: Malte Krug Date: Wed, 4 Feb 2026 13:21:10 +0100 Subject: [PATCH 03/13] Format code and add expv coverage --- ...antumPropagatorsExponentialUtilitiesExt.jl | 23 ++- src/newton.jl | 14 +- test/runtests.jl | 4 +- test/test_grape_exponentialutilities.jl | 26 +-- test/test_grape_exputils_liouvillian.jl | 61 +++---- test/test_prop_interfaces.jl | 164 +++++++++--------- test/test_propagate.jl | 7 + 7 files changed, 157 insertions(+), 142 deletions(-) diff --git a/ext/QuantumPropagatorsExponentialUtilitiesExt.jl b/ext/QuantumPropagatorsExponentialUtilitiesExt.jl index da9c59dc..19833c4a 100644 --- a/ext/QuantumPropagatorsExponentialUtilitiesExt.jl +++ b/ext/QuantumPropagatorsExponentialUtilitiesExt.jl @@ -29,8 +29,11 @@ LinearAlgebra.ishermitian(op::_SizedOp) = LinearAlgebra.ishermitian(op.A) LinearAlgebra.mul!(y, op::_SizedOp, x) = mul!(y, op.A, x) _ensure_size_dim(A) = - hasmethod(size, Tuple{typeof(A), Int}) ? A : - (hasmethod(size, Tuple{typeof(A)}) ? _SizedOp(A) : A) + if hasmethod(size, Tuple{typeof(A),Int}) + A + else + (hasmethod(size, Tuple{typeof(A)}) ? _SizedOp(A) : A) + end """Propagator for Krylov expv propagation via ExponentialUtilities (`method=ExponentialUtilities`). @@ -85,7 +88,7 @@ expv_propagator = init_prop( ) ``` -initializes an [`ExpvPropagator`](@ref). + initializes an `ExpvPropagator`. # Method-specific keyword arguments @@ -121,7 +124,7 @@ function init_prop( t = tlist[1] if backward n = length(tlist) - 1 - t = float(tlist[n + 1]) + t = float(tlist[n+1]) end GT = typeof(generator) OT = typeof(G) @@ -154,16 +157,18 @@ function prop_step!(propagator::ExpvPropagator) @timeit_debug propagator.timing_data "prop_step!" begin if nameof(typeof(propagator.state)) == :GradVector && nameof(parentmodule(typeof(propagator.state))) == :QuantumGradientGenerators - throw(ArgumentError( - "ExponentialUtilities propagation does not support GRAPE `gradient_method=:gradgen`. " * - "Use `gradient_method=:taylor` instead." - )) + throw( + ArgumentError( + "ExponentialUtilities propagation does not support GRAPE `gradient_method=:gradgen`. " * + "Use `gradient_method=:taylor` instead." + ) + ) end H = propagator.genop n = propagator.n tlist = getfield(propagator, :tlist) (0 < n < length(tlist)) || return nothing - dt = tlist[n + 1] - tlist[n] + dt = tlist[n+1] - tlist[n] if propagator.backward dt = -dt end diff --git a/src/newton.jl b/src/newton.jl index b5cafd97..65f755b3 100644 --- a/src/newton.jl +++ b/src/newton.jl @@ -34,7 +34,7 @@ mutable struct NewtonWrk{T} restarts::Int64 m_max::Int64 timing_data::TimerOutput - function NewtonWrk(v0::T; m_max::Int64=10, _timing_data=TimerOutput()) where {T} + function NewtonWrk(v0::T; m_max::Int64 = 10, _timing_data = TimerOutput()) where {T} if m_max <= 2 error("Newton propagation requires m_max > 2") end @@ -60,8 +60,8 @@ mutable struct NewtonWrk{T} end -lbound(array::OffsetArray, dim=1) = first(axes(array)[dim]) -ubound(array::OffsetArray, dim=1) = last(axes(array)[dim]) +lbound(array::OffsetArray, dim = 1) = first(axes(array)[dim]) +ubound(array::OffsetArray, dim = 1) = last(axes(array)[dim]) function leja_radius(z) @@ -281,9 +281,9 @@ function newton!(Ψ, H, dt, wrk; kwargs...) wrk.v, H, _dt; - extended=true, - norm_min=norm_min, - _timing_data=wrk.timing_data + extended = true, + norm_min = norm_min, + _timing_data = wrk.timing_data ) end if m == 1 && s == 0 @@ -294,7 +294,7 @@ function newton!(Ψ, H, dt, wrk; kwargs...) break end @timeit_debug wrk.timing_data "diagonalize_hessenberg_matrix" begin - ritz = diagonalize_hessenberg_matrix(Hess, m, accumulate=true) + ritz = diagonalize_hessenberg_matrix(Hess, m, accumulate = true) end # In the first iteration, the radius will be determined diff --git a/test/runtests.jl b/test/runtests.jl index f2cbd52b..4044c3cb 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -98,7 +98,9 @@ using SafeTestsets include("test_grape_exponentialutilities.jl") end - println("\n* GRAPE ExponentialUtilities Liouvillian (test_grape_exputils_liouvillian.jl):") + println( + "\n* GRAPE ExponentialUtilities Liouvillian (test_grape_exputils_liouvillian.jl):" + ) @time @safetestset "GRAPE ExponentialUtilities Liouvillian" begin include("test_grape_exputils_liouvillian.jl") end diff --git a/test/test_grape_exponentialutilities.jl b/test/test_grape_exponentialutilities.jl index e9801e0d..14e50a3c 100644 --- a/test/test_grape_exponentialutilities.jl +++ b/test/test_grape_exponentialutilities.jl @@ -9,21 +9,21 @@ using ExponentialUtilities @testset "GRAPE with ExponentialUtilities" begin rng = StableRNG(20260204) problem = dummy_control_problem( - N=4, - n_trajectories=1, - n_controls=1, - n_steps=10, - dt=0.1, - density=1.0, - hermitian=true, - complex_operators=true, - prop_method=ExponentialUtilities, - J_T=J_T_sm, - iter_stop=1, - rng=rng + N = 4, + n_trajectories = 1, + n_controls = 1, + n_steps = 10, + dt = 0.1, + density = 1.0, + hermitian = true, + complex_operators = true, + prop_method = ExponentialUtilities, + J_T = J_T_sm, + iter_stop = 1, + rng = rng ) - res = QuantumControl.optimize(problem; method=GRAPE, iter_stop=1, verbose=false) + res = QuantumControl.optimize(problem; method = GRAPE, iter_stop = 1, verbose = false) @test res.iter >= res.iter_start @test isfinite(res.J_T) diff --git a/test/test_grape_exputils_liouvillian.jl b/test/test_grape_exputils_liouvillian.jl index a3229caa..9e54edb4 100644 --- a/test/test_grape_exputils_liouvillian.jl +++ b/test/test_grape_exputils_liouvillian.jl @@ -17,14 +17,14 @@ function qubit_lvn(eps) H0 = 1 / 2 * sigma_z - return liouvillian((H0, (sigma_x, eps[1])); convention=:TDSE) + return liouvillian((H0, (sigma_x, eps[1])); convention = :TDSE) end @testset "GRAPE ExponentialUtilities Liouvillian" begin # Pulse Parameters (Same as unitary) tf = 20ns Nsteps = 100 - tgrid = collect(range(0, tf; length=Nsteps + 1)) + tgrid = collect(range(0, tf; length = Nsteps + 1)) # Target gate on the qubit subspace (embed 2×2 into 3×3) U_tgt = ComplexF64[0 1; 1 0] @@ -38,33 +38,34 @@ end basis_tgt = [U_tgt * b * U_tgt' for b in basis] eps = [ - 0.2 * (1 + 0.05 * rand()) * - QuantumControl.Shapes.flattop.(tgrid, T=tf, t_rise=0.3, func=:blackman) - for _ in 1:2 + 0.2 * + (1 + 0.05 * rand()) * + QuantumControl.Shapes.flattop.(tgrid, T = tf, t_rise = 0.3, func = :blackman) + for _ = 1:2 ] lvn = qubit_lvn(eps) trajectories = [ Trajectory( - initial_state=reshape(basis[i], :), - generator=lvn, - target_state=reshape(basis_tgt[i], :) + initial_state = reshape(basis[i], :), + generator = lvn, + target_state = reshape(basis_tgt[i], :) ) for i in eachindex(basis) ] problem = ControlProblem( - trajectories=trajectories, - tlist=tgrid, - iter_stop=10, - prop_method=:expv, - gradient_method=:taylor, - prop_expv_kwargs=(; ishermitian=false), - J_T=QuantumControl.Functionals.J_T_re + trajectories = trajectories, + tlist = tgrid, + iter_stop = 10, + prop_method = :expv, + gradient_method = :taylor, + prop_expv_kwargs = (; ishermitian = false), + J_T = QuantumControl.Functionals.J_T_re ) println("Starting optimization with QuantumPropagators: ExponentialUtilities...") - result = QuantumControl.optimize(problem; method=GRAPE, print_iters=true) + result = QuantumControl.optimize(problem; method = GRAPE, print_iters = true) display(result) @test isfinite(result.J_T) @@ -73,7 +74,7 @@ end @testset "GRAPE ExponentialUtilities Gradgen Error" begin tf = 1.0 Nsteps = 10 - tgrid = collect(range(0, tf; length=Nsteps + 1)) + tgrid = collect(range(0, tf; length = Nsteps + 1)) basis = Matrix{ComplexF64}[] push!(basis, [1 0; 0 0]) @@ -84,28 +85,28 @@ end U_tgt = ComplexF64[0 1; 1 0] basis_tgt = [U_tgt * b * U_tgt' for b in basis] - eps = [0.2 * QuantumControl.Shapes.flattop.(tgrid, T=tf, t_rise=0.3) for _ in 1:2] + eps = [0.2 * QuantumControl.Shapes.flattop.(tgrid, T = tf, t_rise = 0.3) for _ = 1:2] lvn = qubit_lvn(eps) trajectories = [ Trajectory( - initial_state=reshape(basis[i], :), - generator=lvn, - target_state=reshape(basis_tgt[i], :), + initial_state = reshape(basis[i], :), + generator = lvn, + target_state = reshape(basis_tgt[i], :), ) for i in eachindex(basis) ] problem = ControlProblem( - trajectories=trajectories, - tlist=tgrid, - iter_stop=1, - prop_method=:expv, - prop_expv_kwargs=(; ishermitian=false), - grad_prop_method=:expv, - gradient_method=:gradgen, - J_T=QuantumControl.Functionals.J_T_re + trajectories = trajectories, + tlist = tgrid, + iter_stop = 1, + prop_method = :expv, + prop_expv_kwargs = (; ishermitian = false), + grad_prop_method = :expv, + gradient_method = :gradgen, + J_T = QuantumControl.Functionals.J_T_re ) - res = QuantumControl.optimize(problem; method=GRAPE, print_iters=false) + res = QuantumControl.optimize(problem; method = GRAPE, print_iters = false) @test occursin("gradient_method=:gradgen", res.message) end diff --git a/test/test_prop_interfaces.jl b/test/test_prop_interfaces.jl index 68f0e858..00a05e31 100644 --- a/test/test_prop_interfaces.jl +++ b/test/test_prop_interfaces.jl @@ -12,7 +12,7 @@ using QuantumPropagators.Shapes: flattop @testset "Cheby Propagator Interface" begin N = 10 - tlist = collect(range(0, 10, length=101)) + tlist = collect(range(0, 10, length = 101)) rng = StableRNG(677918056) Ψ = random_state_vector(N; rng) Ĥ = random_dynamic_generator(N, tlist; rng) @@ -21,10 +21,10 @@ using QuantumPropagators.Shapes: flattop Ψ, Ĥ, tlist; - method=:cheby, - backward=false, - inplace=true, - verbose=false + method = :cheby, + backward = false, + inplace = true, + verbose = false ) @test check_propagator(propagator) @@ -33,10 +33,10 @@ using QuantumPropagators.Shapes: flattop Ψ, Ĥ, tlist; - method=:cheby, - backward=false, - inplace=false, - verbose=false + method = :cheby, + backward = false, + inplace = false, + verbose = false ) @test check_propagator(propagator) @@ -45,10 +45,10 @@ using QuantumPropagators.Shapes: flattop Ψ, Ĥ, tlist; - method=:cheby, - backward=true, - inplace=true, - verbose=false + method = :cheby, + backward = true, + inplace = true, + verbose = false ) @test check_propagator(propagator) @@ -59,7 +59,7 @@ end @testset "Newton Propagator Interface" begin N = 10 - tlist = collect(range(0, 10, length=101)) + tlist = collect(range(0, 10, length = 101)) rng = StableRNG(677918057) Ψ = random_state_vector(N; rng) Ĥ = random_dynamic_generator(N, tlist; rng) @@ -68,10 +68,10 @@ end Ψ, Ĥ, tlist; - method=:newton, - backward=false, - inplace=true, - verbose=false + method = :newton, + backward = false, + inplace = true, + verbose = false ) @test check_propagator(propagator) @@ -81,10 +81,10 @@ end Ψ, Ĥ, tlist; - method=:newton, - backward=false, - inplace=false, - verbose=false + method = :newton, + backward = false, + inplace = false, + verbose = false ) end @@ -92,10 +92,10 @@ end Ψ, Ĥ, tlist; - method=:newton, - backward=true, - inplace=true, - verbose=false + method = :newton, + backward = true, + inplace = true, + verbose = false ) @test check_propagator(propagator) @@ -107,7 +107,7 @@ end N = 10 rng = StableRNG(677918057) - tlist = collect(range(0, 10, length=101)) + tlist = collect(range(0, 10, length = 101)) Ψ = random_state_vector(N; rng) Ĥ = random_dynamic_generator(N, tlist; rng) @@ -115,10 +115,10 @@ end Ψ, Ĥ, tlist; - method=:expprop, - backward=false, - inplace=true, - verbose=false + method = :expprop, + backward = false, + inplace = true, + verbose = false ) @test check_propagator(propagator) @@ -127,10 +127,10 @@ end Ψ, Ĥ, tlist; - method=:expprop, - backward=false, - inplace=false, - verbose=false + method = :expprop, + backward = false, + inplace = false, + verbose = false ) @test check_propagator(propagator) @@ -139,10 +139,10 @@ end Ψ, Ĥ, tlist; - method=:expprop, - backward=true, - inplace=true, - verbose=false + method = :expprop, + backward = true, + inplace = true, + verbose = false ) @test check_propagator(propagator) @@ -154,7 +154,7 @@ end N = 10 rng = StableRNG(677918057) - tlist = collect(range(0, 10, length=101)) + tlist = collect(range(0, 10, length = 101)) Ψ = random_state_vector(N; rng) Ĥ = random_dynamic_generator(N, tlist; rng) @@ -162,10 +162,10 @@ end Ψ, Ĥ, tlist; - method=ExponentialUtilities, - backward=false, - inplace=true, - verbose=false + method = ExponentialUtilities, + backward = false, + inplace = true, + verbose = false ) @test check_propagator(propagator) @@ -174,10 +174,10 @@ end Ψ, Ĥ, tlist; - method=:expv, - backward=false, - inplace=false, - verbose=false + method = :expv, + backward = false, + inplace = false, + verbose = false ) @test check_propagator(propagator) @@ -186,10 +186,10 @@ end Ψ, Ĥ, tlist; - method=ExponentialUtilities, - backward=true, - inplace=true, - verbose=false + method = ExponentialUtilities, + backward = true, + inplace = true, + verbose = false ) @test check_propagator(propagator) @@ -201,20 +201,20 @@ end N = 10 T = 10.0 - tlist = collect(range(0, T, length=101)) + tlist = collect(range(0, T, length = 101)) rng = StableRNG(677918057) Ψ = random_state_vector(N; rng) - amplitudes = [t -> flattop(t; T, t_rise=0.3 * T),] + amplitudes = [t -> flattop(t; T, t_rise = 0.3 * T),] Ĥ = random_dynamic_generator(N, tlist; amplitudes, rng) propagator = init_prop( Ψ, Ĥ, tlist; - method=OrdinaryDiffEq, - backward=false, - inplace=true, - verbose=false + method = OrdinaryDiffEq, + backward = false, + inplace = true, + verbose = false ) @test check_propagator(propagator) @@ -222,10 +222,10 @@ end Ψ, Ĥ, tlist; - method=OrdinaryDiffEq, - backward=false, - inplace=false, - verbose=false + method = OrdinaryDiffEq, + backward = false, + inplace = false, + verbose = false ) @test check_propagator(propagator) @@ -233,10 +233,10 @@ end Ψ, Ĥ, tlist; - method=OrdinaryDiffEq, - backward=true, - inplace=true, - verbose=false + method = OrdinaryDiffEq, + backward = true, + inplace = true, + verbose = false ) @test check_propagator(propagator) @@ -247,7 +247,7 @@ end N = 10 T = 10.0 - tlist = collect(range(0, T, length=101)) + tlist = collect(range(0, T, length = 101)) rng = StableRNG(677918057) Ψ = random_state_vector(N; rng) Ĥ = random_dynamic_generator(N, tlist; rng) @@ -256,11 +256,11 @@ end Ψ, Ĥ, tlist; - method=OrdinaryDiffEq, - pwc=true, - backward=false, - inplace=true, - verbose=false + method = OrdinaryDiffEq, + pwc = true, + backward = false, + inplace = true, + verbose = false ) @test check_propagator(propagator) @@ -268,11 +268,11 @@ end Ψ, Ĥ, tlist; - method=OrdinaryDiffEq, - pwc=true, - backward=false, - inplace=false, - verbose=false + method = OrdinaryDiffEq, + pwc = true, + backward = false, + inplace = false, + verbose = false ) @test check_propagator(propagator) @@ -280,11 +280,11 @@ end Ψ, Ĥ, tlist; - method=OrdinaryDiffEq, - pwc=true, - backward=true, - inplace=true, - verbose=false + method = OrdinaryDiffEq, + pwc = true, + backward = true, + inplace = true, + verbose = false ) @test check_propagator(propagator) diff --git a/test/test_propagate.jl b/test/test_propagate.jl index b4b42dff..f4a4bc42 100644 --- a/test/test_propagate.jl +++ b/test/test_propagate.jl @@ -29,6 +29,11 @@ using StaticArrays: @SMatrix, SVector, @SVector Ψ_out = propagate(Ψ0, generator, tlist; method = :expprop, storage = storage) Ψ_out_expv = propagate(Ψ0, generator, tlist; method = ExponentialUtilities) + Ψ_out_expv_ni = + propagate(Ψ0, generator, tlist; method = ExponentialUtilities, inplace = false) + generator_herm = (Hermitian(Ĥ),) + Ψ_out_expv_herm = + propagate(Ψ0, generator_herm, tlist; method = ExponentialUtilities, inplace = true) Ψ_expected = ComplexF64[-1/√2, -1im/√2] # note the phases pop0 = abs.(storage[1, :]) .^ 2 @@ -36,6 +41,8 @@ using StaticArrays: @SMatrix, SVector, @SVector @test norm(Ψ_out - Ψ_expected) < 1e-12 @test norm(Ψ_out_expv - Ψ_expected) < 1e-9 + @test norm(Ψ_out_expv_ni - Ψ_expected) < 1e-9 + @test norm(Ψ_out_expv_herm - Ψ_expected) < 1e-9 @test pop0[end] ≈ 0.5 # Propagating backward in time should exactly reverse the dynamics (since From d855586a46c683e8d95b42e424d5b3abb907ae04 Mon Sep 17 00:00:00 2001 From: Michael Goerz Date: Mon, 16 Feb 2026 21:26:41 -0500 Subject: [PATCH 04/13] Remove unnecessary `_ensure_size_dim` and `_SizedOp` --- ...antumPropagatorsExponentialUtilitiesExt.jl | 19 ------------------- 1 file changed, 19 deletions(-) diff --git a/ext/QuantumPropagatorsExponentialUtilitiesExt.jl b/ext/QuantumPropagatorsExponentialUtilitiesExt.jl index 19833c4a..619c56c6 100644 --- a/ext/QuantumPropagatorsExponentialUtilitiesExt.jl +++ b/ext/QuantumPropagatorsExponentialUtilitiesExt.jl @@ -18,23 +18,6 @@ using QuantumPropagators.Interfaces: supports_inplace using QuantumPropagators.Generators: ScaledOperator import QuantumPropagators: init_prop, prop_step!, set_t! -struct _SizedOp{A} - A::A -end - -Base.size(op::_SizedOp) = size(op.A) -Base.size(op::_SizedOp, dim::Integer) = size(op.A)[dim] -Base.eltype(op::_SizedOp) = eltype(op.A) -LinearAlgebra.ishermitian(op::_SizedOp) = LinearAlgebra.ishermitian(op.A) -LinearAlgebra.mul!(y, op::_SizedOp, x) = mul!(y, op.A, x) - -_ensure_size_dim(A) = - if hasmethod(size, Tuple{typeof(A),Int}) - A - else - (hasmethod(size, Tuple{typeof(A)}) ? _SizedOp(A) : A) - end - """Propagator for Krylov expv propagation via ExponentialUtilities (`method=ExponentialUtilities`). @@ -182,7 +165,6 @@ function prop_step!(propagator::ExpvPropagator) else H = convert(propagator.convert_operator, _pwc_get_genop(propagator, n)) end - H = _ensure_size_dim(H) H = ScaledOperator(-1im, H) @timeit_debug propagator.timing_data "expv" begin Ψ = ExponentialUtilities.expv(dt_expv, H, Ψ; propagator.expv_kwargs...) @@ -190,7 +172,6 @@ function prop_step!(propagator::ExpvPropagator) copyto!(propagator.state, convert(typeof(propagator.state), Ψ)) else H = convert(propagator.convert_operator, _pwc_get_genop(propagator, n)) - H = _ensure_size_dim(H) H = ScaledOperator(-1im, H) @timeit_debug propagator.timing_data "expv" begin Ψ = ExponentialUtilities.expv(dt_expv, H, Ψ; propagator.expv_kwargs...) From 0ff527ba59e51c309b5160f8435a9e1d80a95d4b Mon Sep 17 00:00:00 2001 From: Michael Goerz Date: Tue, 17 Feb 2026 00:29:52 -0500 Subject: [PATCH 05/13] Revise ExponentialUtilities propagator MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Revised PR according to review: * Deleted `ext/QuantumPropagatorsQuantumGradientGeneratorsExt.jl` (type-piracy) * Removed convert_state/convert_operator fields and all conversion machinery * Removed :gradgen guard that rejected GradVector states * Removed :expv symbol alias * Removed ScaledOperator wrapping — instead passes -1im * dt as the time argument to ExponentialUtilities.expv. This avoids the use of `ScaledOperator` * Removed both "GRAPE note" from docs * Add ExponentialUtilities to InterLinks., use @extref link to ExponentialUtilities.expv * Add `test/test_exputils.jl` — new component-level test file replacing GRAPE integration tests --- Project.toml | 3 - docs/make.jl | 1 + docs/src/methods.md | 66 +------- ...antumPropagatorsExponentialUtilitiesExt.jl | 63 ++------ ...PropagatorsQuantumGradientGeneratorsExt.jl | 11 -- src/generators.jl | 1 + test/runtests.jl | 13 +- test/test_exputils.jl | 146 ++++++++++++++++++ test/test_grape_exponentialutilities.jl | 8 +- test/test_prop_interfaces.jl | 2 +- 10 files changed, 183 insertions(+), 131 deletions(-) delete mode 100644 ext/QuantumPropagatorsQuantumGradientGeneratorsExt.jl create mode 100644 test/test_exputils.jl diff --git a/Project.toml b/Project.toml index cb7666ff..bf3f866b 100644 --- a/Project.toml +++ b/Project.toml @@ -18,14 +18,12 @@ TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" [weakdeps] ExponentialUtilities = "d4d017d3-3776-5f7e-afef-a10c40355c18" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" -QuantumGradientGenerators = "a563f35e-61db-434d-8c01-8b9e3ccdfd85" RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" [extensions] QuantumPropagatorsExponentialUtilitiesExt = "ExponentialUtilities" QuantumPropagatorsODEExt = "OrdinaryDiffEq" -QuantumPropagatorsQuantumGradientGeneratorsExt = "QuantumGradientGenerators" QuantumPropagatorsRecursiveArrayToolsExt = "RecursiveArrayTools" QuantumPropagatorsStaticArraysExt = "StaticArrays" @@ -35,7 +33,6 @@ ExponentialUtilities = "1.11" OffsetArrays = "1" OrdinaryDiffEq = "6.59" ProgressMeter = "1" -QuantumGradientGenerators = "0.1" Random = "1" RecursiveArrayTools = "3.12" SpecialFunctions = "1, 2" diff --git a/docs/make.jl b/docs/make.jl index 9830c8cb..9d809358 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -47,6 +47,7 @@ links = InterLinks( "ComponentArrays" => "https://sciml.github.io/ComponentArrays.jl/stable/", "RecursiveArrayTools" => "https://docs.sciml.ai/RecursiveArrayTools/stable/", "ArrayInterface" => "https://docs.sciml.ai/ArrayInterface/stable/", + "ExponentialUtilities" => "https://docs.sciml.ai/ExponentialUtilities/stable/", "qutip" => "https://qutip.readthedocs.io/en/qutip-5.0.x/", ) diff --git a/docs/src/methods.md b/docs/src/methods.md index a553049a..c2e50b01 100644 --- a/docs/src/methods.md +++ b/docs/src/methods.md @@ -75,62 +75,19 @@ package is loaded using ExponentialUtilities ``` -and then passed as `method=ExponentialUtilities` (or `method=:expv`) to +and then passed as `method=ExponentialUtilities` to [`propagate`](@ref) or [`init_prop`](@ref): -```julia -init_prop( - state, - generator, - tlist, - method::Val{:ExponentialUtilities}; - kwargs... -) +```@docs +init_prop(state, generator, tlist, method::Val{:ExponentialUtilities}; kwargs...) ``` -This method evaluates ``\exp(-i \op{H} dt) |Ψ⟩`` via a Krylov expv algorithm -without explicitly forming the matrix exponential. It is therefore often a -good fit for larger systems or matrix-free operators where direct matrix +This method evaluates ``\exp(-i \op{H} dt) |Ψ⟩`` via a Krylov +[`expv`](@extref ExponentialUtilities :jl:function:`ExponentialUtilities.expv`) +algorithm without explicitly forming the matrix exponential. It is therefore +often a good fit for larger systems or matrix-free operators where direct matrix exponentiation is too costly. -Example initialization: - -```julia -using ExponentialUtilities - -expv_propagator = init_prop( - state, - generator, - tlist; - method=ExponentialUtilities, # or :expv - inplace=QuantumPropagators.Interfaces.supports_inplace(state), - backward=false, - verbose=false, - parameters=nothing, - expv_kwargs=(; ishermitian=false), # set for non-Hermitian generators - convert_state=typeof(state), - convert_operator=Matrix{ComplexF64}, -) -``` - -**Method-specific keyword arguments** - -* `expv_kwargs`: NamedTuple of keyword arguments forwarded to - `ExponentialUtilities.expv`. Use this to set `ishermitian=false` for - non-Hermitian generators (e.g., Liouvillians). -* `convert_state`: Type to which to temporarily convert the state before - calling `expv`. -* `convert_operator`: Type to which to convert the operator before calling - `expv`. - -**GRAPE note** - -* When using `method=:expv` with GRAPE, set `gradient_method=:taylor`. - The default `:gradgen` path uses `GradVector`/`GradgenOperator`, which is - not currently supported by the ExponentialUtilities Krylov backend. -* For non-Hermitian generators (e.g., Liouvillians), pass - `prop_expv_kwargs=(; ishermitian=false)` to avoid Hermitian-only paths. - **Advantages** * Avoids explicit construction of ``\op{U}`` @@ -145,15 +102,6 @@ expv_propagator = init_prop( **When to use** * Large, sparse, or matrix-free generators -* Piecewise-constant GRAPE-style workflows that need expv - -**GRAPE note** - -* When using `method=:expv` with GRAPE, set `gradient_method=:taylor`. - The default `:gradgen` path uses `GradVector`/`GradgenOperator`, which is - not currently supported by the ExponentialUtilities Krylov backend. -* For non-Hermitian generators (e.g., Liouvillians), pass - `prop_expv_kwargs=(; ishermitian=false)` to avoid Hermitian-only paths. ## [`Cheby`](@id method_cheby) diff --git a/ext/QuantumPropagatorsExponentialUtilitiesExt.jl b/ext/QuantumPropagatorsExponentialUtilitiesExt.jl index 619c56c6..eac9bdf5 100644 --- a/ext/QuantumPropagatorsExponentialUtilitiesExt.jl +++ b/ext/QuantumPropagatorsExponentialUtilitiesExt.jl @@ -15,7 +15,6 @@ using QuantumPropagators: _pwc_process_parameters using QuantumPropagators.Controls: get_controls using QuantumPropagators.Interfaces: supports_inplace -using QuantumPropagators.Generators: ScaledOperator import QuantumPropagators: init_prop, prop_step!, set_t! @@ -34,8 +33,6 @@ mutable struct ExpvPropagator{GT,OT,ST} <: PWCPropagator genop::OT backward::Bool inplace::Bool - convert_state::Type - convert_operator::Type expv_kwargs::NamedTuple const timing_data::TimerOutput end @@ -44,13 +41,6 @@ end set_t!(propagator::ExpvPropagator, t) = _pwc_set_t!(propagator, t) -_expv_convert_state(state) = typeof(state) -_expv_convert_operator(::Any) = Any -_expv_convert_operator(::QuantumPropagators.Generator) = Matrix{ComplexF64} -_expv_convert_operator(::QuantumPropagators.Operator) = Matrix{ComplexF64} -_expv_convert_operator(::QuantumPropagators.ScaledOperator) = Matrix{ComplexF64} - - """ ```julia using ExponentialUtilities @@ -65,8 +55,6 @@ expv_propagator = init_prop( verbose=false, parameters=nothing, expv_kwargs=(;), - convert_state=_expv_convert_state(state), - convert_operator=_expv_convert_operator(generator), _... ) ``` @@ -77,10 +65,6 @@ expv_propagator = init_prop( * `expv_kwargs`: NamedTuple of keyword arguments forwarded to `ExponentialUtilities.expv`. -* `convert_state`: Type to which to temporarily convert the state before - calling `expv`. -* `convert_operator`: Type to which to convert the operator before calling - `expv`. """ function init_prop( state, @@ -92,8 +76,6 @@ function init_prop( verbose = false, parameters = nothing, expv_kwargs = (;), - convert_state = _expv_convert_state(state), - convert_operator = _expv_convert_operator(generator), _... ) tlist = convert(Vector{Float64}, tlist) @@ -123,31 +105,14 @@ function init_prop( G, backward, inplace, - convert_state, - convert_operator, expv_kwargs, timing_data, ) end -# Aliases -init_prop(state, generator, tlist, method::Val{:expv}; kwargs...) = - init_prop(state, generator, tlist, Val(:ExponentialUtilities); kwargs...) - - function prop_step!(propagator::ExpvPropagator) @timeit_debug propagator.timing_data "prop_step!" begin - if nameof(typeof(propagator.state)) == :GradVector && - nameof(parentmodule(typeof(propagator.state))) == :QuantumGradientGenerators - throw( - ArgumentError( - "ExponentialUtilities propagation does not support GRAPE `gradient_method=:gradgen`. " * - "Use `gradient_method=:taylor` instead." - ) - ) - end - H = propagator.genop n = propagator.n tlist = getfield(propagator, :tlist) (0 < n < length(tlist)) || return nothing @@ -155,28 +120,34 @@ function prop_step!(propagator::ExpvPropagator) if propagator.backward dt = -dt end - dt_expv = complex(dt) - Ψ = convert(propagator.convert_state, propagator.state) if propagator.inplace if supports_inplace(propagator.genop) _pwc_set_genop!(propagator, n) - H = convert(propagator.convert_operator, propagator.genop) + H = propagator.genop else - H = convert(propagator.convert_operator, _pwc_get_genop(propagator, n)) + H = _pwc_get_genop(propagator, n) end - H = ScaledOperator(-1im, H) @timeit_debug propagator.timing_data "expv" begin - Ψ = ExponentialUtilities.expv(dt_expv, H, Ψ; propagator.expv_kwargs...) + Ψ = ExponentialUtilities.expv( + -1im * dt, + H, + propagator.state; + propagator.expv_kwargs... + ) end - copyto!(propagator.state, convert(typeof(propagator.state), Ψ)) + copyto!(propagator.state, Ψ) else - H = convert(propagator.convert_operator, _pwc_get_genop(propagator, n)) - H = ScaledOperator(-1im, H) + H = _pwc_get_genop(propagator, n) @timeit_debug propagator.timing_data "expv" begin - Ψ = ExponentialUtilities.expv(dt_expv, H, Ψ; propagator.expv_kwargs...) + Ψ = ExponentialUtilities.expv( + -1im * dt, + H, + propagator.state; + propagator.expv_kwargs... + ) end - setfield!(propagator, :state, convert(typeof(propagator.state), Ψ)) + setfield!(propagator, :state, Ψ) end _pwc_advance_time!(propagator) diff --git a/ext/QuantumPropagatorsQuantumGradientGeneratorsExt.jl b/ext/QuantumPropagatorsQuantumGradientGeneratorsExt.jl deleted file mode 100644 index 484d5b32..00000000 --- a/ext/QuantumPropagatorsQuantumGradientGeneratorsExt.jl +++ /dev/null @@ -1,11 +0,0 @@ -module QuantumPropagatorsQuantumGradientGeneratorsExt - -using QuantumGradientGenerators: GradgenOperator, GradVector - -Base.size(G::GradgenOperator, dim::Integer) = size(G.G, dim) -Base.similar(::GradVector, ::Type{T}, dims::Tuple{Int,Int}) where {T} = - Matrix{T}(undef, dims...) -Base.similar(::GradVector, ::Type{T}, dims::Tuple{Int}) where {T} = - Vector{T}(undef, dims[1]) - -end diff --git a/src/generators.jl b/src/generators.jl index 7f96e579..7bf61eb1 100644 --- a/src/generators.jl +++ b/src/generators.jl @@ -176,6 +176,7 @@ function Base.copyto!(tgt::Operator, src::Operator) copyto!(tgt.coeffs, src.coeffs) end +# All ops in an Operator must have the same size (they are summed together) Base.size(O::Operator) = size(O.ops[1]) Base.size(O::Operator, dim::Integer) = size(O.ops[1], dim) Base.eltype(::Type{Operator{OT,CT}}) where {OT,CT} = promote_type(eltype(OT), CT) diff --git a/test/runtests.jl b/test/runtests.jl index 4044c3cb..7d7d6aed 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -93,16 +93,9 @@ using SafeTestsets include("test_propagate_sequence.jl") end - println("\n* GRAPE with ExponentialUtilities (test_grape_exponentialutilities.jl):") - @time @safetestset "GRAPE ExponentialUtilities" begin - include("test_grape_exponentialutilities.jl") - end - - println( - "\n* GRAPE ExponentialUtilities Liouvillian (test_grape_exputils_liouvillian.jl):" - ) - @time @safetestset "GRAPE ExponentialUtilities Liouvillian" begin - include("test_grape_exputils_liouvillian.jl") + println("\n* ExponentialUtilities propagation (test_exputils.jl):") + @time @safetestset "ExponentialUtilities" begin + include("test_exputils.jl") end println("\n* Time-dependent observables (test_timedependent_observables.jl):") diff --git a/test/test_exputils.jl b/test/test_exputils.jl new file mode 100644 index 00000000..c2f88540 --- /dev/null +++ b/test/test_exputils.jl @@ -0,0 +1,146 @@ +using Test +using LinearAlgebra +using SparseArrays +using QuantumPropagators +using QuantumPropagators.Generators: liouvillian +using ExponentialUtilities +using StableRNGs: StableRNG + + +@testset "Forward propagation (Hermitian)" begin + # Simple TLS Rabi cycling: compare with exact matrix exponential + Ψ₀ = ComplexF64[1, 0] + Ĥ = ComplexF64[0 0.5; 0.5 0] + tlist = collect(range(0, 1.5π, length = 101)) + generator = (Ĥ,) + + Ψ_out = propagate(Ψ₀, generator, tlist; method = ExponentialUtilities) + Ψ_expected = ComplexF64[-1/√2, -1im/√2] + + @test norm(Ψ_out - Ψ_expected) < 1e-9 +end + + +@testset "Forward propagation (not in-place)" begin + Ψ₀ = ComplexF64[1, 0] + Ĥ = ComplexF64[0 0.5; 0.5 0] + tlist = collect(range(0, 1.5π, length = 101)) + generator = (Ĥ,) + + Ψ_out = propagate(Ψ₀, generator, tlist; method = ExponentialUtilities, inplace = false) + Ψ_expected = ComplexF64[-1/√2, -1im/√2] + + @test norm(Ψ_out - Ψ_expected) < 1e-9 +end + + +@testset "Backward propagation" begin + Ψ₀ = ComplexF64[1, 0] + Ĥ = ComplexF64[0 0.5; 0.5 0] + tlist = collect(range(0, 1.5π, length = 101)) + generator = (Ĥ,) + + Ψ_fw = propagate(Ψ₀, generator, tlist; method = ExponentialUtilities) + Ψ_bw = propagate(Ψ_fw, generator, tlist; method = ExponentialUtilities, backward = true) + + @test norm(Ψ_bw - Ψ₀) < 1e-9 +end + + +@testset "Liouvillian propagation" begin + rng = StableRNG(20260216) + # 2-level system with decay (Liouvillian) + σ_x = ComplexF64[0 1; 1 0] + σ_z = ComplexF64[1 0; 0 -1] + γ = 0.1 # decay rate + c_op = √γ * ComplexF64[0 1; 0 0] # lowering operator + + H₀ = 0.5 * σ_z + ℒ = liouvillian(H₀, [c_op]; convention = :TDSE) + + # Start from |0⟩⟨0| + ρ₀ = ComplexF64[1 0; 0 0] + ρ⃗₀ = reshape(ρ₀, :) + + tlist = collect(range(0, 10.0, length = 201)) + + ρ⃗_out = propagate( + ρ⃗₀, + (ℒ,), + tlist; + method = ExponentialUtilities, + expv_kwargs = (; ishermitian = false) + ) + + ρ_out = reshape(ρ⃗_out, 2, 2) + + # Density matrix must have trace 1 and be positive semidefinite + @test tr(ρ_out) ≈ 1.0 atol = 1e-9 + @test all(eigvals(Hermitian(ρ_out)) .>= -1e-10) + + # Compare with exact solution via matrix exponential + L_full = Array(ℒ) + U = exp(-1im * L_full * tlist[end]) + ρ⃗_expected = U * ρ⃗₀ + @test norm(ρ⃗_out - ρ⃗_expected) < 1e-9 +end + + +@testset "Liouvillian backward propagation" begin + σ_z = ComplexF64[1 0; 0 -1] + γ = 0.01 + c_op = √γ * ComplexF64[0 1; 0 0] + + H₀ = 0.5 * σ_z + ℒ = liouvillian(H₀, [c_op]; convention = :TDSE) + + ρ₀ = ComplexF64[1 0; 0 0] + ρ⃗₀ = reshape(ρ₀, :) + + tlist = collect(range(0, 1.0, length = 51)) + + ρ⃗_fw = propagate( + ρ⃗₀, + (ℒ,), + tlist; + method = ExponentialUtilities, + expv_kwargs = (; ishermitian = false) + ) + + ρ⃗_bw = propagate( + ρ⃗_fw, + (ℒ,), + tlist; + method = ExponentialUtilities, + backward = true, + expv_kwargs = (; ishermitian = false) + ) + + # For non-unitary dynamics, backward propagation won't exactly recover ρ₀, + # but applying exp(-iLdt) backward should still be consistent with exp(+iLdt) + # Compare with exact matrix exponential backward propagation + L_full = Array(ℒ) + U_bw = exp(+1im * L_full * tlist[end]) + ρ⃗_bw_expected = U_bw * ρ⃗_fw + @test norm(ρ⃗_bw - ρ⃗_bw_expected) < 1e-9 +end + + +@testset "Time-dependent generator" begin + rng = StableRNG(20260216) + Ψ₀ = ComplexF64[1, 0] + H₀ = ComplexF64[0 0; 0 1] + H₁ = ComplexF64[0 1; 1 0] + + tlist = collect(range(0, 10.0, length = 201)) + ϵ = 0.2 * ones(length(tlist)) + + generator = hamiltonian(H₀, (H₁, ϵ)) + + Ψ_out_expv = propagate(Ψ₀, generator, tlist; method = ExponentialUtilities) + Ψ_out_exp = propagate(Ψ₀, generator, tlist; method = :expprop) + + @test norm(Ψ_out_expv - Ψ_out_exp) < 1e-9 +end + +# TODO: test a GradGen propagation diff --git a/test/test_grape_exponentialutilities.jl b/test/test_grape_exponentialutilities.jl index 14e50a3c..5e13d105 100644 --- a/test/test_grape_exponentialutilities.jl +++ b/test/test_grape_exponentialutilities.jl @@ -23,7 +23,13 @@ using ExponentialUtilities rng = rng ) - res = QuantumControl.optimize(problem; method = GRAPE, iter_stop = 1, verbose = false) + res = QuantumControl.optimize( + problem; + method = GRAPE, + iter_stop = 10, + verbose = false, + rethrow_exceptions = true + ) @test res.iter >= res.iter_start @test isfinite(res.J_T) diff --git a/test/test_prop_interfaces.jl b/test/test_prop_interfaces.jl index 00a05e31..4c232dbd 100644 --- a/test/test_prop_interfaces.jl +++ b/test/test_prop_interfaces.jl @@ -174,7 +174,7 @@ end Ψ, Ĥ, tlist; - method = :expv, + method = ExponentialUtilities, backward = false, inplace = false, verbose = false From a369d7b3287c6f235078eb214425b6d56b5ee3f9 Mon Sep 17 00:00:00 2001 From: Michael Goerz Date: Thu, 26 Feb 2026 16:08:17 -0500 Subject: [PATCH 06/13] Add thorough testing for ExponentialUtilities.jl This is to explore the functionalities of the package and the guarantee that everythink works reliably. --- CLAUDE.md | 11 + src/generators.jl | 2 +- test/runtests.jl | 5 + test/test_exponential_utilities.jl | 620 +++++++++++++++++++++++++++++ 4 files changed, 637 insertions(+), 1 deletion(-) create mode 100644 test/test_exponential_utilities.jl diff --git a/CLAUDE.md b/CLAUDE.md index 60b0815a..c9a25276 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -76,11 +76,20 @@ Uses SafeTestsets for isolated test execution. Tests are comprehensive and inclu Running `make test` prints out coverage data in a table. +Tips for writing tests: + +- Use functions like `random_state_vector`, `random_dynamic_generator`, `random_matrix` from `QuantumControlTestUtils.RandomObjects` to generate random objects. These should always be given an explicit `StableRNG` as `rng`, with a unique seed + +- When a new seed is required, obtain one with `julia --project=test -e 'using QuantumControlTestUtils.RandomObjects: randseed; print(randseed())'` + + If necessary, detailed line-by-line coverage information can be obtained by running julia --project=test -e 'include("devrepl.jl"); generate_coverage_html()' after `make test`. This will produce html files inside the `coverage` subfolder, with `coverage/src` mirroring the structure of the `src` folder of `.jl` files. Lines with `` are not covered. Ignore the raw tracefiles in the `.coverage` subfolder. ### Development Environment + The project uses a sophisticated development setup: + - Development REPL (devrepl.jl) with Revise.jl for hot reloading - Automatic dependency management via installorg.jl script - Integrated documentation building and serving @@ -100,4 +109,6 @@ Each Julia function that is not explicitly private or has a name starting with a * When adding a new dependency to any `Project.toml` file, run `make distclean`, and then `make test/Manifest.toml`, `make docs/Manifest.toml`, etc. to recreate manifest files as necessary. +* In order to get the documentation (docstring) of any function, e.g., `some_func`, in any package in the test environment, e.g., `SomeLib`, run `julia --project=. -e 'using REPL; using SomeLib; print(Base.doc(SomeLib.some_func))'` + * Never commit any changes or ask to commit. I will always create git commits manually. diff --git a/src/generators.jl b/src/generators.jl index 7bf61eb1..e2e70cfc 100644 --- a/src/generators.jl +++ b/src/generators.jl @@ -256,7 +256,7 @@ function Base.show(io::IO, O::ScaledOperator{CT,OT}) where {CT,OT} print(io, "ScaledOperator($(O.coeff), $(O.operator))") end -function Base.show(io::IO, ::MIME"text/plain", O::ScaledOperator{CT,OT}) where {CT,OT} +function Base.show(io::IO, ::MIME"text/plain", O::ScaledOperator{CT,<:Operator}) where {CT} Base.summary(io, O) println(io, "\n operator.ops::$(typeof(O.operator.ops)):") for op in O.operator.ops diff --git a/test/runtests.jl b/test/runtests.jl index 7d7d6aed..7b8bc702 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -78,6 +78,11 @@ using SafeTestsets include("test_expprop.jl") end + println("\n* ExponentialUtilities (test_exponential_utilities.jl):") + @time @safetestset "ExponentialUtilities" begin + include("test_exponential_utilities.jl") + end + println("\n* Propagator Interfaces (test_prop_interfaces.jl):") @time @safetestset "Propagator Interfaces" begin include("test_prop_interfaces.jl") diff --git a/test/test_exponential_utilities.jl b/test/test_exponential_utilities.jl new file mode 100644 index 00000000..ae1bd1c3 --- /dev/null +++ b/test/test_exponential_utilities.jl @@ -0,0 +1,620 @@ +using Test + +using QuantumPropagators.Controls: evaluate +using QuantumPropagators.Generators: ScaledOperator +using StableRNGs: StableRNG +using QuantumControlTestUtils.RandomObjects: + random_state_vector, random_dynamic_generator, random_matrix +using LinearAlgebra: norm, ishermitian +import ExponentialUtilities +import StaticArrays: SMatrix, SVector +import OffsetArrays: OffsetVector +using BenchmarkTools: @btimed + +RUN_BENCHMARKS = false + + +@testset "expv with dense complex non-Hermitian matrix" begin + + N = 10 + rng = StableRNG(677968056) + A = random_matrix(N; rng, complex = true, hermitian = false) + Ψ₀ = random_state_vector(N; rng) + dt = 0.1 + t = -1im * dt + Ψ = ExponentialUtilities.expv(t, A, Ψ₀) + Ψ_expected = exp(t * A) * Ψ₀ + @test norm(Ψ - Ψ_expected) < 1e-12 + +end + + +@testset "expv with dense real non-Hermitian matrix" begin + + N = 10 + rng = StableRNG(2945201349) + A = random_matrix(N; rng, complex = false, hermitian = false) + Ψ₀ = random_state_vector(N; rng) + dt = 0.1 + t = -1im * dt + Ψ = ExponentialUtilities.expv(t, A, Ψ₀) + Ψ_expected = exp(t * A) * Ψ₀ + @test norm(Ψ - Ψ_expected) < 1e-12 + +end + + +@testset "expv with dense complex Hermitian matrix and complex/real time" begin + + # Hermitian A triggers Lanczos instead of Arnoldi + N = 10 + rng = StableRNG(697691674) + A = random_matrix(N; rng, hermitian = true) + @test ishermitian(A) + Ψ₀ = random_state_vector(N; rng) + dt = 0.1 + t = -1im * dt + Ψ = ExponentialUtilities.expv(t, A, Ψ₀) + Ψ_expected = exp(t * A) * Ψ₀ + @test norm(Ψ - Ψ_expected) < 1e-12 + if RUN_BENCHMARKS + b1 = @btimed ExponentialUtilities.expv($t, $A, $Ψ₀) + end + + # We can also check against the code path where we explicitly disable + # `ishermitian` detection + Ψ = ExponentialUtilities.expv(t, A, Ψ₀; ishermitian = false) + @test norm(Ψ - Ψ_expected) < 1e-12 + if RUN_BENCHMARKS + b2 = @btimed ExponentialUtilities.expv($t, $A, $Ψ₀; ishermitian = false) + @test b2.time > b1.time + end + + # Instead of using complex time, we can also use real time and absorb the + # factor `-1im` in `A` + O = ScaledOperator(-1im, A) + @test !ishermitian(O) + Ψ = ExponentialUtilities.expv(dt, O, Ψ₀) + @test norm(Ψ - Ψ_expected) < 1e-12 + if RUN_BENCHMARKS + b3 = @btimed ExponentialUtilities.expv($dt, $O, $Ψ₀) + @test b3.time > b1.time + end + +end + + +@testset "expv with dense real symmetric matrix" begin + + # Symmetric (Hermitian) `A` triggers Lanczos instead of Arnoldi + N = 10 + rng = StableRNG(1056344208) + A = random_matrix(N; rng, complex = false, hermitian = true) + @test ishermitian(A) + Ψ₀ = random_state_vector(N; rng) + dt = 0.1 + t = -1im * dt + Ψ = ExponentialUtilities.expv(t, A, Ψ₀) + Ψ_expected = exp(t * A) * Ψ₀ + @test norm(Ψ - Ψ_expected) < 1e-12 + if RUN_BENCHMARKS + b1 = @btimed ExponentialUtilities.expv($t, $A, $Ψ₀) + end + # We can also check against the code path where we explicitly disable + # `ishermitian` detection + Ψ = ExponentialUtilities.expv(t, A, Ψ₀; ishermitian = false) + @test norm(Ψ - Ψ_expected) < 1e-12 + if RUN_BENCHMARKS + b2 = @btimed ExponentialUtilities.expv($t, $A, $Ψ₀; ishermitian = false) + @test b2.time > b1.time + end + +end + + +@testset "expv with sparse complex non-Hermitian matrix" begin + + N = 50 + rng = StableRNG(1788065498) + A = random_matrix(N; rng, density = 0.3, hermitian = false) + @test !ishermitian(A) + Ψ₀ = random_state_vector(N; rng) + dt = 0.1 + t = -1im * dt + Ψ = ExponentialUtilities.expv(t, A, Ψ₀) + Ψ_expected = exp(t * Matrix(A)) * Ψ₀ + @test norm(Ψ - Ψ_expected) < 1e-12 + +end + + +@testset "expv with sparse complex Hermitian matrix" begin + + # Sparse Hermitian A triggers Lanczos + N = 50 + rng = StableRNG(2648750134) + A = random_matrix(N; rng, density = 0.3, hermitian = true) + @test ishermitian(A) + Ψ₀ = random_state_vector(N; rng) + dt = 0.1 + t = -1im * dt + Ψ = ExponentialUtilities.expv(t, A, Ψ₀) + Ψ_expected = exp(t * Matrix(A)) * Ψ₀ + @test norm(Ψ - Ψ_expected) < 1e-12 + +end + + +@testset "expv with Operator" begin + + N = 10 + rng = StableRNG(677968056) + tlist = collect(range(0, 10, length = 101)) + H = random_dynamic_generator(N, tlist; rng) + Ψ₀ = random_state_vector(N; rng) + A = evaluate(H, tlist, 1) + dt = tlist[2] - tlist[1] + t = -1im * dt + Ψ = ExponentialUtilities.expv(t, A, Ψ₀) + Ψ_expected = exp(t * Array(A)) * Ψ₀ + @test norm(Ψ - Ψ_expected) < 1e-12 + +end + + +@testset "expv with StaticArrays SMatrix and SVector" begin + + N = 4 + rng = StableRNG(814242294) + A = SMatrix{N,N}(random_matrix(N; rng, hermitian = true)) + Ψ₀ = SVector{N}(random_state_vector(N; rng)) + dt = 0.1 + t = -1im * dt + Ψ = ExponentialUtilities.expv(t, A, Ψ₀) + Ψ_expected = exp(t * Matrix(A)) * Vector(Ψ₀) + @test Ψ isa SVector{N} + @test norm(Ψ - Ψ_expected) < 1e-12 + +end + + +@testset "expv with pre-built KrylovSubspace" begin + + # arnoldi builds the Krylov subspace once; expv(t, Ks) reuses it, + # which is efficient when applying the same operator at many time steps + N = 10 + rng = StableRNG(3750746111) + A = random_matrix(N; rng, hermitian = true) + Ψ₀ = random_state_vector(N; rng) + dt = 0.1 + t = -1im * dt + Ks = ExponentialUtilities.arnoldi(A, Ψ₀) + Ψ = ExponentialUtilities.expv(t, Ks) + Ψ_expected = exp(t * A) * Ψ₀ + @test norm(Ψ - Ψ_expected) < 1e-12 + +end + + +@testset "expv! in-place with pre-built KrylovSubspace" begin + + N = 10 + rng = StableRNG(1187394246) + A = random_matrix(N; rng, hermitian = true) + Ψ₀ = random_state_vector(N; rng) + dt = 0.1 + t = -1im * dt + Ks = ExponentialUtilities.arnoldi(A, Ψ₀) + Ψ = similar(Ψ₀, promote_type(typeof(t), eltype(A), eltype(Ψ₀))) + ExponentialUtilities.expv!(Ψ, t, Ks) + Ψ_expected = exp(t * A) * Ψ₀ + @test norm(Ψ - Ψ_expected) < 1e-12 + +end + + +@testset "expv! in-place with dense complex Hermitian matrix and complex/real time" begin + + N = 10 + rng = StableRNG(3503257612) + A = random_matrix(N; rng, hermitian = true) + @test ishermitian(A) + Ψ₀ = random_state_vector(N; rng) + dt = 0.1 + t = -1im * dt + Ψ_expected = exp(t * A) * Ψ₀ + + # Baseline: expv + Ψ = ExponentialUtilities.expv(t, A, Ψ₀) + @test norm(Ψ - Ψ_expected) < 1e-12 + if RUN_BENCHMARKS + b0 = @btimed ExponentialUtilities.expv($t, $A, $Ψ₀) + end + + # Lanczos path: arnoldi detects ishermitian(A) and builds a symmetric + # tridiagonal H; expv! then uses the fast SymTridiagonal eigen path. + # For Hermitian A the H matrix is real: ExpvCache element type = real(eltype(A)). + Ks1 = ExponentialUtilities.arnoldi(A, Ψ₀) + w = similar(Ψ₀) + ExponentialUtilities.expv!(w, t, Ks1) + @test norm(w - Ψ_expected) < 1e-12 + if RUN_BENCHMARKS + cache1 = ExponentialUtilities.ExpvCache{real(eltype(A))}(Ks1.m) + b1 = @btimed ExponentialUtilities.expv!($w, $t, $Ks1; cache = $cache1) + @test b1.time < b0.time + end + + # Arnoldi path: force full Arnoldi; H is upper Hessenberg, so expv! falls + # back to the generic Higham matrix exponential, which is slower. + # Arnoldi (including when forced on a Hermitian A) builds complex H: + # ExpvCache element type = eltype(A). + Ks2 = ExponentialUtilities.arnoldi(A, Ψ₀; ishermitian = false) + ExponentialUtilities.expv!(w, t, Ks2) + @test norm(w - Ψ_expected) < 1e-12 + if RUN_BENCHMARKS + cache2 = ExponentialUtilities.ExpvCache{eltype(A)}(Ks2.m) + b2 = @btimed ExponentialUtilities.expv!($w, $t, $Ks2; cache = $cache2) + # The timings here are a lot flakier, so we don't + # @test b2.time > b1.time + end + + # Real-time path: absorb `-1im` into the operator via ScaledOperator; + # the scaled operator is not Hermitian so arnoldi uses Arnoldi iteration, + # giving an upper Hessenberg H (ComplexF64), same as Ks2. + O = ScaledOperator(-1im, A) + @test !ishermitian(O) + Ks3 = ExponentialUtilities.arnoldi(O, Ψ₀) + ExponentialUtilities.expv!(w, dt, Ks3) + @test norm(w - Ψ_expected) < 1e-12 + if RUN_BENCHMARKS + cache3 = ExponentialUtilities.ExpvCache{eltype(O)}(Ks3.m) + b3 = @btimed ExponentialUtilities.expv!($w, $dt, $Ks3; cache = $cache3) + # The timings here are a lot flakier, so we don't + # @test b3.time > b1.time + end + +end + + +@testset "expv propagation loop over time steps" begin + + # Propagate Ψ through the first few time steps of a random dynamic + # generator, reusing a single KrylovSubspace. + # At each step: arnoldi! updates Ks in-place for the current (A, Ψ), + # then expv(t, Ks) returns the new Ψ without mutating the input. + N = 10 + rng = StableRNG(1373376837) + tlist = collect(range(0, 10, length = 101)) + H = random_dynamic_generator(N, tlist; rng) + Ψ₀ = random_state_vector(N; rng) + dt = tlist[2] - tlist[1] + t = -1im * dt + n_steps = 5 + + # Reference: sequential dense matrix exponentials + Ψ_expected = copy(Ψ₀) + for n = 1:n_steps + A = evaluate(H, tlist, n) + Ψ_expected = exp(t * Array(A)) * Ψ_expected + end + + # Loop with arnoldi! + expv(t, Ks), reusing Ks across all steps + A₀ = evaluate(H, tlist, 1) + Ks = ExponentialUtilities.arnoldi(A₀, Ψ₀) + Ψ = copy(Ψ₀) + for n = 1:n_steps + A = evaluate(H, tlist, n) + ExponentialUtilities.arnoldi!(Ks, A, Ψ) + Ψ = ExponentialUtilities.expv(t, Ks) + end + @test norm(Ψ - Ψ_expected) < 1e-12 + +end + + +@testset "expv! propagation loop over time steps" begin + + # Propagate Ψ through the first few time steps of a random dynamic + # generator, reusing a single KrylovSubspace and ExpvCache. + # At each step n: arnoldi! updates Ks in-place for the current (A, Ψ), + # then expv! overwrites Ψ. arnoldi! reads b before expv! writes w, so + # passing the same vector for both is safe. + # Ks and cache must remain type-consistent across the loop: the H element + # type is fixed at construction (real for Lanczos, complex for Arnoldi), + # so all arnoldi! calls must use the same ishermitian setting as the + # initial arnoldi call. + N = 10 + rng = StableRNG(2266238742) + tlist = collect(range(0, 10, length = 101)) + H = random_dynamic_generator(N, tlist; rng) + Ψ₀ = random_state_vector(N; rng) + dt = tlist[2] - tlist[1] + t = -1im * dt + n_steps = 5 + + # Reference: sequential dense matrix exponentials + Ψ_expected = copy(Ψ₀) + for n = 1:n_steps + A = evaluate(H, tlist, n) + Ψ_expected = exp(t * Array(A)) * Ψ_expected + end + + # Loop with arnoldi! + expv!, reusing Ks and cache across all steps. + # The H element type is real(T) for Lanczos and T for Arnoldi, where + # T = promote_type(eltype(A), eltype(Ψ)) is the Krylov basis vector type. + # Hermitian A₀ → Lanczos → real H → ExpvCache{real(T)} + A₀ = evaluate(H, tlist, 1) + @test ishermitian(A₀) + Ks = ExponentialUtilities.arnoldi(A₀, Ψ₀) + T = promote_type(eltype(A₀), eltype(Ψ₀)) + cache = ExponentialUtilities.ExpvCache{real(T)}(Ks.m) + Ψ = copy(Ψ₀) + if RUN_BENCHMARKS + # One might think that these in-place versions are completely + # non-allocating. Unfortunately, that's not the case. The allocations + # don't disappear if this is wrapped into a function. + b1 = @btimed ExponentialUtilities.arnoldi!($Ks, $A₀, $Ψ) + @test_broken b1.alloc == 0 + b2 = @btimed ExponentialUtilities.expv!($Ψ, $t, $Ks; cache = $cache) + @test_broken b2.alloc == 0 + end + Ψ = copy(Ψ₀) + for n = 1:n_steps + A = evaluate(H, tlist, n) + ExponentialUtilities.arnoldi!(Ks, A, Ψ) + ExponentialUtilities.expv!(Ψ, t, Ks; cache = cache) + end + @test norm(Ψ - Ψ_expected) < 1e-12 + +end + + +@testset "expv! propagation loop over time steps with non-Hermitian generator" begin + + # Non-Hermitian generators require Arnoldi iteration (not Lanczos), giving + # an upper Hessenberg H matrix (ComplexF64). The ExpvCache element type + # must be ComplexF64 accordingly. Otherwise the pattern is the same as for + # Hermitian generators. + N = 10 + rng = StableRNG(1573710086) + tlist = collect(range(0, 10, length = 101)) + H = random_dynamic_generator(N, tlist; rng, hermitian = false) + Ψ₀ = random_state_vector(N; rng) + dt = tlist[2] - tlist[1] + t = -1im * dt + n_steps = 5 + + # Reference: sequential dense matrix exponentials + Ψ_expected = copy(Ψ₀) + for n = 1:n_steps + A = evaluate(H, tlist, n) + Ψ_expected = exp(t * Array(A)) * Ψ_expected + end + + # Loop with arnoldi! + expv!, reusing Ks and cache across all steps. + # The H element type is T for Arnoldi, where + # T = promote_type(eltype(A), eltype(Ψ)) is the Krylov basis vector type. + # Non-Hermitian A₀ → Arnoldi → H element type T → ExpvCache{T} + A₀ = evaluate(H, tlist, 1) + @test !ishermitian(A₀) + Ks = ExponentialUtilities.arnoldi(A₀, Ψ₀) + T = promote_type(eltype(A₀), eltype(Ψ₀)) + cache = ExponentialUtilities.ExpvCache{T}(Ks.m) + Ψ = copy(Ψ₀) + for n = 1:n_steps + A = evaluate(H, tlist, n) + @test !ishermitian(A) + ExponentialUtilities.arnoldi!(Ks, A, Ψ) + ExponentialUtilities.expv!(Ψ, t, Ks; cache = cache) + end + @test norm(Ψ - Ψ_expected) < 1e-12 + +end + + +@testset "expv! propagation loop using error-estimate interface" begin + + # Variation of the propagation loop using expv!(w, t, A, b, Ks, cache): + # this interface builds the Krylov subspace internally via Lanczos and + # terminates adaptively when Saad's error estimate drops below the + # requested tolerance (default atol=1e-8, rtol=1e-4). + # Only Hermitian operators are supported: get_subspace_cache requires a + # Lanczos (real H) subspace, which arnoldi builds automatically when + # ishermitian(A) is true. + # Passing w == b (same Ψ vector) is safe: b is fully consumed into + # V[:,1] before w is written at the end of the call. + N = 10 + rng = StableRNG(106614078) + tlist = collect(range(0, 10, length = 101)) + H = random_dynamic_generator(N, tlist; rng) + Ψ₀ = random_state_vector(N; rng) + dt = tlist[2] - tlist[1] + t = -1im * dt + n_steps = 5 + + # Reference: sequential dense matrix exponentials + Ψ_expected = copy(Ψ₀) + for n = 1:n_steps + A = evaluate(H, tlist, n) + Ψ_expected = exp(t * Array(A)) * Ψ_expected + end + + A₀ = evaluate(H, tlist, 1) + @test ishermitian(A₀) + Ks = ExponentialUtilities.arnoldi(A₀, Ψ₀) # auto-detects Hermitian → Lanczos + cache = ExponentialUtilities.get_subspace_cache(Ks) + Ψ = copy(Ψ₀) + for n = 1:n_steps + A = evaluate(H, tlist, n) + @test ishermitian(A) + ExponentialUtilities.expv!(Ψ, t, A, Ψ, Ks, cache) + end + # accuracy is governed by the per-step rtol=1e-4 default + @test norm(Ψ - Ψ_expected) < 1e-3 + +end + + +@testset "expv! propagation loop using error-estimate interface with tight tolerance" begin + + # Same error-estimate interface, but with atol=rtol=1e-14. The termination + # condition ε = atol + rtol*‖b‖ is now so tight that the algorithm uses + # all m = min(Ks.maxiter, N) Lanczos steps, spanning the full N-dimensional + # Krylov space and giving near-machine-precision results. + N = 10 + rng = StableRNG(1653744503) + tlist = collect(range(0, 10, length = 101)) + H = random_dynamic_generator(N, tlist; rng) + Ψ₀ = random_state_vector(N; rng) + dt = tlist[2] - tlist[1] + t = -1im * dt + n_steps = 5 + + # Reference: sequential dense matrix exponentials + Ψ_expected = copy(Ψ₀) + for n = 1:n_steps + A = evaluate(H, tlist, n) + Ψ_expected = exp(t * Array(A)) * Ψ_expected + end + + A₀ = evaluate(H, tlist, 1) + Ks = ExponentialUtilities.arnoldi(A₀, Ψ₀) # auto-detects Hermitian → Lanczos + cache = ExponentialUtilities.get_subspace_cache(Ks) + Ψ = copy(Ψ₀) + for n = 1:n_steps + A = evaluate(H, tlist, n) + ExponentialUtilities.expv!(Ψ, t, A, Ψ, Ks, cache; atol = 1e-14, rtol = 1e-14) + end + @test norm(Ψ - Ψ_expected) < 1e-12 + +end + + +@testset "expv! propagation loop using error-estimate interface with non-Hermitian generator" begin + + # The error-estimate interface expv!(w, t, A, b, Ks, cache) only supports + # Hermitian operators. For non-Hermitian generators, get_subspace_cache + # throws an ErrorException because HermitianSubspaceCache requires a + # Lanczos subspace with real H. Additionally, the method signature + # requires KrylovSubspace{T,B,B} (second and third type params equal), + # which only holds for Lanczos (Float64, Float64) not Arnoldi + # (ComplexF64, Float64). + N = 10 + rng = StableRNG(208995058) + tlist = collect(range(0, 10, length = 101)) + H = random_dynamic_generator(N, tlist; rng, hermitian = false) + Ψ₀ = random_state_vector(N; rng) + dt = tlist[2] - tlist[1] + t = -1im * dt + n_steps = 5 + + # Reference: sequential dense matrix exponentials + Ψ_expected = copy(Ψ₀) + for n = 1:n_steps + A = evaluate(H, tlist, n) + Ψ_expected = exp(t * Array(A)) * Ψ_expected + end + + # get_subspace_cache throws ErrorException for non-Hermitian Ks; + # a MethodError on expv! dispatch would also occur since the method + # requires KrylovSubspace{T,B,B} but non-Hermitian Ks has B=ComplexF64 ≠ Float64 + A₀ = evaluate(H, tlist, 1) + @test !ishermitian(A₀) + Ks = ExponentialUtilities.arnoldi(A₀, Ψ₀) + @test_broken begin + cache = ExponentialUtilities.get_subspace_cache(Ks) + Ψ = copy(Ψ₀) + for n = 1:n_steps + A = evaluate(H, tlist, n) + ExponentialUtilities.expv!(Ψ, t, A, Ψ, Ks, cache) + end + norm(Ψ - Ψ_expected) < 1e-3 + end + +end + + +@testset "expv with mode=:error_estimate" begin + + # mode=:error_estimate uses adaptive Krylov iteration with an on-the-fly + # error estimate; default rtol ≈ 3e-4 (= sqrt(tol) with tol=1e-7) + N = 10 + rng = StableRNG(3050745481) + A = random_matrix(N; rng, hermitian = true) + Ψ₀ = random_state_vector(N; rng) + dt = 0.1 + t = -1im * dt + Ψ = ExponentialUtilities.expv(t, A, Ψ₀; mode = :error_estimate) + Ψ_expected = exp(t * A) * Ψ₀ + @test norm(Ψ - Ψ_expected) < 1e-4 + +end + + +@testset "expv with custom Krylov subspace size m" begin + + N = 20 + rng = StableRNG(1375059820) + A = random_matrix(N; rng, hermitian = true) + Ψ₀ = random_state_vector(N; rng) + dt = 0.1 + t = -1im * dt + Ψ = ExponentialUtilities.expv(t, A, Ψ₀; m = 8) + Ψ_expected = exp(t * A) * Ψ₀ + @test norm(Ψ - Ψ_expected) < 1e-12 + +end + + +@testset "expv with iop (incomplete orthogonalization)" begin + + # iop=k uses only the last k vectors for re-orthogonalization, + # reducing cost for large operators at the price of some accuracy + N = 20 + rng = StableRNG(1918914874) + A = random_matrix(N; rng, hermitian = false) + Ψ₀ = random_state_vector(N; rng) + dt = 0.1 + t = -1im * dt + Ψ = ExponentialUtilities.expv(t, A, Ψ₀; iop = 2) + Ψ_expected = exp(t * A) * Ψ₀ + @test norm(Ψ - Ψ_expected) < 1e-10 + +end + + +@testset "expv with custom opnorm" begin + + # opnorm can be provided to avoid recomputing it internally, + # useful when the norm is expensive or already known + N = 10 + rng = StableRNG(1355318934) + A = random_matrix(N; rng, hermitian = false) + Ψ₀ = random_state_vector(N; rng) + dt = 0.1 + t = -1im * dt + Ψ = ExponentialUtilities.expv(t, A, Ψ₀; opnorm = norm(A, Inf)) + Ψ_expected = exp(t * A) * Ψ₀ + @test norm(Ψ - Ψ_expected) < 1e-12 + +end + + +@testset "expv with OffsetArray state vector (broken)" begin + + # expv does not support OffsetVector as state vector b: arnoldi builds + # V = similar(b, T, (n, m+1)) which yields an OffsetMatrix, and the + # subsequent broadcast in expv! throws a DimensionMismatch + N = 9 + rng = StableRNG(2841259413) + A = random_matrix(N; rng, hermitian = true) + Ψ₀ = OffsetVector(random_state_vector(N; rng), (-(N÷2)):(N÷2)) + dt = 0.1 + t = -1im * dt + @test_broken begin + Ψ = ExponentialUtilities.expv(t, A, Ψ₀) + Ψ_expected = exp(t * A) * collect(Ψ₀) + norm(Ψ - Ψ_expected) < 1e-12 + end + +end From 1e6c8c3c6d7e283966201f293205cc55d4d6bdd9 Mon Sep 17 00:00:00 2001 From: Malte Krug Date: Sat, 28 Feb 2026 01:30:23 +0100 Subject: [PATCH 07/13] Optimize in-place using expv! with Krylov subspace caching --- ...antumPropagatorsExponentialUtilitiesExt.jl | 59 ++++++++++++++++--- 1 file changed, 50 insertions(+), 9 deletions(-) diff --git a/ext/QuantumPropagatorsExponentialUtilitiesExt.jl b/ext/QuantumPropagatorsExponentialUtilitiesExt.jl index eac9bdf5..bd112e5a 100644 --- a/ext/QuantumPropagatorsExponentialUtilitiesExt.jl +++ b/ext/QuantumPropagatorsExponentialUtilitiesExt.jl @@ -22,7 +22,7 @@ import QuantumPropagators: init_prop, prop_step!, set_t! This is a [`PWCPropagator`](@ref). """ -mutable struct ExpvPropagator{GT,OT,ST} <: PWCPropagator +mutable struct ExpvPropagator{GT,OT,ST,KST,CT} <: PWCPropagator const generator::GT state::ST t::Float64 # time at which current `state` is defined @@ -31,6 +31,8 @@ mutable struct ExpvPropagator{GT,OT,ST} <: PWCPropagator parameters::AbstractDict controls genop::OT + Ks::KST + cache::CT backward::Bool inplace::Bool expv_kwargs::NamedTuple @@ -91,10 +93,29 @@ function init_prop( n = length(tlist) - 1 t = float(tlist[n+1]) end + + if inplace + A₀ = QuantumPropagators.Controls.evaluate(generator, tlist, n) + Ks = ExponentialUtilities.arnoldi(A₀, state; expv_kwargs...) + if haskey(expv_kwargs, :mode) && expv_kwargs[:mode] == :error_estimate + cache = ExponentialUtilities.get_subspace_cache(Ks) + else + T = promote_type(eltype(A₀), eltype(state)) + cache_type = ishermitian(A₀) ? real(T) : T + cache = ExponentialUtilities.ExpvCache{cache_type}(Ks.m) + end + else + Ks = nothing + cache = nothing + end + GT = typeof(generator) OT = typeof(G) ST = typeof(state) - return ExpvPropagator{GT,OT,ST}( + KST = typeof(Ks) + CT = typeof(cache) + + return ExpvPropagator{GT,OT,ST,KST,CT}( generator, inplace ? copy(state) : state, t, @@ -103,6 +124,8 @@ function init_prop( parameters, controls, G, + Ks, + cache, backward, inplace, expv_kwargs, @@ -129,14 +152,32 @@ function prop_step!(propagator::ExpvPropagator) H = _pwc_get_genop(propagator, n) end @timeit_debug propagator.timing_data "expv" begin - Ψ = ExponentialUtilities.expv( - -1im * dt, - H, - propagator.state; - propagator.expv_kwargs... - ) + if haskey(propagator.expv_kwargs, :mode) && + propagator.expv_kwargs[:mode] == :error_estimate + ExponentialUtilities.expv!( + propagator.state, + -1im * dt, + H, + propagator.state, + propagator.Ks, + propagator.cache; + propagator.expv_kwargs..., + ) + else + ExponentialUtilities.arnoldi!( + propagator.Ks, + H, + propagator.state; + propagator.expv_kwargs..., + ) + ExponentialUtilities.expv!( + propagator.state, + -1im * dt, + propagator.Ks; + cache = propagator.cache, + ) + end end - copyto!(propagator.state, Ψ) else H = _pwc_get_genop(propagator, n) @timeit_debug propagator.timing_data "expv" begin From cc2bdd149db7896fba421780d0f53d52fbf08039 Mon Sep 17 00:00:00 2001 From: Malte Krug Date: Sat, 28 Feb 2026 01:38:01 +0100 Subject: [PATCH 08/13] Removed GRAPE specific tests. Will moved to GRAPE.jl later for testing ExponentialUtilites --- test/test_grape_exponentialutilities.jl | 36 -------- test/test_grape_exputils_liouvillian.jl | 112 ------------------------ 2 files changed, 148 deletions(-) delete mode 100644 test/test_grape_exponentialutilities.jl delete mode 100644 test/test_grape_exputils_liouvillian.jl diff --git a/test/test_grape_exponentialutilities.jl b/test/test_grape_exponentialutilities.jl deleted file mode 100644 index 5e13d105..00000000 --- a/test/test_grape_exponentialutilities.jl +++ /dev/null @@ -1,36 +0,0 @@ -using Test -using GRAPE -using QuantumControl -using QuantumControl.Functionals: J_T_sm -using QuantumControlTestUtils.DummyOptimization: dummy_control_problem -using StableRNGs: StableRNG -using ExponentialUtilities - -@testset "GRAPE with ExponentialUtilities" begin - rng = StableRNG(20260204) - problem = dummy_control_problem( - N = 4, - n_trajectories = 1, - n_controls = 1, - n_steps = 10, - dt = 0.1, - density = 1.0, - hermitian = true, - complex_operators = true, - prop_method = ExponentialUtilities, - J_T = J_T_sm, - iter_stop = 1, - rng = rng - ) - - res = QuantumControl.optimize( - problem; - method = GRAPE, - iter_stop = 10, - verbose = false, - rethrow_exceptions = true - ) - - @test res.iter >= res.iter_start - @test isfinite(res.J_T) -end diff --git a/test/test_grape_exputils_liouvillian.jl b/test/test_grape_exputils_liouvillian.jl deleted file mode 100644 index 9e54edb4..00000000 --- a/test/test_grape_exputils_liouvillian.jl +++ /dev/null @@ -1,112 +0,0 @@ -using QuantumControl -using GRAPE -using LinearAlgebra -using QuantumPropagators.Controls: get_controls -using ExponentialUtilities -using Test - -# Scale -const GHz = 1.0 -const MHz = 0.001 * GHz -const ns = 1.0 - -function qubit_lvn(eps) - # Basis vectors - sigma_x = [0 1; 1 0] - sigma_z = [1 0; 0 -1] - - H0 = 1 / 2 * sigma_z - - return liouvillian((H0, (sigma_x, eps[1])); convention = :TDSE) -end - -@testset "GRAPE ExponentialUtilities Liouvillian" begin - # Pulse Parameters (Same as unitary) - tf = 20ns - Nsteps = 100 - tgrid = collect(range(0, tf; length = Nsteps + 1)) - - # Target gate on the qubit subspace (embed 2×2 into 3×3) - U_tgt = ComplexF64[0 1; 1 0] - - # Basis for process verification (3×3 density operators) - basis = Matrix{ComplexF64}[] - push!(basis, [1 0; 0 0]) - push!(basis, [0 0; 0 1]) - push!(basis, 1 / 2 * [1 1; 1 1]) - push!(basis, 1 / 2 * [1 im; -im 1]) - basis_tgt = [U_tgt * b * U_tgt' for b in basis] - - eps = [ - 0.2 * - (1 + 0.05 * rand()) * - QuantumControl.Shapes.flattop.(tgrid, T = tf, t_rise = 0.3, func = :blackman) - for _ = 1:2 - ] - - lvn = qubit_lvn(eps) - - trajectories = [ - Trajectory( - initial_state = reshape(basis[i], :), - generator = lvn, - target_state = reshape(basis_tgt[i], :) - ) for i in eachindex(basis) - ] - - problem = ControlProblem( - trajectories = trajectories, - tlist = tgrid, - iter_stop = 10, - prop_method = :expv, - gradient_method = :taylor, - prop_expv_kwargs = (; ishermitian = false), - J_T = QuantumControl.Functionals.J_T_re - ) - - println("Starting optimization with QuantumPropagators: ExponentialUtilities...") - result = QuantumControl.optimize(problem; method = GRAPE, print_iters = true) - display(result) - - @test isfinite(result.J_T) -end - -@testset "GRAPE ExponentialUtilities Gradgen Error" begin - tf = 1.0 - Nsteps = 10 - tgrid = collect(range(0, tf; length = Nsteps + 1)) - - basis = Matrix{ComplexF64}[] - push!(basis, [1 0; 0 0]) - push!(basis, [0 0; 0 1]) - push!(basis, 1 / 2 * [1 1; 1 1]) - push!(basis, 1 / 2 * [1 im; -im 1]) - - U_tgt = ComplexF64[0 1; 1 0] - basis_tgt = [U_tgt * b * U_tgt' for b in basis] - - eps = [0.2 * QuantumControl.Shapes.flattop.(tgrid, T = tf, t_rise = 0.3) for _ = 1:2] - lvn = qubit_lvn(eps) - - trajectories = [ - Trajectory( - initial_state = reshape(basis[i], :), - generator = lvn, - target_state = reshape(basis_tgt[i], :), - ) for i in eachindex(basis) - ] - - problem = ControlProblem( - trajectories = trajectories, - tlist = tgrid, - iter_stop = 1, - prop_method = :expv, - prop_expv_kwargs = (; ishermitian = false), - grad_prop_method = :expv, - gradient_method = :gradgen, - J_T = QuantumControl.Functionals.J_T_re - ) - - res = QuantumControl.optimize(problem; method = GRAPE, print_iters = false) - @test occursin("gradient_method=:gradgen", res.message) -end From 26a34d1e20aa0f2c3265e65d12b15ba4a76f9253 Mon Sep 17 00:00:00 2001 From: Malte Krug Date: Sat, 28 Feb 2026 03:15:24 +0100 Subject: [PATCH 09/13] Fix ExpvCache type inference to respect ishermitian kwarg If A0 was accidentally hermitian, it caused an InexactError in the arnoldi method. --- ext/QuantumPropagatorsExponentialUtilitiesExt.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/ext/QuantumPropagatorsExponentialUtilitiesExt.jl b/ext/QuantumPropagatorsExponentialUtilitiesExt.jl index bd112e5a..b945e3b8 100644 --- a/ext/QuantumPropagatorsExponentialUtilitiesExt.jl +++ b/ext/QuantumPropagatorsExponentialUtilitiesExt.jl @@ -101,7 +101,8 @@ function init_prop( cache = ExponentialUtilities.get_subspace_cache(Ks) else T = promote_type(eltype(A₀), eltype(state)) - cache_type = ishermitian(A₀) ? real(T) : T + is_herm = get(expv_kwargs, :ishermitian, ishermitian(A₀)) + cache_type = is_herm ? real(T) : T cache = ExponentialUtilities.ExpvCache{cache_type}(Ks.m) end else From 48431a4adf3c79d8a6b8b37a43bf60c61c717b2f Mon Sep 17 00:00:00 2001 From: Malte Krug Date: Sun, 1 Mar 2026 15:02:13 +0100 Subject: [PATCH 10/13] Move ExponentialUtilitiesPropagator into core module Define ExponentialUtilitiesPropagator and set_t! in src so the type is publicly accessible from QuantumPropagators. Keep the extension limited to ExponentialUtilities dependent methods. --- ...antumPropagatorsExponentialUtilitiesExt.jl | 36 +++---------------- src/QuantumPropagators.jl | 1 + src/exponential_utilities_propagator.jl | 27 ++++++++++++++ test/test_prop_interfaces.jl | 3 +- 4 files changed, 35 insertions(+), 32 deletions(-) create mode 100644 src/exponential_utilities_propagator.jl diff --git a/ext/QuantumPropagatorsExponentialUtilitiesExt.jl b/ext/QuantumPropagatorsExponentialUtilitiesExt.jl index b945e3b8..c19032a5 100644 --- a/ext/QuantumPropagatorsExponentialUtilitiesExt.jl +++ b/ext/QuantumPropagatorsExponentialUtilitiesExt.jl @@ -6,41 +6,15 @@ using ExponentialUtilities using QuantumPropagators: QuantumPropagators, - PWCPropagator, + ExponentialUtilitiesPropagator, _pwc_get_max_genop, _pwc_get_genop, _pwc_set_genop!, - _pwc_set_t!, _pwc_advance_time!, _pwc_process_parameters using QuantumPropagators.Controls: get_controls using QuantumPropagators.Interfaces: supports_inplace -import QuantumPropagators: init_prop, prop_step!, set_t! - - -"""Propagator for Krylov expv propagation via ExponentialUtilities (`method=ExponentialUtilities`). - -This is a [`PWCPropagator`](@ref). -""" -mutable struct ExpvPropagator{GT,OT,ST,KST,CT} <: PWCPropagator - const generator::GT - state::ST - t::Float64 # time at which current `state` is defined - n::Int64 # index of next interval to propagate - const tlist::Vector{Float64} - parameters::AbstractDict - controls - genop::OT - Ks::KST - cache::CT - backward::Bool - inplace::Bool - expv_kwargs::NamedTuple - const timing_data::TimerOutput -end - - -set_t!(propagator::ExpvPropagator, t) = _pwc_set_t!(propagator, t) +import QuantumPropagators: init_prop, prop_step! """ @@ -61,7 +35,7 @@ expv_propagator = init_prop( ) ``` - initializes an `ExpvPropagator`. + initializes an [`ExponentialUtilitiesPropagator`](@ref). # Method-specific keyword arguments @@ -116,7 +90,7 @@ function init_prop( KST = typeof(Ks) CT = typeof(cache) - return ExpvPropagator{GT,OT,ST,KST,CT}( + return ExponentialUtilitiesPropagator{GT,OT,ST,KST,CT}( generator, inplace ? copy(state) : state, t, @@ -135,7 +109,7 @@ function init_prop( end -function prop_step!(propagator::ExpvPropagator) +function prop_step!(propagator::ExponentialUtilitiesPropagator) @timeit_debug propagator.timing_data "prop_step!" begin n = propagator.n tlist = getfield(propagator, :tlist) diff --git a/src/QuantumPropagators.jl b/src/QuantumPropagators.jl index 4d5573ac..67a9f565 100644 --- a/src/QuantumPropagators.jl +++ b/src/QuantumPropagators.jl @@ -58,6 +58,7 @@ end #! format: on include("pwc_utils.jl") +include("exponential_utilities_propagator.jl") include("cheby_propagator.jl") include("newton_propagator.jl") include("exp_propagator.jl") diff --git a/src/exponential_utilities_propagator.jl b/src/exponential_utilities_propagator.jl new file mode 100644 index 00000000..038c3229 --- /dev/null +++ b/src/exponential_utilities_propagator.jl @@ -0,0 +1,27 @@ +using TimerOutputs: TimerOutput + + +"""Propagator for Krylov `expv` propagation (`method=ExponentialUtilities`). + +This is a [`PWCPropagator`](@ref). Methods that depend on +`ExponentialUtilities` are loaded via Julia package extensions. +""" +mutable struct ExponentialUtilitiesPropagator{GT,OT,ST,KST,CT} <: PWCPropagator + const generator::GT + state::ST + t::Float64 # time at which current `state` is defined + n::Int64 # index of next interval to propagate + const tlist::Vector{Float64} + parameters::AbstractDict + controls + genop::OT + Ks::KST + cache::CT + backward::Bool + inplace::Bool + expv_kwargs::NamedTuple + const timing_data::TimerOutput +end + + +set_t!(propagator::ExponentialUtilitiesPropagator, t) = _pwc_set_t!(propagator, t) diff --git a/test/test_prop_interfaces.jl b/test/test_prop_interfaces.jl index 4c232dbd..bc36191f 100644 --- a/test/test_prop_interfaces.jl +++ b/test/test_prop_interfaces.jl @@ -1,7 +1,7 @@ using Test using QuantumControlTestUtils.RandomObjects: random_state_vector, random_dynamic_generator -using QuantumPropagators: QuantumPropagators, init_prop +using QuantumPropagators: QuantumPropagators, ExponentialUtilitiesPropagator, init_prop using QuantumPropagators.Interfaces: check_propagator using StableRNGs: StableRNG using ExponentialUtilities @@ -168,6 +168,7 @@ end verbose = false ) + @test propagator isa ExponentialUtilitiesPropagator @test check_propagator(propagator) propagator = init_prop( From ba552ac8bd53ff15f87b4c506dac58f981832ebe Mon Sep 17 00:00:00 2001 From: Malte Krug Date: Fri, 20 Mar 2026 21:25:21 +0100 Subject: [PATCH 11/13] Address PR review feedback for ExponentialUtilities propagator - Document expv_kwargs parameters (m, ishermitian, tol, opnorm, iop, mode) with accurate forwarding description (arnoldi!/expv!, not just expv) - Expand methods.md with Arnoldi/Lanczos details, AbstractArray interface requirement, and supports_vector/matrix_interface trait references - Disambiguate test set names to avoid ambiguous failure reports - Convert expv return type before setfield! in non-inplace path to handle element type promotion (e.g. real state + complex time) - Add test verifying expv type promotion behavior - Add early ArgumentError for mode=:error_estimate with non-Hermitian generators instead of deferring to obscure ExponentialUtilities error - Fix struct indentation to match project conventions (4-space) Co-Authored-By: Claude Opus 4.6 --- docs/src/methods.md | 20 ++++++-- ...antumPropagatorsExponentialUtilitiesExt.jl | 48 +++++++++++++++++-- src/exponential_utilities_propagator.jl | 28 +++++------ test/runtests.jl | 8 ++-- test/test_exponential_utilities.jl | 31 ++++++++++++ 5 files changed, 108 insertions(+), 27 deletions(-) diff --git a/docs/src/methods.md b/docs/src/methods.md index c2e50b01..a875e326 100644 --- a/docs/src/methods.md +++ b/docs/src/methods.md @@ -84,24 +84,36 @@ init_prop(state, generator, tlist, method::Val{:ExponentialUtilities}; kwargs... This method evaluates ``\exp(-i \op{H} dt) |Ψ⟩`` via a Krylov [`expv`](@extref ExponentialUtilities :jl:function:`ExponentialUtilities.expv`) -algorithm without explicitly forming the matrix exponential. It is therefore -often a good fit for larger systems or matrix-free operators where direct matrix -exponentiation is too costly. +algorithm without explicitly forming the matrix exponential. It builds a +Krylov subspace (via Arnoldi or Lanczos iteration) and then computes the +action of the matrix exponential on the state within that subspace. +This is often a good fit for larger systems or matrix-free operators where +direct matrix exponentiation is too costly. + +The propagator requires that states and operators satisfy the `AbstractArray` +interface (specifically, +[`supports_vector_interface`](@ref QuantumPropagators.Interfaces.supports_vector_interface) +for states and +[`supports_matrix_interface`](@ref QuantumPropagators.Interfaces.supports_matrix_interface) +for operators). **Advantages** * Avoids explicit construction of ``\op{U}`` * Works with matrix-free operators that support `mul!` * Good scaling for large sparse systems +* Supports both Hermitian (Lanczos) and non-Hermitian (Arnoldi) generators **Disadvantages** * Requires ExponentialUtilities.jl -* Performance depends on Krylov parameters and operator structure +* Performance depends on Krylov subspace dimension and operator structure +* Requires operators and states to implement the `AbstractArray` interface **When to use** * Large, sparse, or matrix-free generators +* Systems where ``\op{U}`` is too expensive to form explicitly ## [`Cheby`](@id method_cheby) diff --git a/ext/QuantumPropagatorsExponentialUtilitiesExt.jl b/ext/QuantumPropagatorsExponentialUtilitiesExt.jl index c19032a5..1a63ece4 100644 --- a/ext/QuantumPropagatorsExponentialUtilitiesExt.jl +++ b/ext/QuantumPropagatorsExponentialUtilitiesExt.jl @@ -35,12 +35,41 @@ expv_propagator = init_prop( ) ``` - initializes an [`ExponentialUtilitiesPropagator`](@ref). +initializes an [`ExponentialUtilitiesPropagator`](@ref). # Method-specific keyword arguments -* `expv_kwargs`: NamedTuple of keyword arguments forwarded to - `ExponentialUtilities.expv`. +* `expv_kwargs`: NamedTuple of keyword arguments forwarded to the underlying + [`ExponentialUtilities`](https://docs.sciml.ai/ExponentialUtilities/stable/) + routines. For in-place propagation, these are passed to both `arnoldi!` + (Krylov subspace construction) and `expv!` (exponentiation). For + out-of-place propagation, they are passed to `expv`, which forwards them + internally to `arnoldi`. The most useful keyword arguments are: + + * `m::Int`: Maximum Krylov subspace dimension. Defaults to + `min(30, size(A, 1))`. Larger values improve accuracy but increase memory + and computation cost per step. + * `ishermitian::Bool`: Whether the operator is Hermitian. Defaults to + `LinearAlgebra.ishermitian(A)`. When `true`, the cheaper Lanczos iteration + is used instead of the full Arnoldi process. Note that the propagator + passes `-𝕚 dt` as the time argument, so a Hermitian generator ``Ĥ`` still + results in a non-Hermitian product ``-𝕚 Ĥ dt``; however, + `ExponentialUtilities` handles the complex time scaling internally. + * `tol::Real`: Tolerance for happy breakdown. Defaults to `1e-7`. The + Arnoldi iteration stops early when the norm of the next Krylov vector + drops below `tol * opnorm`. + * `opnorm`: Operator norm of `A`. By default, this is computed + automatically. Supplying it avoids redundant norm computations. + * `iop::Int`: Incomplete orthogonalization procedure length. Defaults to `0` + (full orthogonalization). A positive value limits the number of previous + Krylov vectors used for orthogonalization, reducing cost at the expense of + numerical stability. + * `mode::Symbol`: Termination strategy for `expv`. Either + `:happy_breakdown` (default) or `:error_estimate`. In `:happy_breakdown` + mode, the iteration relies on early termination when the Krylov subspace + captures the relevant dynamics. In `:error_estimate` mode, an adaptive + step-size strategy is used with additional tolerance parameters `rtol` + (relative tolerance, defaults to `√tol`). """ function init_prop( state, @@ -70,12 +99,21 @@ function init_prop( if inplace A₀ = QuantumPropagators.Controls.evaluate(generator, tlist, n) + is_herm = get(expv_kwargs, :ishermitian, ishermitian(A₀)) + if haskey(expv_kwargs, :mode) && expv_kwargs[:mode] == :error_estimate && !is_herm + throw( + ArgumentError( + "`mode=:error_estimate` in `expv_kwargs` requires a " * + "Hermitian generator. Use the default " * + "`mode=:happy_breakdown` for non-Hermitian generators." + ) + ) + end Ks = ExponentialUtilities.arnoldi(A₀, state; expv_kwargs...) if haskey(expv_kwargs, :mode) && expv_kwargs[:mode] == :error_estimate cache = ExponentialUtilities.get_subspace_cache(Ks) else T = promote_type(eltype(A₀), eltype(state)) - is_herm = get(expv_kwargs, :ishermitian, ishermitian(A₀)) cache_type = is_herm ? real(T) : T cache = ExponentialUtilities.ExpvCache{cache_type}(Ks.m) end @@ -163,7 +201,7 @@ function prop_step!(propagator::ExponentialUtilitiesPropagator) propagator.expv_kwargs... ) end - setfield!(propagator, :state, Ψ) + setfield!(propagator, :state, convert(typeof(propagator.state), Ψ)) end _pwc_advance_time!(propagator) diff --git a/src/exponential_utilities_propagator.jl b/src/exponential_utilities_propagator.jl index 038c3229..d0d25fc7 100644 --- a/src/exponential_utilities_propagator.jl +++ b/src/exponential_utilities_propagator.jl @@ -7,20 +7,20 @@ This is a [`PWCPropagator`](@ref). Methods that depend on `ExponentialUtilities` are loaded via Julia package extensions. """ mutable struct ExponentialUtilitiesPropagator{GT,OT,ST,KST,CT} <: PWCPropagator - const generator::GT - state::ST - t::Float64 # time at which current `state` is defined - n::Int64 # index of next interval to propagate - const tlist::Vector{Float64} - parameters::AbstractDict - controls - genop::OT - Ks::KST - cache::CT - backward::Bool - inplace::Bool - expv_kwargs::NamedTuple - const timing_data::TimerOutput + const generator::GT + state::ST + t::Float64 # time at which current `state` is defined + n::Int64 # index of next interval to propagate + const tlist::Vector{Float64} + parameters::AbstractDict + controls + genop::OT + Ks::KST + cache::CT + backward::Bool + inplace::Bool + expv_kwargs::NamedTuple + const timing_data::TimerOutput end diff --git a/test/runtests.jl b/test/runtests.jl index 7b8bc702..60f5c5aa 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -78,8 +78,8 @@ using SafeTestsets include("test_expprop.jl") end - println("\n* ExponentialUtilities (test_exponential_utilities.jl):") - @time @safetestset "ExponentialUtilities" begin + println("\n* ExponentialUtilities (baseline) (test_exponential_utilities.jl):") + @time @safetestset "ExponentialUtilities (baseline)" begin include("test_exponential_utilities.jl") end @@ -98,8 +98,8 @@ using SafeTestsets include("test_propagate_sequence.jl") end - println("\n* ExponentialUtilities propagation (test_exputils.jl):") - @time @safetestset "ExponentialUtilities" begin + println("\n* ExponentialUtilities (propagation) (test_exputils.jl):") + @time @safetestset "ExponentialUtilities (propagation)" begin include("test_exputils.jl") end diff --git a/test/test_exponential_utilities.jl b/test/test_exponential_utilities.jl index ae1bd1c3..b1bcb818 100644 --- a/test/test_exponential_utilities.jl +++ b/test/test_exponential_utilities.jl @@ -600,6 +600,37 @@ end end +@testset "expv return type with real state and complex time" begin + + # When the state vector has real eltype but t is complex, expv returns a + # complex vector due to promote_type(typeof(t), eltype(A), eltype(b)). + # This verifies the type promotion behavior that motivates the convert() + # call in the non-inplace propagator path. + N = 10 + rng = StableRNG(2296342638) + A = random_matrix(N; rng, complex = false, hermitian = true) + b_real = real.(random_state_vector(N; rng)) + @test eltype(b_real) == Float64 + dt = 0.1 + t = -1im * dt + + # expv with complex t and real b returns Vector{ComplexF64} + w = ExponentialUtilities.expv(t, A, b_real) + @test eltype(w) == ComplexF64 + @test typeof(w) != typeof(b_real) + w_expected = exp(t * A) * b_real + @test norm(w - w_expected) < 1e-12 + + # expv with complex t and complex b returns the same type as the input + b_complex = complex.(b_real) + w2 = ExponentialUtilities.expv(t, A, b_complex) + @test eltype(w2) == ComplexF64 + @test typeof(w2) == typeof(b_complex) + @test norm(w2 - w_expected) < 1e-12 + +end + + @testset "expv with OffsetArray state vector (broken)" begin # expv does not support OffsetVector as state vector b: arnoldi builds From a11a0d18603eca893370cec0621a6378cee251c8 Mon Sep 17 00:00:00 2001 From: Michael Goerz Date: Thu, 16 Apr 2026 11:54:46 -0400 Subject: [PATCH 12/13] Minor tweaks to the documentation --- docs/src/methods.md | 75 +++++++++++++++++++++++--------------------- docs/src/overview.md | 12 ++++--- test/runtests.jl | 4 +-- 3 files changed, 50 insertions(+), 41 deletions(-) diff --git a/docs/src/methods.md b/docs/src/methods.md index a875e326..f1cf31a1 100644 --- a/docs/src/methods.md +++ b/docs/src/methods.md @@ -15,14 +15,17 @@ As discussed in the [Overview](@ref overview_approaches), time propagation can b We consider this especially in the piecewise-constant case (`pwc=true` in [`propagate`](@ref)/[`init_prop`](@ref)), which is required for the traditional optimization methods [GRAPE](https://juliaquantumcontrol.github.io/GRAPE.jl/stable/) and [Krotov](https://juliaquantumcontrol.github.io/Krotov.jl/stable/). In these propagations, the time-dependent generator ``\op{H}(t)`` is [evaluated](@ref QuantumPropagators.Controls.evaluate) to a constant operator ``\op{H}`` on each interval of the time grid. The analytical solution to the Schrödinger or Liouville equation is well known, and propagation step simply has to evaluate the application of the time evolution operator ``\op{U} = \exp[-i \op{H} dt]`` to the state ``|Ψ⟩``. The following methods are built in to `QuantumPropagators`: * [`ExpProp`](@ref method_expprop) – constructs ``\op{U}`` explicitly and then applies it to ``|Ψ⟩`` - * [`ExponentialUtilities`](@ref method_exponentialutilities) – applies ``\op{U} |Ψ⟩`` using Krylov methods without explicitly forming ``\op{U}`` * [`Cheby`](@ref method_cheby) — expansion of ``\op{U} |Ψ⟩`` into Chebychev polynomials, valid if ``\op{H}`` has real eigenvalues * [`Newton`](@ref method_newton) – expansion of ``\op{U} |Ψ⟩`` into Newton polynomials, valid if ``\op{H}`` has complex eigenvalues (non-Hermitian Hamiltonian, Liouvillian) - The `ExpProp` method is generally not numerically efficient, but works well for small system for for debugging. The two "core" methods based on a polynomials series expansions are more suitable for bigger systems and provide both efficiency and high precision (in general, the is truncated as soon as some desired precision is reached, which is machine precision by default). + The `ExpProp` method is generally not numerically efficient, but works well for small system for for debugging. The two "core" methods based on a polynomial series expansions are more suitable for bigger systems and provide both efficiency and high precision (in general, the series is truncated as soon as some desired precision is reached, which is machine precision by default). Note that this high precision is *within the piecewise-constant approximation*. The discretization itself may introduce a non-negligible error compared to the time-continuous dynamics. There is tradeoff: A smaller `dt` decreases the discretization error, but the polynomial expansions are more effective with larger time steps. + There is also support for _external_ packages to efficiently evaluate the application of the piecewise-constant time evolution operator: + + * [`ExponentialUtilities`](@ref method_exponentialutilities) – applies ``\op{U} |Ψ⟩`` using Krylov methods without explicitly forming ``\op{U}`` + 2. By solving the equation of motion explicitly with an ODE solver. We support the use of any of the ODE solvers in [OrdinaryDiffEq.jl](https://docs.sciml.ai/OrdinaryDiffEq/stable/): @@ -66,6 +69,38 @@ init_prop(state, generator, tlist, method::Val{:ExpProp}; kwargs...) * Comparing against another propagator +## [`Cheby`](@id method_cheby) + +The method should be loaded with + +``` +using QuantumPropagators: Cheby +``` + +and then passed as `method=Cheby` to [`propagate`](@ref) or [`init_prop`](@ref): + +```@docs +init_prop(state, generator, tlist, method::Val{:Cheby}; kwargs...) +``` + +The time evolution operator of the piecewise-constant Schrödinger equation ``|Ψ(t)⟩ = e^{-i Ĥ dt} |Ψ(0)⟩`` is evaluated by an expansion into Chebychev polynomials [Tal-EzerJCP1984, KosloffJCP1988](@cite). This requires ``Ĥ`` to be Hermitian (have real eigenvalues) and to have a known spectral range, so that it can be normalized to the domain ``[-1, 1]`` on which the Chebychev polynomials are defined. + +See [GoerzPhd2015; Chapter 3.2.1](@cite) for a detailed description of the method. + +**Advantages** + +* Very efficient for high precision and moderately large time steps + +**Disadvantages** + +* Only valid for Hermitian operators +* Must be able to estimate the spectral envelope + +**When to use** + +* Closed quantum systems with piecewise constant dynamics + + ## [`ExponentialUtilities`](@id method_exponentialutilities) The method requires that the [ExponentialUtilities.jl](https://docs.sciml.ai/ExponentialUtilities/stable/) @@ -108,7 +143,9 @@ for operators). * Requires ExponentialUtilities.jl * Performance depends on Krylov subspace dimension and operator structure -* Requires operators and states to implement the `AbstractArray` interface +* Requires operators and states to implement an array interface, see + [`QuantumPropagators.Interfaces.supports_matrix_interface`](@ref) and + [`QuantumPropagators.Interfaces.supports_vector_interface`](@ref), respectively **When to use** @@ -116,38 +153,6 @@ for operators). * Systems where ``\op{U}`` is too expensive to form explicitly -## [`Cheby`](@id method_cheby) - -The method should be loaded with - -``` -using QuantumPropagators: Cheby -``` - -and then passed as `method=Cheby` to [`propagate`](@ref) or [`init_prop`](@ref): - -```@docs -init_prop(state, generator, tlist, method::Val{:Cheby}; kwargs...) -``` - -The time evolution operator of the piecewise-constant Schrödinger equation ``|Ψ(t)⟩ = e^{-i Ĥ dt} |Ψ(0)⟩`` is evaluated by an expansion into Chebychev polynomials [Tal-EzerJCP1984, KosloffJCP1988](@cite). This requires ``Ĥ`` to be Hermitian (have real eigenvalues) and to have a known spectral range, so that it can be normalized to the domain ``[-1, 1]`` on which the Chebychev polynomials are defined. - -See [GoerzPhd2015; Chapter 3.2.1](@cite) for a detailed description of the method. - -**Advantages** - -* Very efficient for high precision and moderately large time steps - -**Disadvantages** - -* Only valid for Hermitian operators -* Must be able to estimate the spectral envelope - -**When to use** - -* Closed quantum systems with piecewise constant dynamics - - ## [`Newton`](@id method_newton) The method should be loaded with diff --git a/docs/src/overview.md b/docs/src/overview.md index 25d9fe9f..c231f172 100644 --- a/docs/src/overview.md +++ b/docs/src/overview.md @@ -77,7 +77,7 @@ This form differs from most textbooks by a factor of ``i``, but has the benefit There are two fundamental approaches to solving the Schrödinger equation (or any equation of motion): -1. We can analytically solve the Schrödinger equation and then numerically evaluate the solution. Mathematically, this is the application of the time evolution operator ``\op{H}`` as ``|Ψ(t+dt)⟩ = \op{U}(t) |Ψ(t)⟩``. For a piecewise-constant ``\op{H}(t)``where there is a time-independent ``\op{H}`` in the interval ``[t, t+dt]``, the time evolution operator is well-known to be +1. We can analytically solve the Schrödinger equation and then numerically evaluate the solution. Mathematically, this is the application of the time evolution operator ``\op{U}`` as ``|Ψ(t+dt)⟩ = \op{U}(t) |Ψ(t)⟩``. For a piecewise-constant ``\op{H}(t)``where there is a time-independent ``\op{H}`` in the interval ``[t, t+dt]``, the time evolution operator is well-known to be ```math \op{U} = \exp[-i \op{H} dt]\,. @@ -112,7 +112,7 @@ Furthermore, as a fallback for very small system or for debugging, * `method=ExpProp`: Explicitly construct the time evolution operator by matrix exponentiation and apply it to the state. -If the [OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl) or (equivalently) th [DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/) package is loaded, `QuantumPropagators` can delegate to it: +If the [OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl) or (equivalently) the [DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/) package is loaded, `QuantumPropagators` can delegate to it: ```julia using OrdinaryDiffEq @@ -122,7 +122,11 @@ allows to pass * `method=OrdinaryDiffEq`: Use [OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl) as a backend with any of the [algorithms available for ODEs in DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/), `alg=Tsit5()` by default. -Unlike any of the built-in methods, `OrdinaryDiffEq` is able to propagate for time-continuous generators. This is the default for that propagator (`pwc=false`). By setting `pwc=true` or `piecewise=true`) the ODE solvers can also be used for piecewise-constant Hamiltonians or Liouvillians, providing an alternative to the built-in `method=Cheby` and `method=Newton`. +Unlike any of the built-in methods, `OrdinaryDiffEq` is able to propagate for time-continuous generators. This is the default for that propagator (`pwc=false`). By setting `pwc=true` or `piecewise=true`, the ODE solvers can also be used for piecewise-constant Hamiltonians or Liouvillians, providing an alternative to the built-in `method=Cheby` and `method=Newton`. + +If the [ExponentialUtilities.jl](https://docs.sciml.ai/ExponentialUtilities/stable/expv/) package is loaded, this enables + +* `method=ExponentialUtilities`: Evaluate the application of the piecewise-constant time evolution operator via [`ExponentialUtilities.expv`](@extref). See the more extended discussion of [Propagation Methods](@ref) for more details. @@ -182,7 +186,7 @@ Most propagators support both an in-place and a not-in-place mode. These modes c In-place operations can be dramatically more efficient for large Hilbert space dimensions. On the other hand, not-in-place operations can be more efficient for small Hilbert spaces, in particular when a [static vector](@extref StaticArrays `SVector`) can be used to represent the state. Moreover, frameworks for automatic differentiation such as [Zygote](https://fluxml.ai/Zygote.jl/stable/) do not support in-place operations. -When using custom structs for states, operators, or generators, the struct itself not not need to be mutable (according to [`Base.ismutable`](@extref Julia)) in order to support `inplace=true`. It only must support the in-place operations defined in the [formal interface](@ref QuantumPropagatorsInterfacesAPI) and indicate that support by defining [`QuantumPropagators.Interfaces.supports_inplace`](@ref). +When using custom structs for states, operators, or generators, the struct itself does not need to be mutable (according to [`Base.ismutable`](@extref Julia)) in order to support `inplace=true`. It only must support the in-place operations defined in the [formal interface](@ref QuantumPropagatorsInterfacesAPI) and indicate that support by defining [`QuantumPropagators.Interfaces.supports_inplace`](@ref). Typically, in-place operations on immutable custom structs involve mutating the mutable properties of that struct. diff --git a/test/runtests.jl b/test/runtests.jl index 60f5c5aa..3268b7ea 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -78,7 +78,7 @@ using SafeTestsets include("test_expprop.jl") end - println("\n* ExponentialUtilities (baseline) (test_exponential_utilities.jl):") + println("\n* ExponentialUtilities baseline (test_exponential_utilities.jl):") @time @safetestset "ExponentialUtilities (baseline)" begin include("test_exponential_utilities.jl") end @@ -98,7 +98,7 @@ using SafeTestsets include("test_propagate_sequence.jl") end - println("\n* ExponentialUtilities (propagation) (test_exputils.jl):") + println("\n* ExponentialUtilities propagation (test_exputils.jl):") @time @safetestset "ExponentialUtilities (propagation)" begin include("test_exputils.jl") end From a1397692ad390c89a265e986a99a1637f0463433 Mon Sep 17 00:00:00 2001 From: Michael Goerz Date: Thu, 16 Apr 2026 13:04:40 -0400 Subject: [PATCH 13/13] Add test for Gradgen propagation --- test/test_exputils.jl | 28 +++++++++++++++++++++++++++- 1 file changed, 27 insertions(+), 1 deletion(-) diff --git a/test/test_exputils.jl b/test/test_exputils.jl index c2f88540..68cde686 100644 --- a/test/test_exputils.jl +++ b/test/test_exputils.jl @@ -3,7 +3,9 @@ using LinearAlgebra using SparseArrays using QuantumPropagators using QuantumPropagators.Generators: liouvillian +using QuantumPropagators: Cheby using ExponentialUtilities +using QuantumGradientGenerators: GradGenerator, GradVector using StableRNGs: StableRNG @@ -143,4 +145,28 @@ end @test norm(Ψ_out_expv - Ψ_out_exp) < 1e-9 end -# TODO: test a GradGen propagation +@testset "GradGen propagation" begin + Ψ₀ = ComplexF64[1, 0] + H₀ = ComplexF64[0 0; 0 1] + H₁ = ComplexF64[0 1; 1 0] + + tlist = collect(range(0, 10.0, length = 201)) + ϵ = 0.2 * ones(length(tlist)) + + generator = hamiltonian(H₀, (H₁, ϵ)) + G̃ = GradGenerator(generator) + Ψ̃₀ = GradVector(Ψ₀, 1) + + Ψ̃_out = propagate( + Ψ̃₀, + G̃, + tlist; + method = ExponentialUtilities, + inplace = false, + expv_kwargs = (; ishermitian = false, tol = 1e-12) + ) + Ψ_cheby = propagate(Ψ₀, generator, tlist; method = Cheby) + + @test norm(Ψ̃_out.state - Ψ_cheby) < 1e-7 + @test norm(Ψ̃_out.grad_states[1]) > 0 +end