GML part M4: one row-pivoted Cayley kernel at every sequence length, with the closed forms for N = 2 to 5 kept

CHANGELOG.md

@@ -73,7 +73,7 @@ breaking release).
 
 - **`assign_q_and_p` splits its input by slicing, and a loss that reads one half no longer raises.** It takes `x[1:N, …]` and `x[N+1:2N, …]` with `N::Integer`, where it took `N::Int`. Its six kernels, its rrule and `_vcat_qp` are gone. The old rrule built a cotangent of `N + 1` rows, `vcat(ZeroTangent(), p̄)`, and raised `DimensionMismatch` when a loss read only one of the two halves. The test is `test/kernels/assign_q_and_p.jl`.
 
-- **The kernel rrules accept a zero seed and a seed of the other float type.** Each remaining kernel rrule (`tensor_tensor_mul`, `tensor_transpose_tensor_mul`, `tensor_mat_skew_sym_assign`, `cpu_inverse`, `tensor_inverse2` to `tensor_inverse5`, and `_mul(::Adjoint, x)`) answers a zero seed (`ZeroTangent`, `NoTangent` or a thunk of one) with `NoTangent()` or `ZeroTangent()`, where it raised a `MethodError`. It converts a seed of the other float type to the element type of the primal: a Float64 seed on a Float32 primal raised a `MethodError`, and `cpu_inverse` returned a Float32 `dA` for a Float64 primal. `init_output` already did this. The style is the same in every rule: eager, with no `@thunk`. The four generated `tensor_inverseN` rrules call one hand-written `_tensor_inverse_pullback`, in `src/kernels/inverses/tensor_inverse_pullback.jl`, and `scripts/utilities/inverse_generator.jl` emits that call. The tests are in `test/kernels/kernel_pullbacks.jl`: `test_rrule` in Float64, Float32 against a Float64 finite-difference reference, zero seeds, and seeds of the other precision.
+- **The kernel rrules accept a zero seed and a seed of the other float type.** Each remaining kernel rrule (`tensor_tensor_mul`, `tensor_transpose_tensor_mul`, `tensor_mat_skew_sym_assign`, `tensor_inverse2` to `tensor_inverse5`, and `_mul(::Adjoint, x)`) answers a zero seed (`ZeroTangent`, `NoTangent` or a thunk of one) with `NoTangent()` or `ZeroTangent()`, where it raised a `MethodError`. It converts a seed of the other float type to the element type of the primal, where a Float64 seed on a Float32 primal raised a `MethodError`. `init_output` already did this. The style is the same in every rule: eager, with no `@thunk`. The four generated `tensor_inverseN` rrules call one hand-written `_tensor_inverse_pullback`, in `src/kernels/inverses/tensor_inverse_pullback.jl`, and `scripts/utilities/inverse_generator.jl` emits that call. The tests are in `test/kernels/kernel_pullbacks.jl`: `test_rrule` in Float64, Float32 against a Float64 finite-difference reference, zero seeds, and seeds of the other precision.
 
 - **A network with a Stiefel leaf trains on data of the other precision.** Training on data whose element type differs from the network's raised `AssertionError: eltype(A) == eltype(B)` in `mat_tensor_mul`. `_mul` checks no element type, because GEMM promotes. This closes known issue K7, which leaves `KNOWN_ISSUES.md`. The test is `test/integration/mixed_precision_training.jl`.
 
@@ -98,6 +98,16 @@ breaking release).
 
 ### Changed
 
+- **`VolumePreservingAttention` computes the Cayley transform with the kernel `tensor_cayley` at every sequence length, 0 included, and on device arrays; the batch inverses are deleted.** `tensor_cayley(C)` is internal and not exported. It computes `(I + Cₖ)⁻¹(I − Cₖ)` for every slice `k` by row-pivoted Gauss–Jordan elimination, one work item per slice. Its `rrule` computes `C̄ = −(I + Λ)ᵀ Λ̄ (I + Λ)ᵀ / 2` from `Λ` alone, with no second solve and no identity tensor.
+
+  `cpu_inverse`, its `rrule` and `cpu_tensor_cayley` are deleted, and so are the closed forms `tensor_cayley2` to `tensor_cayley5`, the symbolic inverses `tensor_inverse2` to `tensor_inverse5` with their pullback, the identity tensor `init_output`, and `scripts/utilities/inverse_generator.jl`. This closes known issues C18 and K16, which leave `KNOWN_ISSUES.md`.
+
+  New tests: `test/kernels/cayley.jl` checks the kernel's orthogonality and value, and its pullback against a central difference, for N = 2 to 8. `test/layers/volume_preserving_attention.jl` checks the layer's forward pass and gradients for N = 2 to 8, on the host and on a `JLArray`. `test/devices/jlarray.jl`, `test/devices/metal.jl` and `test/quality/jet.jl` reach the kernel, the last in `Float16`, `Float32` and `Float64`.
+
+  The closed forms went because the kernel is not meaningfully slower. Against them, on 1000 skew slices of size N = 2 to 5 on an Apple M-series CPU under load, the kernel's forward time lay between 0.83 and 1.42 of theirs and its backward time between 0.40 and 0.57. On Metal, in `Float32`, the forward ratio lay between 0.28 and 1.28 and the backward ratio between 0.44 and 0.99. These are two cold CPU runs and four warm Metal runs, and one Metal ratio moved by up to 4.6× at N = 2, so the figures are noisy; no ratio exceeded 2.
+
+  The orthogonality check in `test/kernels/cayley.jl` compares the kernel's `‖ΛₖᵀΛₖ − I‖∞` with that of the pivoted reference `(I + Cₖ) \ (I − Cₖ)`, for N = 2 to 8, both element types and two scales, on three input sets. The ratio of the kernel's error to the reference's is held to a limit of 4 on all three sets. On the two exactly skew sets the kernel's error must also lie within `10·n·eps(T)·(1 + ‖Cₖ‖₂²)`; on the `:skew_sym` set, where `C` is skew only to rounding, the ratio is the only check. The value rows compare each `Λₖ` with a `BigFloat` solve rounded to `T` on the same three sets.
+
 - `Pkg.test()` no longer runs the doctests. `test/quality/doctests.jl` is now the group `doctests`, which an empty `ARGS` does not run; `Pkg.test(test_args = ["doctests"])` runs it. In CI the Doctests job stays their runner, so the test matrix no longer runs them a second time.
 - **CI uploads coverage from the `Julia 1 - ubuntu-latest` job instead of `Julia min`, and a test job saves the Julia cache only when it succeeds.**
 - **The package needs Julia 1.12 or later** (`[compat] julia = "1.12"`, was `"1.11"`). On Julia 1.11.9, the gradient of an `LNNLoss` with respect to the network parameters stops with `StackOverflowError` in type inference ("Internal error: during type inference of … Encountered stack overflow") on macOS and Windows; the `Julia min` job on Linux passes. The overflow comes from SymbolicUtils 4.49.0, whose `rrule` for `SymbolicUtils.Code.create_array` returns an exactly typed tuple of tangents, one per array element; Zygote's reverse pass through the generated Hessian of `lagrangian_acceleration` is then inferred far deeper. With the 4.48.0 rule restored, or in a task with a 64 MB stack, the same gradient passes on 1.11.9; Julia 1.12.7 and 1.13.1 pass it on the root task. This closes #327: the `Julia min` jobs on macOS and Windows had failed in the *Lagrangian Neural Network* testset since 2026-10-01, and every testset after it was never run there.
@@ -108,7 +118,7 @@ breaking release).
 - **The optimizer reads `show_progress` as a `Bool`** instead of comparing it with `== true`. A non-`Bool` value now raises a `TypeError`: before, `nothing`, `0` or a string hid the progress bar with no error, and `1` or `1.0` showed it. Every call in the tests and the manual passes `false`, and every call in `scripts/` passes `true`.
 - **fatou reports no finding in `src/` or `test/`.** Only tests and one lint suppression changed; no behaviour of the package changed. `test/layers/multi_head_attention_stiefel_setup.jl` defined `check_setup(ps::NetworkParameters)` twice, and the later method silently replaced the `mapparameters` one, which never ran. The later method is deleted, and the test now asserts that the matrix method of `check_setup` runs once per leaf of the parameters, counted with `foldparameters`. The unused destructured values `Σ`, `Vt` and `Ũ₂` in `test/optimizers/psd_optim.jl` and `test/optimizers/svd_optim.jl` are now `_`; those files check `Ũ₁` and `Ũ₃` and never checked `Ũ₂`. The unreachable `new{…}(…)` after the `error` in the `GeneralizedHamiltonianArchitecture` constructor stays, because it records the field order for the future constructor, under `# fatou-ignore unreachable-code`.
 - **Three test files left the top level of `test/`**, because a test file sits there only where it mirrors `src/<name>.jl`. `test/norm_eltype.jl` tests `_norm` of `src/utils.jl` and is now `test/utils.jl`; `test/map_to_cpu_tests.jl` tests `src/map_to_cpu.jl` and is now `test/map_to_cpu.jl`; `test/documented_defaults.jl` tests the keyword defaults of `src/layers/symplectic_attention.jl`, `src/architectures/symplectic_transformer.jl` and `src/architectures/standard_transformer_integrator.jl`, whose deepest common directory is `src/` itself, and is now `test/integration/documented_defaults.jl`. Only the paths in `test/runtests.jl` changed; the labels, groups and order are the same.
-- **The matrix and tensor products are one internal product, `_mul`.** `src/kernels/mul.jl` defines `_mul(W, x)` as `W * x` for a vector or matrix `x`, so a structured weight is multiplied by its own `*`, and as one GEMM over `reshape(x, size(x, 1), :)` for a 3-tensor `x`; `_mul(x3, W)` as one GEMM between two `permutedims` of `x3`; and `_mul(x3, y3) = tensor_tensor_mul`. It replaces `mat_tensor_mul` and `mat_tensor_mul!`, `tensor_mat_mul` and `tensor_mat_mul!`, `symmetric_mat_mul`, `symmetric_mat_right_mul`, `lo_mat_mul`, `up_mat_mul`, `skew_mat_mul`, `tensor_transpose_mat_mul`, `custom_mat_mul`, `_custom_mul`, `_matrix_cotangent` and the Stiefel methods of `mat_tensor_mul`, together with their kernels and rrules. `custom_vec_mul` is `.*` now, and `vec_tensor_mul` and `tensor_scalar_product` are gone; a broadcast, and Zygote's rule for it, replaces them. A `SymmetricMatrix`, `SkewSymMatrix`, `StrictlyLowerTriangular` or `StrictlyUpperTriangular` weight is made dense once per call, as `W * one(W)` on its own backend, in a 3-tensor left product and in every right product; a vector or matrix `x` in a left product takes the weight's own `*` instead. A structured matrix as the right argument of a left product is made dense the same way. Its gradient reaches the storage through GeometricOptimizers' `ProjectTo` and NeuralNetworkParameters' `storage_gradient`. A short rrule on `_mul(::Adjoint, x)` keeps the cotangent of a `Matrix` parameter that is used as `W'` a `Matrix`, at every rank of `x` and for a Stiefel leaf. `ManifoldLayer`, `ResNetLayer` and `VolumePreservingFeedForwardLayer` each have one `x::AbstractArray` method in place of a matrix method and a tensor method. `tensor_transpose`, with its kernel and rrule, is gone: callers write `permutedims(x, (2, 1, 3))`. `cpu_tensor_cayley` loses its shadowed duplicate in `cpu_inverse.jl`, and the four kernel files that survive lose the file-top strings that described the removed kernels. No exported name changes, so this is not a breaking change.
+- **The matrix and tensor products are one internal product, `_mul`.** `src/kernels/mul.jl` defines `_mul(W, x)` as `W * x` for a vector or matrix `x`, so a structured weight is multiplied by its own `*`, and as one GEMM over `reshape(x, size(x, 1), :)` for a 3-tensor `x`; `_mul(x3, W)` as one GEMM between two `permutedims` of `x3`; and `_mul(x3, y3) = tensor_tensor_mul`. It replaces `mat_tensor_mul` and `mat_tensor_mul!`, `tensor_mat_mul` and `tensor_mat_mul!`, `symmetric_mat_mul`, `symmetric_mat_right_mul`, `lo_mat_mul`, `up_mat_mul`, `skew_mat_mul`, `tensor_transpose_mat_mul`, `custom_mat_mul`, `_custom_mul`, `_matrix_cotangent` and the Stiefel methods of `mat_tensor_mul`, together with their kernels and rrules. `custom_vec_mul` is `.*` now, and `vec_tensor_mul` and `tensor_scalar_product` are gone; a broadcast, and Zygote's rule for it, replaces them. A `SymmetricMatrix`, `SkewSymMatrix`, `StrictlyLowerTriangular` or `StrictlyUpperTriangular` weight is made dense once per call, as `W * one(W)` on its own backend, in a 3-tensor left product and in every right product; a vector or matrix `x` in a left product takes the weight's own `*` instead. A structured matrix as the right argument of a left product is made dense the same way. Its gradient reaches the storage through GeometricOptimizers' `ProjectTo` and NeuralNetworkParameters' `storage_gradient`. A short rrule on `_mul(::Adjoint, x)` keeps the cotangent of a `Matrix` parameter that is used as `W'` a `Matrix`, at every rank of `x` and for a Stiefel leaf. `ManifoldLayer`, `ResNetLayer` and `VolumePreservingFeedForwardLayer` each have one `x::AbstractArray` method in place of a matrix method and a tensor method. `tensor_transpose`, with its kernel and rrule, is gone: callers write `permutedims(x, (2, 1, 3))`. The four kernel files that survive lose the file-top strings that described the removed kernels. No exported name changes, so this is not a breaking change.
 
   Measured with `scripts/utilities/product_timings.jl` at 10×10×1000, in µs, before → after. The script's `x` is a 3-tensor, so this table does not cover a vector or matrix `x`. The before column is the same script run at `2a63f26a` with the old names, `mat_tensor_mul` for the left product and `tensor_mat_mul` for the right, with base and branch processes alternated. BLAS ran on one thread. Each variant ran in five fresh processes, each doing one warm-up and taking the minimum of 100 calls, and the table gives the median of the five. The backward column is `Zygote.gradient` through a `NetworkParameters`. The machine was an Apple M-series with Julia 1.13.1, GeometricOptimizers 0.9.0, NeuralNetworkParameters 0.5.2 and Zygote 0.7.13, shared with other work, at a load average between 8 and 37.
 

KNOWN_ISSUES.md

@@ -89,16 +89,6 @@ What is known to be wrong in GeometricMachineLearning and is not fixed.
   deliberately, and the resolution is upstream declaring them `public`. Each check's reason is
   written beside it in the test file.
 
-### C18 · `src/kernels/exponentials/tensor_exponential.jl` no longer defines `tensor_exponential`; it holds the identity tensor and its `rrule`, used by `tensor_cayley.jl` and `cpu_inverse.jl`.
-
-- **location:** `src/kernels/exponentials/tensor_exponential.jl`
-- **kind:** defect
-- **found:** 2026-09-19
-- **evidence:**
-  The file name, and the `exponentials/` directory around it, are wider than what they hold. A
-  comment in the file itself cites "issue C18" by number, so this entry exists to keep that reference
-  valid. Renaming both the file and the directory is the fix.
-
 ### C19 · `_make_block_for_initialization` has two methods
 
 - **location:** `src/architectures/symplectic_transformer.jl`
@@ -208,32 +198,6 @@ What is known to be wrong in GeometricMachineLearning and is not fixed.
   belongs with the data signatures of the losses as a whole: either each loss accepts only the data
   its vector field can take, or `(q, p)` data is converted before it reaches the field.
 
-### K16 · `scripts/utilities/inverse_generator.jl` reproduces the wrappers of the committed inverse files byte for byte, but not their kernels.
-
-- **location:** `scripts/utilities/inverse_generator.jl:1`
-- **kind:** found late
-- **found:** 2026-10-15, by the builder of part M3; the same output comes from the generator of `origin/main` at `2a63f26a`
-- **evidence:**
-  Formatted with `style = "sciml"`, the generator writes the kernel of each of the four
-  `inverse_NxN.jl` files, the part before `function tensor_inverseN(`, with the factors of some
-  products in another order than the committed file, from line 6 or 7 on:
-
-  ```
-  <                 t2 = (*)((*)(t1, (getindex)(A, 2, 1, k)), (getindex)(A, 1, 2, k))
-  ---
-  >                 t2 = (*)((*)(t1, (getindex)(A, 1, 2, k)), (getindex)(A, 2, 1, k))
-  ```
-
-  The order depends on the Julia version, not on the package versions. With Symbolics 7.39.2,
-  SymbolicUtils 4.46.6 and JuliaFormatter 2.14.0 fixed, Julia 1.11.9, 1.12.6, 1.12.7 and 1.13.1
-  each give a different order, and none of them reproduces all four committed kernels; 1.12.6
-  reproduces the 2×2 one. On Julia 1.13.1, SymbolicUtils 4.46.6 and 4.51.0 give identical output,
-  and two runs on one Julia version give identical output. The Julia build that made the committed
-  kernels on 2026-09-14 is no longer installed, so that run cannot be repeated.
-
-  The rest of each file, from `function tensor_inverseN(` to the end, is byte for byte what the
-  generator's `wrapper_source(n)` writes, and the script's header claims no more than that.
-
 ### K17 · `assign_q_and_p` and `PoissonTensor * v` drop the rows from `2N + 1` on without an error.
 
 - **location:** `src/kernels/assign_q_and_p.jl`
@@ -247,6 +211,33 @@ What is known to be wrong in GeometricMachineLearning and is not fixed.
   Closing it means an error for `size(x, 1) != 2N`. That changes the public behaviour of
   `assign_q_and_p` and of `PoissonTensor * v`, so it is a change of its own and not a cleanup.
 
+### K19 · In `Float32` the Cayley kernel returns non-finite slices for a skew `C` of norm 1e8, where no method is accurate.
+
+- **location:** `src/kernels/cayley.jl`
+- **kind:** found late
+- **found:** 2026-10-17, by the critics of part M4
+- **evidence:**
+  On 200 skew `Float32` slices of `N = 5` whose largest `‖Cₖ‖₂` is 6.55e8, `tensor_cayley` returns
+  non-finite values in 23 of them and raises no error. The scale is far outside check 3's `{1, 1e2}`,
+  and no method is accurate there: the dense solve `(I + Cₖ) \ (I − Cₖ)` is finite on all 200 but has an
+  orthogonality error `‖ΛₖᵀΛₖ − I‖∞` of 7.3e4 at worst, which is no answer at that scale. The cause
+  is cancellation in the last pivot, not overflow: in all 23 non-finite slices the last `Float32`
+  pivot is exactly `0`. The same elimination in `Float64` gives a last pivot of magnitude 1.7 to 15
+  on those slices, and no intermediate entry above 6.8e8. In `Float32` the rounding of the earlier
+  steps, about `eps(Float32) · ‖Cₖ‖₂ ≈ 78`, is larger than that pivot; for odd `N` the skew `Cₖ` is
+  singular, so the pivot of `I + Cₖ` is small at that scale. The kernel's error on its 177 finite
+  slices is at most 284.
+
+  ```julia
+  using GeometricMachineLearning: tensor_cayley
+  import Random
+  Random.seed!(3)
+  R = 1.0f8 * randn(Float32, 5, 5, 200)
+  C = R - permutedims(R, (2, 1, 3))
+  nonfinite(L) = count(k -> !all(isfinite, L[:, :, k]), 1:200)
+  nonfinite(tensor_cayley(C))    # 23 of 200
+  ```
+
 ## Upstream
 
 ### K2 · Revise prints EMFILE errors in the test log

docs/src/layers/attention_layer.md

@@ -258,6 +258,7 @@ Attention was used before [the transformer](@ref "Standard Transformer") was int
 
 ```@docs
 GeometricMachineLearning.tensor_mat_skew_sym_assign
+GeometricMachineLearning.tensor_cayley
 VolumePreservingAttention
 ```