diff --git a/CHANGELOG.md b/CHANGELOG.md
index 7a84549a8..124844eb5 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -17,6 +17,7 @@ breaking release).
 
 ### Added
 
+- **`LatentVectorFieldLoss` trains an autoencoder on states and on the full vector field there.** For states ``x`` with the vector field ``X_H(x)`` of a Hamiltonian `H` that the user supplies as a function of a matrix of states, and ``\xi=\Psi^\mathrm{enc}(x)``, the loss is ``\lVert\Psi^\mathrm{dec}(\xi)-x\rVert/\lVert x\rVert + \lambda\lVert\mathbb{J}_{2n}\nabla_\xi(H\circ\Psi^\mathrm{dec})(\xi)-\nabla\Psi^\mathrm{enc}(x)X_H(x)\rVert/\lVert\nabla\Psi^\mathrm{enc}(x)X_H(x)\rVert``: the reduced vector field is compared in the latent space with the full one pushed forward by the encoder, so encoder and decoder are trained on the dynamics without latent targets. Both derivatives are central differences (step `h`), so the parameter gradient needs no nested automatic differentiation. It is built from the `AutoEncoder` architecture, which fixes where the encoder ends in the chain, and accepts matrices and the 3-tensors of a `DataLoader` built from an input and an output matrix. With `reconstruction = false` the reconstruction term is left out and only the reduced dynamics are fitted. With `energy = true` the loss also has the relative error of ``H\circ\Psi^\mathrm{dec}\circ\Psi^\mathrm{enc}`` against ``H``, which keeps its weight where the vector field vanishes; with `reconstruction = false` it replaces the reconstruction term. `test/loss/latent_vector_field_loss.jl` checks the loss against the same expression with Zygote's derivatives, the loss without reconstruction, the energy term, the weight `λ`, the Zygote gradient, and that 200 epochs of training halve it.
 - **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.
diff --git a/docs/src/reduced_order_modeling/losses.md b/docs/src/reduced_order_modeling/losses.md
index 145233b2a..73793da3c 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)).
 
+[`LatentVectorFieldLoss`](@ref) trains encoder and decoder on states of the full system and on its vector field there: besides the reconstruction it compares the reduced vector field ``\mathbb{J}_{2n}\nabla_\xi(H\circ\Psi^\mathrm{dec})`` at ``\xi=\Psi^\mathrm{enc}(x)`` with the full vector field pushed forward by the encoder.
+
 ## Library Functions
 
 ```@docs
 TransformerLoss
 AutoEncoderLoss
 ReducedLoss
+LatentVectorFieldLoss
 projection_error
 reduction_error
 ```
diff --git a/src/GeometricMachineLearning.jl b/src/GeometricMachineLearning.jl
index bb474b69b..93fe9587e 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, LatentVectorFieldLoss
 
 #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/latent_vector_field_loss.jl")
 
 export AbstractPullback, ZygotePullback, SymbolicPullback
 include("pullbacks/zygote_pullback.jl")
diff --git a/src/loss/latent_vector_field_loss.jl b/src/loss/latent_vector_field_loss.jl
new file mode 100644
index 000000000..869d3fe21
--- /dev/null
+++ b/src/loss/latent_vector_field_loss.jl
@@ -0,0 +1,131 @@
+@doc raw"""
+    LatentVectorFieldLoss(arch, hamiltonian; λ = 1, h = 1e-4, reconstruction = true, energy = false)
+
+Make an instance of `LatentVectorFieldLoss` for an [`AutoEncoder`](@ref) architecture `arch` and a Hamiltonian ``H`` on the full space.
+
+The loss trains encoder and decoder together on states ``x`` of the full system and on its vector field ``X_H(x) = \mathbb{J}_{2N}\nabla{}H(x)`` there. With ``\xi = \Psi^\mathrm{enc}(x)`` it is
+
+```math
+L = \frac{\lVert\Psi^\mathrm{dec}(\xi) - x\rVert}{\lVert x\rVert}
+  + \lambda\,\frac{\lVert\mathbb{J}_{2n}\nabla_\xi(H\circ\Psi^\mathrm{dec})(\xi) - \nabla\Psi^\mathrm{enc}(x)X_H(x)\rVert}{\lVert\nabla\Psi^\mathrm{enc}(x)X_H(x)\rVert},
+```
+
+the reconstruction error plus ``\lambda`` times the error of the *reduced vector field* against the full vector field pushed forward by the encoder, both relative and with the norms taken over the batch.
+
+!!! warning "Reconstruction and dynamics can contradict each other"
+    The reconstruction term asks that ``\Psi^\mathrm{dec}(\Psi^\mathrm{enc}(x)) = x``; the vector-field and energy terms ask that the reduced Hamiltonian system reproduce the full system's trajectories and energies. Where no symplectic chart can do both (for instance where the data contain a separatrix that a single chart cannot reconstruct), a network trained on both must fail one of them, and the loss decides which. All terms are relative errors summed with unit weights (and `λ`), so nothing depends on units; but that is a convention for where the network gives way, not a resolution of the conflict.
+
+!!! warning "The vector field has to be the one the data follow"
+    When the states lie on a submanifold, as the states of a constrained system do, the target is the system's vector field there, which is tangent to it, and not ``\mathbb{J}_{2N}\nabla{}H`` of a Hamiltonian that merely restricts to the right energy on it. The push-forward of a field that leaves the submanifold involves the derivatives of the encoder off the data, which do not describe how encoded trajectories move. The second term compares the two vector fields in the latent space: it vanishes when the encoder maps the trajectories of the full system onto those of the reduced one, which for a symplectic decoder are the level sets of ``H\circ\Psi^\mathrm{dec}``. No target decoder and no latent targets are needed, only states and the full vector field. Where the pushed-forward field vanishes on the whole batch, the absolute error replaces the relative one.
+
+Both derivatives are central differences with step `h`: ``\nabla_\xi(H\circ\Psi^\mathrm{dec})`` from ``2\cdot2n`` decoder evaluations, and ``\nabla\Psi^\mathrm{enc}(x)X_H(x) \approx (\Psi^\mathrm{enc}(x + hX_H) - \Psi^\mathrm{enc}(x - hX_H))/(2h)`` from two encoder evaluations. 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
+
+- `arch`: the [`AutoEncoder`](@ref) whose network is trained; it fixes where the encoder ends in the chain.
+- `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 vector-field term.
+- `h = 1e-4`: the step of the central differences.
+- `energy = false`: with `true` the loss also has the relative energy error
+  ``\lVert H(\Psi^\mathrm{dec}(\xi)) - H(x)\rVert/\lVert H(x)\rVert``: the reduced Hamiltonian
+  ``H\circ\Psi^\mathrm{dec}`` has to take the energy of the state it encodes. Unlike the
+  reconstruction it does not tie ``\Psi^\mathrm{dec}(\xi)`` to ``x``, and unlike the vector-field
+  term it keeps its weight where the vector field vanishes. Together with `reconstruction = false`
+  it replaces the reconstruction term.
+- `reconstruction = true`: with `false` the reconstruction term is left out and the loss is the vector-field term alone, ``\lambda`` times its relative error. Only the reduced dynamics are then fitted; nothing ties ``\Psi^\mathrm{dec}(\xi)`` to ``x`` except through ``H\circ\Psi^\mathrm{dec}``.
+
+# Functor
+
+```julia
+loss(model, ps, input, output)
+loss(nn, input, output)
+```
+
+`input` holds the states ``x``, one per column (``2N`` rows), and `output` the vector field ``X_H(x)`` at them (``2N`` rows). The model is the whole autoencoder.
+
+# Examples
+
+The loss is zero when the vector-field targets are what the network itself predicts, here for the linear Hamiltonian ``H(x) = c\cdot x`` with a constant vector field:
+
+```jldoctest
+using GeometricMachineLearning
+import Random
+Random.seed!(123)
+
+arch = SymplecticAutoencoder(4, 2; n_encoder_blocks = 2, n_decoder_blocks = 2)
+nn = NeuralNetwork(arch)
+loss = LatentVectorFieldLoss(arch, X -> vec(sum(X; dims = 1)))
+x = randn(4, 10)
+v = ones(4, 10)
+loss(nn, x, v) > 0
+
+# output
+
+true
+```
+"""
+struct LatentVectorFieldLoss{N, HT, T} <: NetworkLoss
+    hamiltonian::HT
+    λ::T
+    h::T
+    reconstruction::Bool
+    energy::Bool
+    function LatentVectorFieldLoss{N}(hamiltonian::HT, λ::T, h::T, reconstruction::Bool = true, energy::Bool = false) where {N, HT, T}
+        new{N, HT, T}(hamiltonian, λ, h, reconstruction, energy)
+    end
+end
+
+function LatentVectorFieldLoss(arch::AutoEncoder, hamiltonian; λ::Real = 1, h::Real = 1e-4, reconstruction::Bool = true, energy::Bool = false)
+    T = float(promote_type(typeof(λ), typeof(h)))
+    N = length(encoder_model(arch).layers)
+    LatentVectorFieldLoss{N}(hamiltonian, T(λ), T(h), reconstruction, energy)
+end
+
+n_encoder_layers(::LatentVectorFieldLoss{N}) where {N} = N
+
+# The relative error, and the absolute one where the target vanishes.
+_relative_vf_error(prediction, target) = (n = norm(target); iszero(n) ? norm(prediction) : norm(prediction - target) / n)
+
+# The encoder and the decoder of the chain, each with its parameters.
+function _split_autoencoder(loss::LatentVectorFieldLoss{N}, model::Chain, ps) where {N}
+    vals = values(ps)
+    M = length(model.layers)
+    enc = Chain(ntuple(i -> model.layers[i], Val(N))...)
+    dec = Chain(ntuple(i -> model.layers[N + i], Val(M - N))...)
+    enc, ntuple(i -> vals[i], Val(N)), dec, ntuple(i -> vals[N + i], Val(M - N))
+end
+
+# ∇_ξ(H∘Ψ^dec) for each column of ξ, by central differences.
+function _latent_hamiltonian_gradient(loss::LatentVectorFieldLoss, dec, ps_dec, ξ::AbstractMatrix)
+    n = size(ξ, 1)
+    h = convert(eltype(ξ), 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(dec(ξ .+ δ, ps_dec))
+        Hm = loss.hamiltonian(dec(ξ .- δ, ps_dec))
+        reshape((Hp .- Hm) ./ (2h), 1, :)
+    end
+    reduce(vcat, rows)
+end
+
+# 𝕁_{2n} g = (g_p, -g_q) for g = (g_q, g_p).
+_poisson_latent(g::AbstractMatrix) = (n = size(g, 1) ÷ 2; vcat(g[(n + 1):end, :], -g[1:n, :]))
+
+function (loss::LatentVectorFieldLoss)(model::Chain, ps, input::AbstractMatrix, output::AbstractMatrix)
+    enc, ps_enc, dec, ps_dec = _split_autoencoder(loss, model, ps)
+    h = convert(eltype(input), loss.h)
+    ξ = enc(input, ps_enc)
+    decoded = dec(ξ, ps_dec)
+    reconstruction_error = loss.reconstruction ? _relative_vf_error(decoded, input) : zero(eltype(input))
+    energy_error = loss.energy ? _relative_vf_error(loss.hamiltonian(decoded), loss.hamiltonian(input)) : zero(eltype(input))
+    pushed_forward = (enc(input .+ h .* output, ps_enc) .- enc(input .- h .* output, ps_enc)) ./ (2h)
+    reduced = _poisson_latent(_latent_hamiltonian_gradient(loss, dec, ps_dec, ξ))
+    reconstruction_error + energy_error + loss.λ * _relative_vf_error(reduced, pushed_forward)
+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 states either way.
+function (loss::LatentVectorFieldLoss)(model::Chain, 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/latent_vector_field_loss.jl b/test/loss/latent_vector_field_loss.jl
new file mode 100644
index 000000000..9cea5532d
--- /dev/null
+++ b/test/loss/latent_vector_field_loss.jl
@@ -0,0 +1,78 @@
+using GeometricMachineLearning
+using GeometricMachineLearning: 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 ℝ⁴
+XH(X) = vcat(X[3:4, :], vcat(zeros(1, size(X, 2)), -ones(1, size(X, 2))))   # 𝕁∇H = (p, -∇_q H)
+arch = SymplecticAutoencoder(4, 2; n_encoder_blocks = 2, n_decoder_blocks = 2)
+
+# The loss agrees with the same expression evaluated with Zygote's derivatives.
+function test_against_autodiff()
+    nn = NeuralNetwork(arch)
+    enc, dec = encoder(nn), decoder(nn)
+    x = randn(4, 30); v = XH(x)
+    ξ = enc(x)
+    pushed = reduce(hcat, [Zygote.jacobian(enc, x[:, i])[1] * v[:, i] for i in axes(x, 2)])
+    g = reduce(hcat, [Zygote.gradient(ζ -> H(reshape(dec(ζ), 4, 1))[1], ξ[:, i])[1] for i in axes(x, 2)])
+    reduced = vcat(g[2:2, :], -g[1:1, :])
+    expected = norm(dec(ξ) - x) / norm(x) + norm(reduced - pushed) / norm(pushed)
+    @test LatentVectorFieldLoss(arch, H)(nn, x, v) ≈ expected rtol = 1e-6
+end
+
+# λ weighs the vector-field term only.
+function test_weight()
+    nn = NeuralNetwork(arch)
+    x = randn(4, 30); v = XH(x)
+    rec = LatentVectorFieldLoss(arch, H; λ = 0)(nn, x, v)
+    @test LatentVectorFieldLoss(arch, H; λ = 3)(nn, x, v) - rec ≈ 3 * (LatentVectorFieldLoss(arch, H)(nn, x, v) - rec) rtol = 1e-8
+end
+
+# Reverse mode through the loss gives a finite gradient for encoder and decoder parameters.
+function test_differentiable()
+    nn = NeuralNetwork(arch)
+    x = randn(4, 30); v = XH(x)
+    loss = LatentVectorFieldLoss(arch, H)
+    dp = Zygote.gradient(ps -> loss(nn.model, ps, x, v), params(nn))[1]
+    @test dp !== nothing
+end
+
+# Training with the optimizer reduces the loss.
+function test_training()
+    Random.seed!(2026)
+    nn = NeuralNetwork(arch)
+    x = randn(4, 200); v = XH(x)
+    loss = LatentVectorFieldLoss(arch, H)
+    before = loss(nn, x, v)
+    o = Optimizer(Adam(), nn)
+    o(nn, DataLoader(x, v; suppress_info = true), Batch(50), 200, loss; show_progress = false)
+    @test loss(nn, x, v) < before / 2
+end
+
+# Without the reconstruction term the loss is the vector-field term alone.
+function test_without_reconstruction()
+    nn = NeuralNetwork(arch)
+    x = randn(4, 30); v = XH(x)
+    rec = LatentVectorFieldLoss(arch, H; λ = 0)(nn, x, v)
+    @test LatentVectorFieldLoss(arch, H; reconstruction = false)(nn, x, v) ≈ LatentVectorFieldLoss(arch, H)(nn, x, v) - rec rtol = 1e-10
+end
+
+# The energy term is the relative error of H∘Ψ^dec∘Ψ^enc against H.
+function test_energy()
+    nn = NeuralNetwork(arch)
+    x = randn(4, 30); v = XH(x)
+    vf = LatentVectorFieldLoss(arch, H; reconstruction = false)(nn, x, v)
+    Hd = H(decoder(nn)(encoder(nn)(x)))
+    @test LatentVectorFieldLoss(arch, H; reconstruction = false, energy = true)(nn, x, v) ≈ vf + norm(Hd - H(x)) / norm(H(x)) rtol = 1e-8
+end
+
+test_against_autodiff()
+test_without_reconstruction()
+test_energy()
+test_weight()
+test_differentiable()
+test_training()
diff --git a/test/runtests.jl b/test/runtests.jl
index 9d1e773bc..c19087972 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 "LatentVectorFieldLoss" include("loss/latent_vector_field_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")