Add ExponentialUtilities propagator

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 `<span class="tlaUNC">` 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.

Project.toml

@@ -16,17 +16,20 @@ 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"
 RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
 [extensions]
+QuantumPropagatorsExponentialUtilitiesExt = "ExponentialUtilities"
 QuantumPropagatorsODEExt = "OrdinaryDiffEq"
 QuantumPropagatorsRecursiveArrayToolsExt = "RecursiveArrayTools"
 QuantumPropagatorsStaticArraysExt = "StaticArrays"
 
 [compat]
 ArrayInterface = "7.0"
+ExponentialUtilities = "1.11"
 OffsetArrays = "1"
 OrdinaryDiffEq = "6.59"
 ProgressMeter = "1"

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"

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
@@ -46,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/",
 )
 

docs/src/methods.md

@@ -18,10 +18,14 @@ As discussed in the [Overview](@ref overview_approaches), time propagation can b
    * [`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/):
@@ -97,6 +101,58 @@ See [GoerzPhd2015; Chapter 3.2.1](@cite) for a detailed description of the metho
 * Closed quantum systems with piecewise constant dynamics
 
 
+## [`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` 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`](@extref ExponentialUtilities :jl:function:`ExponentialUtilities.expv`)
+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 subspace dimension and operator structure
+* 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**
+
+* Large, sparse, or matrix-free generators
+* Systems where ``\op{U}`` is too expensive to form explicitly
+
+
 ## [`Newton`](@id method_newton)
 
 The method should be loaded with

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.