From 69e7af8f8e2b80a45f7c97a99152e7c9a257f080 Mon Sep 17 00:00:00 2001 From: benedict-96 Date: Tue, 6 Oct 2026 18:49:06 +0200 Subject: [PATCH 1/3] Add ReducedHamiltonianLoss: fit a decoder to states and reduced-Hamiltonian gradients MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit For latent points z with target states x* and target gradients g* of the reduced Hamiltonian, the loss is the relative state error plus λ times the relative error of ∇_z(H∘Ψ^dec). For a symplectic decoder this gradient is the reduced vector field up to 𝕁, so the loss fits the reduced dynamics, critical points and separatrices included, which a fit of the states alone does not control. The z-gradient is taken by central differences, so the parameter gradient needs no nested AD. Co-Authored-By: Claude Opus 5.5 (1M context) --- CHANGELOG.md | 1 + docs/src/reduced_order_modeling/losses.md | 3 + src/GeometricMachineLearning.jl | 3 +- src/loss/reduced_hamiltonian_loss.jl | 96 +++++++++++++++++++++++ test/loss/reduced_hamiltonian_loss.jl | 68 ++++++++++++++++ test/runtests.jl | 1 + 6 files changed, 171 insertions(+), 1 deletion(-) create mode 100644 src/loss/reduced_hamiltonian_loss.jl create mode 100644 test/loss/reduced_hamiltonian_loss.jl diff --git a/CHANGELOG.md b/CHANGELOG.md index 7a84549a8..10a530d7b 100644 --- a/CHANGELOG.md +++ b/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. 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 ForwardDiff, 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_` 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. diff --git a/docs/src/reduced_order_modeling/losses.md b/docs/src/reduced_order_modeling/losses.md index 145233b2a..3560f9a18 100644 --- a/docs/src/reduced_order_modeling/losses.md +++ b/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 ``` diff --git a/src/GeometricMachineLearning.jl b/src/GeometricMachineLearning.jl index bb474b69b..d98e46703 100644 --- a/src/GeometricMachineLearning.jl +++ b/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") diff --git a/src/loss/reduced_hamiltonian_loss.jl b/src/loss/reduced_hamiltonian_loss.jl new file mode 100644 index 000000000..3e40f1b7c --- /dev/null +++ b/src/loss/reduced_hamiltonian_loss.jl @@ -0,0 +1,96 @@ +@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}. +``` + +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 ForwardDiff: 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], ζ) 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 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 = norm(model(input, ps) - x_target) / norm(x_target) + gradient_error = norm(_reduced_hamiltonian_gradient(loss, model, ps, input) - g_target) / + norm(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 diff --git a/test/loss/reduced_hamiltonian_loss.jl b/test/loss/reduced_hamiltonian_loss.jl new file mode 100644 index 000000000..fc2510059 --- /dev/null +++ b/test/loss/reduced_hamiltonian_loss.jl @@ -0,0 +1,68 @@ +using GeometricMachineLearning +using GeometricMachineLearning: _reduced_hamiltonian_gradient, params +using ForwardDiff: gradient +using Zygote: Zygote +using Test +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, [gradient(ζ -> H(reshape(dec(ζ), 4, 1))[1], ζ) for ζ in eachcol(z)]) +targets(dec, z) = vcat(dec(z), reduced_gradient(dec, z)) + +# The central differences agree with forward-mode AD. +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) + dp = Zygote.gradient(ps -> loss(dec.model, ps, z, targets(teacher, z)), 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 + +test_gradient() +test_self_consistency() +test_weight() +test_differentiable() +test_training() diff --git a/test/runtests.jl b/test/runtests.jl index 9d1e773bc..5d3605bdc 100644 --- a/test/runtests.jl +++ b/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") From af0925fa6dec3fd53f4d1e07b3d3aeef3135cbc5 Mon Sep 17 00:00:00 2001 From: benedict-96 Date: Tue, 6 Oct 2026 21:49:12 +0200 Subject: [PATCH 2/3] ReducedHamiltonianLoss: absolute error where a target vanishes; tests without ForwardDiff A batch whose target gradients all vanish (every latent point a critical point of the target) made the relative gradient error 0/0; that term is now the absolute error there (from a CodeRabbit review comment), with a test. The tests took their reference gradients from ForwardDiff, which test/Project.toml does not list, so every test job failed; they use Zygote now. Co-Authored-By: Claude Opus 5.5 (1M context) --- CHANGELOG.md | 2 +- src/loss/reduced_hamiltonian_loss.jl | 13 +++++++++---- test/loss/reduced_hamiltonian_loss.jl | 20 ++++++++++++++++---- 3 files changed, 26 insertions(+), 9 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 10a530d7b..753b6bbdc 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -17,7 +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. 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 ForwardDiff, 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. +- **`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_` 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. diff --git a/src/loss/reduced_hamiltonian_loss.jl b/src/loss/reduced_hamiltonian_loss.jl index 3e40f1b7c..0bbcd1b44 100644 --- a/src/loss/reduced_hamiltonian_loss.jl +++ b/src/loss/reduced_hamiltonian_loss.jl @@ -7,9 +7,11 @@ The loss fits a [`Decoder`](@ref) ``\Psi^\mathrm{dec}:\mathbb{R}^{2n}\to\mathbb{ ```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}. + + \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. @@ -63,6 +65,10 @@ function ReducedHamiltonianLoss(hamiltonian; λ::Real = 1, h::Real = 1e-4) 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) @@ -82,9 +88,8 @@ function (loss::ReducedHamiltonianLoss)(model::Union{Chain, AbstractExplicitLaye n = size(input, 1) x_target = output[1:(end - n), :] g_target = output[(end - n + 1):end, :] - state_error = norm(model(input, ps) - x_target) / norm(x_target) - gradient_error = norm(_reduced_hamiltonian_gradient(loss, model, ps, input) - g_target) / - norm(g_target) + 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 diff --git a/test/loss/reduced_hamiltonian_loss.jl b/test/loss/reduced_hamiltonian_loss.jl index fc2510059..d1fd1698f 100644 --- a/test/loss/reduced_hamiltonian_loss.jl +++ b/test/loss/reduced_hamiltonian_loss.jl @@ -1,17 +1,17 @@ using GeometricMachineLearning using GeometricMachineLearning: _reduced_hamiltonian_gradient, params -using ForwardDiff: gradient 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, [gradient(ζ -> H(reshape(dec(ζ), 4, 1))[1], ζ) for ζ in eachcol(z)]) +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 forward-mode AD. +# 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) @@ -43,7 +43,8 @@ function test_differentiable() teacher = decoder(NeuralNetwork(SymplecticAutoencoder(4, 2; n_encoder_blocks = 2, n_decoder_blocks = 2))) z = randn(2, 30) loss = ReducedHamiltonianLoss(H) - dp = Zygote.gradient(ps -> loss(dec.model, ps, z, targets(teacher, z)), params(dec))[1] + out = targets(teacher, z) + dp = Zygote.gradient(ps -> loss(dec.model, ps, z, out), params(dec))[1] @test dp !== nothing end @@ -61,7 +62,18 @@ function test_training() @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() From 94cfd16b346011c00c6017bf92ae540cd202128e Mon Sep 17 00:00:00 2001 From: benedict-96 Date: Tue, 6 Oct 2026 22:11:30 +0200 Subject: [PATCH 3/3] ReducedHamiltonianLoss doctest: take the reference gradients from Zygote The test suite runs the doctests in test/Project.toml, which has Zygote but not ForwardDiff, so the doctest's `using ForwardDiff` failed every test job. The docs environment has both. Co-Authored-By: Claude Opus 5.5 (1M context) --- src/loss/reduced_hamiltonian_loss.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/loss/reduced_hamiltonian_loss.jl b/src/loss/reduced_hamiltonian_loss.jl index 0bbcd1b44..69d77cf5a 100644 --- a/src/loss/reduced_hamiltonian_loss.jl +++ b/src/loss/reduced_hamiltonian_loss.jl @@ -37,7 +37,7 @@ A decoder fitted to its own states and reduced gradients has a loss of the size ```jldoctest using GeometricMachineLearning -using ForwardDiff: gradient +using Zygote: gradient import Random Random.seed!(123) @@ -45,7 +45,7 @@ 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], ζ) for ζ in eachcol(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