Add ReducedHamiltonianLoss: fit a decoder to states and reduced-Hamiltonian gradients

CHANGELOG.md

@@ -17,6 +17,7 @@ breaking release).
 
 ### Added
 
+- **`ReducedHamiltonianLoss` fits a decoder to states and to the gradient of the reduced Hamiltonian.** For latent points ``z`` with targets ``x^*`` and ``g^*``, the loss is ``\lVert\Psi^\mathrm{dec}(z)-x^*\rVert/\lVert x^*\rVert + \lambda\lVert\nabla_z(H\circ\Psi^\mathrm{dec})(z)-g^*\rVert/\lVert g^*\rVert``, for a Hamiltonian `H` on the full space that the user supplies as a function of a matrix of states. For a symplectic decoder the second term is the reduced vector field up to ``\mathbb{J}``, so it fits the reduced dynamics, critical points and separatrices included, which a fit of the states alone does not control. The gradient in ``z`` is taken by central differences (step `h`), so the parameter gradient needs no nested automatic differentiation. Where a target vanishes on the whole batch (``g^* = 0`` where every latent point is a critical point of the target), its term is the absolute error instead of a relative one. It accepts matrices and the 3-tensors of a `DataLoader` built from an input and an output matrix. `test/loss/reduced_hamiltonian_loss.jl` checks the central differences against Zygote, that a vanishing target gradient gives a finite loss, that a decoder has a vanishing loss on its own targets, the weight `λ`, the Zygote gradient, and that training one `SymplecticAutoencoder` decoder on another's targets halves the loss in 200 epochs.
 - **The test suite runs on an Apple GPU in its `metal` group.** `test/devices/metal.jl` checks the two device paths that were measured by hand and not tested before: the products of a `PoissonTensor` with a wrapped device array, `𝕁 * view(A, [1, 2, 3, 4], :)` and `𝕁 * B'`, which `ext/GPUArraysCoreExt.jl` carries; and the `tensor_mat_mul!` kernel on `MetalBackend()`. Each result must stay on the device and match the host in `Float32`, with scalar indexing off. `Pkg.test(test_args = ["metal"])` runs them without the subject tests, and a missing device records a skip, which happens inside a sandbox. Metal (1.10 or later) is a test dependency on every platform; it installs and precompiles on Linux and Windows, and nothing there loads it. On an Apple M4 Max with Metal 1.11.1 all five assertions pass.
 - **A `Metal` workflow runs the Metal tests on GitHub's `macos-15` runner**, for the `min` and `1` Julia versions. The runner's GPU is Apple's paravirtual device, which Metal.jl supports from 1.10 and on macOS 15 or later. The workflow includes a step "Metal is functional" before the tests that fails the job where the device is unavailable, so it cannot pass without having run the tests. The job is not a required check. However, the Metal tests also run in two required checks—CI.yml's `Julia min - macOS-latest - default` and `Julia 1 - macOS-latest - default`—which execute `test/runtests.jl` with no test arguments on Apple silicon, selecting the `metal` group.
 - **`test/quality/jet.jl` checks the kernel launchers with JET.** It runs in `core`, directly after `quality/aqua.jl`, and takes about 14 s. Each function that launches a kernel of `src/` on the `CPU()` backend, the pullbacks among them, and the `Batch` functor that `test/data_loader/batch_index_set.jl` measures with `@allocated` get one `JET.report_opt` line per element type that the tests use, with `target_modules = (GeometricMachineLearning,)`. Each `@kernel` method also gets a line that analyses its generated `cpu_<kernel>` function at the launcher's context, because JET drops the reports of a kernel body when it analyses the launcher. 142 lines pass on Julia 1.11.9 (JET 0.9.20) and 1.13.1 (JET 0.12.2). Three lines report a runtime dispatch on both versions and are `@test_broken` with issue #325: `PoissonTensor(::CPU, ::Int, ::DataType)`, and the input-and-output method of `convert_input_and_batch_indices_to_array` at `Float32` and `Float64`. Where JET does not work, the file records one `@test_skip`. JET joins `test/Project.toml` with no `[compat]` bound. `KNOWN_ISSUES.md` K2 records the Revise `EMFILE` messages that JET brings into the test log, and K3 the dispatches that the lines do not see: an unstable argument of a kernel launch, and two barriers in the Poisson-tensor kernel.

docs/src/reduced_order_modeling/losses.md

@@ -40,12 +40,15 @@ e_\mathrm{red}(\mu) := \sqrt{
 ```
 where ``\mathbf{x}^{(t)}`` is the solution of the FOM at point ``t`` and ``\mathbf{x}^{(t)}_r`` is the solution of the ROM (in the reduced basis) at point ``t``. The reduction error, as opposed to the projection error, not only measures how well the solution manifold is represented by the reduced basis, but also measures how well the FOM dynamics are approximated by the ROM dynamics (via the induced vector field on the reduced basis). The corresponding function in `GeometricMachineLearning` is [`reduction_error`](@ref). The reduction error is, in contract to the projection error, typically not used during training (even though some authors are using a similar error to do so [lee2020model](@cite)).
 
+For a decoder whose reduced dynamics matter more than the pointwise reconstruction, [`ReducedHamiltonianLoss`](@ref) fits the states together with the gradient of the reduced Hamiltonian ``H\circ\Psi^\mathrm{dec}``, which for a symplectic decoder is the reduced vector field up to ``\mathbb{J}``.
+
 ## Library Functions
 
 ```@docs
 TransformerLoss
 AutoEncoderLoss
 ReducedLoss
+ReducedHamiltonianLoss
 projection_error
 reduction_error
 ```

src/GeometricMachineLearning.jl

@@ -222,7 +222,7 @@ export AdamOptimizerWithDecay, DecayingStatic
 export NeuralNetwork
 
 export NetworkLoss, TransformerLoss, FeedForwardLoss, AutoEncoderLoss, ReducedLoss, HNNLoss,
-       LNNLoss, SymplecticEulerLoss, VariationalMidpointLoss
+       LNNLoss, SymplecticEulerLoss, VariationalMidpointLoss, ReducedHamiltonianLoss
 
 #INCLUDE ARCHITECTURES
 include("architectures/neural_network_integrator.jl")
@@ -257,6 +257,7 @@ include("loss/hnn_loss.jl")
 include("loss/lnn_loss.jl")
 include("loss/symplectic_euler_loss.jl")
 include("loss/variational_midpoint_loss.jl")
+include("loss/reduced_hamiltonian_loss.jl")
 
 export AbstractPullback, ZygotePullback, SymbolicPullback
 include("pullbacks/zygote_pullback.jl")

src/loss/reduced_hamiltonian_loss.jl

@@ -0,0 +1,101 @@
+@doc raw"""
+    ReducedHamiltonianLoss(hamiltonian; λ = 1, h = 1e-4)
+
+Make an instance of `ReducedHamiltonianLoss` for a Hamiltonian ``H`` on the full space.
+
+The loss fits a [`Decoder`](@ref) ``\Psi^\mathrm{dec}:\mathbb{R}^{2n}\to\mathbb{R}^{2N}`` on latent points ``z`` to target states ``x^*`` and to target gradients ``g^*`` of the *reduced Hamiltonian* ``H\circ\Psi^\mathrm{dec}``:
+
+```math
+L = \frac{\lVert\Psi^\mathrm{dec}(z) - x^*\rVert}{\lVert x^*\rVert}
+  + \lambda\,\frac{\lVert\nabla_z(H\circ\Psi^\mathrm{dec})(z) - g^*\rVert}{\lVert g^*\rVert},
+```
+
+with the absolute error in place of a relative one whose target vanishes on the batch (``g^* = 0`` where every latent point of the batch is a critical point of the target).
+
+For a symplectic decoder, as in a [`SymplecticAutoencoder`](@ref), the reduced vector field is ``\mathbb{J}_{2n}\nabla_z(H\circ\Psi^\mathrm{dec})``, so the second term fits the reduced dynamics directly: their orbits are the level sets of ``H\circ\Psi^\mathrm{dec}``, and a small error in its gradient gives nearly the same phase portrait, its critical points and separatrices included. A fit of the states alone can be small while this gradient is not, in particular where ``H\circ\Psi^\mathrm{dec}`` is steep.
+
+The gradient with respect to ``z`` is computed by central differences with step `h`, from ``2\cdot2n + 1`` evaluations of the decoder on the batch. This keeps the loss differentiable with respect to the parameters by the same reverse-mode pass as any other loss, without nesting automatic differentiation.
+
+# Arguments
+
+- `hamiltonian`: a function that maps a matrix of states in ``\mathbb{R}^{2N}``, one per column, to the vector of their energies. It has to be differentiable by Zygote.
+- `λ = 1`: the weight of the gradient term.
+- `h = 1e-4`: the step of the central differences, in latent coordinates.
+
+# Functor
+
+```julia
+loss(model, ps, input, output)
+loss(nn, input, output)
+```
+
+`input` holds the latent points ``z``, one per column (``2n`` rows). `output` stacks the targets: the first ``2N`` rows are ``x^*``, the last ``2n`` rows are ``g^*``.
+
+# Examples
+
+A decoder fitted to its own states and reduced gradients has a loss of the size of the difference error:
+
+```jldoctest
+using GeometricMachineLearning
+using Zygote: gradient
+import Random
+Random.seed!(123)
+
+H(X) = vec(sum(abs2, X[3:4, :]; dims = 1)) ./ 2 .+ X[2, :]   # |p|²/2 + q₂ on ℝ⁴
+dec = decoder(NeuralNetwork(SymplecticAutoencoder(4, 2; n_encoder_blocks = 2, n_decoder_blocks = 2)))
+z = randn(2, 20)
+x = dec(z)
+g = reduce(hcat, [gradient(ζ -> H(reshape(dec(ζ), 4, 1))[1], collect(ζ))[1] for ζ in eachcol(z)])
+loss = ReducedHamiltonianLoss(H)
+loss(dec, z, vcat(x, g)) < 1e-6
+
+# output
+
+true
+```
+"""
+struct ReducedHamiltonianLoss{HT, T} <: NetworkLoss
+    hamiltonian::HT
+    λ::T
+    h::T
+end
+
+function ReducedHamiltonianLoss(hamiltonian; λ::Real = 1, h::Real = 1e-4)
+    T = float(promote_type(typeof(λ), typeof(h)))
+    ReducedHamiltonianLoss(hamiltonian, T(λ), T(h))
+end
+
+# The relative error, and the absolute one where the target vanishes: a batch of latent points that
+# are all critical points of the target has g* = 0.
+_relative_error(prediction, target) = (n = norm(target); iszero(n) ? norm(prediction) : norm(prediction - target) / n)
+
+# The gradient of `H∘model` with respect to each column of `z`, by central differences.
+function _reduced_hamiltonian_gradient(loss::ReducedHamiltonianLoss, model, ps, z::AbstractMatrix)
+    n = size(z, 1)
+    h = convert(eltype(z), loss.h)
+    rows = map(1:n) do i
+        # the shift is a column broadcast over the batch; built without mutation for Zygote
+        δ = map(j -> j == i ? h : zero(h), 1:n)
+        Hp = loss.hamiltonian(model(z .+ δ, ps))
+        Hm = loss.hamiltonian(model(z .- δ, ps))
+        reshape((Hp .- Hm) ./ (2h), 1, :)
+    end
+    reduce(vcat, rows)
+end
+
+function (loss::ReducedHamiltonianLoss)(model::Union{Chain, AbstractExplicitLayer}, ps,
+        input::AbstractMatrix, output::AbstractMatrix)
+    n = size(input, 1)
+    x_target = output[1:(end - n), :]
+    g_target = output[(end - n + 1):end, :]
+    state_error = _relative_error(model(input, ps), x_target)
+    gradient_error = _relative_error(_reduced_hamiltonian_gradient(loss, model, ps, input), g_target)
+    state_error + loss.λ * gradient_error
+end
+
+# A `DataLoader` built from an input and an output matrix holds them as tensors with a trailing
+# dimension of one; the columns are independent latent points either way.
+function (loss::ReducedHamiltonianLoss)(model::Union{Chain, AbstractExplicitLayer}, ps,
+        input::AbstractArray{T, 3}, output::AbstractArray{T, 3}) where {T}
+    loss(model, ps, reshape(input, size(input, 1), :), reshape(output, size(output, 1), :))
+end

test/loss/reduced_hamiltonian_loss.jl

@@ -0,0 +1,80 @@
+using GeometricMachineLearning
+using GeometricMachineLearning: _reduced_hamiltonian_gradient, params
+using Zygote: Zygote
+using Test
+using LinearAlgebra: norm
+import Random
+
+Random.seed!(1234)
+
+H(X) = vec(sum(abs2, X[3:4, :]; dims = 1)) ./ 2 .+ X[2, :]          # |p|²/2 + q₂ on ℝ⁴
+reduced_gradient(dec, z) = reduce(hcat, [Zygote.gradient(ζ -> H(reshape(dec(ζ), 4, 1))[1], collect(ζ))[1] for ζ in eachcol(z)])
+targets(dec, z) = vcat(dec(z), reduced_gradient(dec, z))
+
+# The central differences agree with automatic differentiation.
+function test_gradient()
+    dec = decoder(NeuralNetwork(SymplecticAutoencoder(4, 2; n_encoder_blocks = 2, n_decoder_blocks = 2)))
+    z = randn(2, 30)
+    loss = ReducedHamiltonianLoss(H)
+    G = _reduced_hamiltonian_gradient(loss, dec.model, params(dec), z)
+    @test G ≈ reduced_gradient(dec, z) rtol = 1e-6
+end
+
+# A decoder fitted to its own targets has a loss of the size of the difference error.
+function test_self_consistency()
+    dec = decoder(NeuralNetwork(SymplecticAutoencoder(4, 2; n_encoder_blocks = 2, n_decoder_blocks = 2)))
+    z = randn(2, 30)
+    @test ReducedHamiltonianLoss(H)(dec, z, targets(dec, z)) < 1e-6
+end
+
+# λ weighs the gradient term only: with targets whose states are exact, the loss is λ times the
+# gradient error.
+function test_weight()
+    dec = decoder(NeuralNetwork(SymplecticAutoencoder(4, 2; n_encoder_blocks = 2, n_decoder_blocks = 2)))
+    z = randn(2, 30)
+    out = targets(dec, z)
+    out[5:6, :] .*= 1.1
+    @test ReducedHamiltonianLoss(H; λ = 3)(dec, z, out) ≈ 3 * ReducedHamiltonianLoss(H)(dec, z, out) rtol = 1e-6
+end
+
+# Reverse mode through the loss gives a finite gradient with respect to the parameters.
+function test_differentiable()
+    dec = decoder(NeuralNetwork(SymplecticAutoencoder(4, 2; n_encoder_blocks = 2, n_decoder_blocks = 2)))
+    teacher = decoder(NeuralNetwork(SymplecticAutoencoder(4, 2; n_encoder_blocks = 2, n_decoder_blocks = 2)))
+    z = randn(2, 30)
+    loss = ReducedHamiltonianLoss(H)
+    out = targets(teacher, z)
+    dp = Zygote.gradient(ps -> loss(dec.model, ps, z, out), params(dec))[1]
+    @test dp !== nothing
+end
+
+# Training a decoder on the targets of another one with the optimizer reduces the loss.
+function test_training()
+    Random.seed!(2026)
+    teacher = decoder(NeuralNetwork(SymplecticAutoencoder(4, 2; n_encoder_blocks = 2, n_decoder_blocks = 2)))
+    student = decoder(NeuralNetwork(SymplecticAutoencoder(4, 2; n_encoder_blocks = 2, n_decoder_blocks = 2)))
+    z = randn(2, 200)
+    out = targets(teacher, z)
+    loss = ReducedHamiltonianLoss(H)
+    before = loss(student, z, out)
+    o = Optimizer(Adam(), student)
+    o(student, DataLoader(z, out; suppress_info = true), Batch(50), 200, loss; show_progress = false)
+    @test loss(student, z, out) < before / 2
+end
+
+# Where every target gradient of the batch vanishes, the gradient term is the absolute error, not NaN.
+function test_zero_target_gradient()
+    dec = decoder(NeuralNetwork(SymplecticAutoencoder(4, 2; n_encoder_blocks = 2, n_decoder_blocks = 2)))
+    z = randn(2, 10)
+    out = vcat(dec(z), zeros(2, 10))
+    l = ReducedHamiltonianLoss(H)(dec, z, out)
+    @test isfinite(l)
+    @test l ≈ norm(reduced_gradient(dec, z)) rtol = 1e-4
+end
+
+test_gradient()
+test_zero_target_gradient()
+test_self_consistency()
+test_weight()
+test_differentiable()
+test_training()

test/runtests.jl

@@ -50,6 +50,7 @@ if "core" in GROUPS
     @safetestset "Check reduced model" include("reduced_system/reduced_system.jl")
     @safetestset "_norm keeps the element type of its argument" include("norm_eltype.jl")
     @safetestset "Symplectic Euler and variational midpoint losses" include("loss/training_method_losses.jl")
+    @safetestset "ReducedHamiltonianLoss" include("loss/reduced_hamiltonian_loss.jl")
     @safetestset "Test NetworkLoss + Optimizer" include("integration/losses_and_optimization.jl")
     @safetestset "ReducedLoss through the Optimizer functor" include("integration/reduced_loss_optimization.jl")
     @safetestset "Optimizer #2" include("optimizers/optimization_step.jl")