CTFlows.jl

Appearance

Flow(Hamiltonian) compatibility ​

Which state/costate types and call styles the Hamiltonian flow built from a scalar Data.Hamiltonian accepts, on CPU with the default AD backend (AutoForwardDiff). Every cell has an executed example; see Compatibility for how the pages stay true. For devices, see GPU flows.

Unlike Flow(VectorField) and Flow(HamiltonianVectorField), a scalar Hamiltonian has no in-place variant — there is no mutability trait on Data.Hamiltonian at all. The flow computes  ,    by automatic differentiation instead of from a user-supplied vector field, so the axis that varies on this page is not OOP/IP but the input type accepted by the AD backend.

All examples use the harmonic oscillator   (autonomous, fixed) — the same dynamics as Flow(HamiltonianVectorField), so the analytic solution is identical:    at  .

The flow is built once — no ad_backend keyword needed, since AutoForwardDiff is already the CPU default:

julia
h     = Data.Hamiltonian((x, p) -> 0.5 * (sum(abs2, x) + sum(abs2, p)))
hflow = Flows.Flow(h; reltol=1e-8)
Flow
├─ system: HamiltonianSystem
│  ├─ time_dependence: Autonomous
│  ├─ variable_dependence: Fixed
│  ├─ Hamiltonian: autonomous, fixed (no variable)
│  │    natural call: h(x, p)
│  │    uniform call: h(t, x, p, v)
│  └─ backend: DifferentiationInterface{CPU}(ad_backend=AutoForwardDiff)
└─ integrator: SciML{CPU} (instance, id=:sciml)
   ├─ internalnorm = real_norm  [default]
   ├─ alg = Tsit5  [default]
   ├─ reltol = 1.0e-8  [user]
   ├─ save_everystep = auto  [default]
   ├─ abstol = 1.0e-8  [default]
   ├─ save_start = auto  [default]
   └─ dense = auto  [default]
   Tip: use describe(SciML{CPU}) to see all available options.

An explicit backend, or a GPU one, can be selected the same way: Flows.Flow(h; ad_backend=ADTypes.AutoForwardDiff()) — see GPU flows for method=:gpu and AutoMooncake/AutoZygote. On CPU, AutoMooncake also works and pays off for large states (a hundred dimensions and more), but it does not support the nested differentiation shown below.

Compatibility table ​

The two call styles are the point call hflow(t0, x0, p0, tf) → final (xf, pf), and the trajectory call hflow((t0, tf), x0, p0) → Trajectories.HamiltonianVectorFieldTrajectory. There is no OOP/IP split.

State/costate typepointtraj
Scalar Real✓✓
Vector Real✓✓
MVector Real✓✓
SVector Real✓✓
Matrix Real (batch)✓✓
MMatrix Real✓✓
SMatrix Real✓✓
Complex (any container)✗ (a)✗ (a)

Legend ​

  • ✓ — works; the result matches the analytic    and is executed on this page.

  • ✗ — raises an error; see the corresponding note.

(a) Complex is not supported by the default backend

AutoForwardDiff cannot differentiate through a Complex input: constructing a dual number over ComplexF64 raises ArgumentError: Cannot create a dual over scalar type ComplexF64. This is a ForwardDiff.jl limitation, not specific to CTFlows.

If your Hamiltonian naturally involves complex state (e.g. from a Schrödinger-type equation), it can typically be rewritten as a real system by splitting the state into real/imaginary (or conjugate) components — this is not a limitation of CTFlows, just a modeling choice that avoids differentiating through Complex directly.

The table above is about the flow's state/costate input. A separate question — how to differentiate the flow itself with respect to (x0, p0), e.g. for a shooting method — is covered next.

Real states ​

Scalar ​

julia
julia> hflow(0.0, 1.0, 0.0, pi/2)   # (xf, pf) ≈ (0.0, -1.0)
(-1.554201241359274e-10, -0.9999999997009265)
julia
julia> sol = hflow((0.0, pi/2), 1.0, 0.0);

julia> (Trajectories.state(sol)(pi/2), Trajectories.costate(sol)(pi/2))   # both scalars
(-1.554197798594726e-10, -0.9999999997009265)

Vector ​

julia
julia> hflow(0.0, [1.0, 0.0], [0.0, 1.0], pi/2)
([-1.554201241359274e-10, 0.9999999997009265], [-0.9999999997009265, -1.554201241359274e-10])

MVector and SVector ​

julia
julia> hflow(0.0, MVector{2}(1.0, 0.0), MVector{2}(0.0, 1.0), pi/2)
([-1.554201241359274e-10, 0.9999999997009265], [-0.9999999997009265, -1.554201241359274e-10])

julia> hflow(0.0, SA[1.0, 0.0], SA[0.0, 1.0], pi/2)
([-1.554201241359274e-10, 0.9999999997009265], [-0.9999999997009265, -1.554201241359274e-10])

Matrix (batch) ​

H is separable and quadratic (sum(abs2, x)), so it reduces a matrix state/costate to a scalar the same way it does a vector, and the induced dynamics are independent per column:

julia
julia> hflow(0.0, [1.0 2.0; 3.0 4.0], [0.0 0.0; 1.0 1.0], pi/2)
([-9.327918151763014e-11 -1.8655836303526027e-10; 0.9999999995245169 0.9999999994312374], [-0.999999999804355 -1.99999999960871; -2.999999999506345 -3.9999999993106985])

julia> X0m = MMatrix{2,2}(1.0, 3.0, 2.0, 4.0)
2×2 MMatrix{2, 2, Float64, 4} with indices SOneTo(2)×SOneTo(2):
 1.0  2.0
 3.0  4.0

julia> P0m = MMatrix{2,2}(0.0, 1.0, 0.0, 1.0)
2×2 MMatrix{2, 2, Float64, 4} with indices SOneTo(2)×SOneTo(2):
 0.0  0.0
 1.0  1.0

julia> X0s = SMatrix{2,2}(1.0, 3.0, 2.0, 4.0)
2×2 SMatrix{2, 2, Float64, 4} with indices SOneTo(2)×SOneTo(2):
 1.0  2.0
 3.0  4.0

julia> P0s = SMatrix{2,2}(0.0, 1.0, 0.0, 1.0)
2×2 SMatrix{2, 2, Float64, 4} with indices SOneTo(2)×SOneTo(2):
 0.0  0.0
 1.0  1.0

julia> hflow(0.0, X0m, P0m, pi/2)   # mutable
([-9.327918151763014e-11 -1.8655836303526027e-10; 0.9999999995245169 0.9999999994312374], [-0.999999999804355 -1.99999999960871; -2.999999999506345 -3.9999999993106985])

julia> hflow(0.0, X0s, P0s, pi/2)   # immutable
([-9.327918151763014e-11 -1.8655836303526027e-10; 0.9999999995245169 0.9999999994312374], [-0.999999999804355 -1.99999999960871; -2.999999999506345 -3.9999999993106985])

Complex states — not supported ​

julia
julia> hflow(0.0, 1.0 + 2.0im, 0.0 + 0.0im, pi/2)
ERROR: ArgumentError: Cannot create a dual over scalar type ComplexF64. If the type behaves as a scalar, define ForwardDiff.can_dual(::Type{ComplexF64}) = true.

See note (a) above — for a system that is naturally complex, rewrite it as a real system on the real/imaginary (or conjugate) components instead.

Automatic differentiation: sensitivities of the flow ​

Since Flow(Hamiltonian) is itself AutoForwardDiff-backed internally, differentiating it with respect to (x0, p0) — e.g. for a shooting method — means nested AD. Do this by differentiating the flow call, never by constructing a ForwardDiff.Dual by hand and passing it as x0/p0: a hand-built, untagged Dual collides with the flow's own internal tag (DualMismatchError), and a hand-supplied custom tag does not fix it either — an arbitrary marker type lacks ForwardDiff.Tag's internal ordering machinery. Wrapping the call in an outer ForwardDiff.jacobian / ForwardDiff.gradient sidesteps this entirely, because ForwardDiff then generates its own properly ordered tag for the perturbation — this is exactly how a NonlinearSolve-style shooting method differentiates a flow that is itself ForwardDiff-backed.

This computes the exact Jacobian of with respect to — for the linear harmonic oscillator it is the constant matrix :

julia
julia> shoot(z) = collect(hflow(0.0, z[1], z[2], pi/2))
shoot (generic function with 1 method)

julia> ForwardDiff.jacobian(shoot, [1.0, 0.0])
2×2 Matrix{Float64}:
 -1.5542e-10   1.0
 -1.0         -1.5542e-10

The same closure-wrapping pattern works with ForwardDiff.gradient for a scalar shooting residual. No internal backend change (e.g. switching to AutoZygote/AutoMooncake to dodge the nested Dual) is needed for this to work. The other way round does not work: a flow built with ad_backend=AutoMooncake() refuses an outer ForwardDiff with a PreconditionError, since Mooncake cannot differentiate through Dual numbers.

See also ​