Deflate a singular assembly instead of rejecting it

CHANGELOG.md

@@ -5,6 +5,80 @@ All notable changes to SimpleSplines.jl are documented here.
 The format follows [Keep a Changelog](https://keepachangelog.com/en/1.1.0/), and the project
 adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
+## [Unreleased] — targeting 0.2.0
+
+### New Features
+
+A `kernel` keyword on `mass_operator`, `CirculantMass` and `FactorizedMass` handles singular
+assemblies on periodic bases. The three `MassOperator` representations carry not only the mass
+matrix but also other assemblies from a `SplineQuadrature`, and a stiffness matrix on a periodic
+basis is singular: the constants lie in its kernel because `-φ'' = ρ` says nothing about them.
+
+`kernel = :reject` keeps the previous behaviour — a singular matrix raises an error.
+`kernel = :project` states that the kernel is the constants and asks for the mean-free
+solution. Both periodic representations implement it in the way their own structure allows, and
+both verify the assertion rather than trusting it: `CirculantMass` checks that the constant
+Fourier mode is the one that vanishes and that no other does, and `FactorizedMass` checks that
+`M𝟙 ≈ 0`.
+
+For a caller solving a periodic equation, this replaces the textbook rank-one shift
+`S + 𝟙𝟙ᵀ/N`, which removes the singularity by making the matrix structurally full. At degree 4
+on 1024 cells the shift takes a periodic stiffness matrix `S` from 9 216 stored entries
+(155 816 B) to 1 048 576 (`S .+ inv(N)`, which stays a `SparseMatrixCSC` and so carries a row
+index beside every entry: 16 785 576 B; as the dense `Matrix(S) .+ inv(N)` it is 8 388 656 B),
+and construction from `O(N)` to `O(N²)`. The deflation leaves the matrix untouched. Measured
+over 32 to 1024 cells: the deflated solve and the shifted one agree to 1e-15 relative, and the
+mean of the result is zero to 1e-15.
+
+A bounded basis raises rather than accepting `kernel = :project`; the banded representation
+carries no deflation.
+
+### Bug Fixes
+
+The circulant invertibility check used an absolute tolerance, and so accepted a singular
+matrix. The zero eigenvalue of a periodic stiffness matrix comes out of the transform at about
+1e-14 rather than at zero, so `abs(λ) > eps(T)` passed it, the solve then divided by it, and
+what came back was amplified noise rather than an error. The tolerance is now relative —
+`N * eps(T)` times the largest eigenvalue, where `pinv` puts the same cutoff — so a matrix
+singular in this sense raises.
+
+Every tolerance in the mass operators is now relative and depends on the element type. Three
+of them were absolute, and each accepted something it should have refused:
+
+  - The constant-kernel check of `kernel = :project` could not fall below `sqrt(eps(T))`,
+    which made it absolute for every matrix whose entries are smaller than that. An
+    *invertible* mass matrix scaled below 1e-8 was accepted as having the constants in its
+    kernel, and the solve then returned a vector with relative residual of order one rather
+    than raising. The bound is now a decade above `N * eps(T)` times `‖M‖∞`, with no floor:
+    the residual of a genuine singular assembly saturates half that scale in the worst case,
+    and the smallest violation an invertible one produces sits 1e4 to 1e12 above it.
+  - The circulance check used `atol = 1e-10`, which is below `Float32` rounding: **no
+    `Float32` assembly was circulant at all**, so `CirculantMass` could not be built in that
+    element type and every periodic uniform `Float32` quadrature raised. The tolerance is now
+    `rtol`, relative to the largest entry of the first column and defaulting to
+    `sqrt(eps(T))`.
+  - `FactorizedMass` with `kernel = :project` took CHOLMOD's `issuccess` as proof that the
+    minor is positive definite. CHOLMOD calls a positive *semi*definite factorisation a
+    success, so a matrix whose kernel is one dimension larger than the constants was accepted,
+    and the solve returned amplified noise — for a 4×4 example, entries of 8.5e15 at a
+    residual of 2.0. The pivots of the factor are now read: the smallest must exceed
+    `N * eps(T)` of the largest, where a genuine assembly holds the ratio above 0.19.
+
+### Breaking Changes
+
+`FactorizedMass` gained a fourth type parameter, which carries the kernel mode
+(`{T, MT, FT, K}`). Code that spells the type out as `FactorizedMass{T, MT, FT}` in a signature
+needs the extra parameter. The constructor is unchanged.
+
+A `CirculantMass` solve multiplies by precomputed reciprocal eigenvalues instead of dividing by
+the eigenvalues themselves. It moves `N ÷ 2 + 1` divisions out of every solve and into the
+construction, and results can differ in the last ULP.
+
+`CirculantMass`'s `atol` keyword is now `rtol`, and it is relative rather than absolute: the
+tolerance is multiplied by the largest entry of the first column, and it defaults to
+`sqrt(eps(T))` rather than to `1e-10`. A caller that passed `atol` needs the new name, and a
+caller that passed a figure chosen against entries of a particular size needs a new figure.
+
 ## [0.1.0] — 2026-09-07
 
 The first release with any implementation in it. Before this the package was a

docs/src/usage/quadrature.md

@@ -232,7 +232,7 @@ that are not translates of anything.
 
 | basis | mesh | representation | a solve is |
 |:--|:--|:--|:--|
-| [`PeriodicBSplineBasis`](@ref) | [`UniformMesh`](@ref) | [`CirculantMass`](@ref) | two planned transforms and a pointwise division, ``O(N\log N)``, no allocation |
+| [`PeriodicBSplineBasis`](@ref) | [`UniformMesh`](@ref) | [`CirculantMass`](@ref) | two planned transforms and a pointwise multiplication, ``O(N\log N)``, no allocation |
 | [`PeriodicBSplineBasis`](@ref) | any other | [`FactorizedMass`](@ref) | a sparse Cholesky; the only representation whose solve allocates |
 | bounded — clamped or recombined | any | [`BandedMass`](@ref) | a banded Cholesky, ``O(Np)``, no allocation |
 
@@ -246,9 +246,11 @@ for (name, bb) in (("periodic, uniform", BSplineBasis(UniformMesh(16, 0 .. 2π),
 end
 ```
 
-The circulant construction *verifies* circulance rather than assuming it, to `atol = 1e-10` by
-default. A matrix that is banded but not circulant would still produce plausible numbers
-through the transform, and the failure would surface much later as a wrong conservation law.
+The circulant construction *verifies* circulance rather than assuming it, to `rtol =
+sqrt(eps(T))` by default — *relative* to the largest entry of the first column, so that the
+verdict survives a rescaling of the assembly and holds in every element type. A matrix that is
+banded but not circulant would still produce plausible numbers through the transform, and the
+failure would surface much later as a wrong conservation law.
 
 ```@example asm
 qg = SplineQuadrature(BSplineBasis(GradedMesh(16, 0 .. 2π), 3, Periodic()))
@@ -259,6 +261,56 @@ catch err
 end
 ```
 
+### A singular assembly
+
+The mass matrix is positive definite, but the same three representations carry the other
+assemblies of a quadrature, and [`stiffness_matrix`](@ref) on a periodic basis is *singular*:
+``-\phi'' = \rho`` says nothing about the constants, and a periodic basis represents them.
+
+```@example asm
+bp = BSplineBasis(UniformMesh(64, 0 .. 2π), 3, Periodic())
+qp = SplineQuadrature(bp)
+Sp = stiffness_matrix(qp)
+N = nbasis(bp)
+maximum(abs, Sp * ones(N))
+```
+
+By default [`mass_operator`](@ref) refuses it, since a singular assembly is usually too coarse
+a quadrature rather than an intended one:
+
+```@example asm
+try
+    mass_operator(Sp, bp)
+catch err
+    println(err.msg)
+end
+```
+
+`kernel = :project` asserts that the kernel is the constants, and asks for the solution that
+has no constant component. Each representation deflates that subspace in the way its own
+structure allows — [`CirculantMass`](@ref) gives the constant mode a zero factor, which is the
+Moore–Penrose pseudoinverse, and [`FactorizedMass`](@ref) factorises the minor that drops the
+first degree of freedom and takes the mean out afterwards. Both verify the assertion, and
+neither changes the matrix:
+
+```@example asm
+op = mass_operator(Sp, bp; kernel = :project)
+rhs = sin.(2π .* (0:(N - 1)) ./ N)     # one whole period, so its mean is already zero
+φ = op \ rhs
+maximum(abs, Sp * φ - rhs), sum(φ)     # it solves the system, in the mean-free gauge
+```
+
+The textbook cure is the rank-one shift by the mean projector,
+``\mathbb{S} + \mathbb{1}\mathbb{1}^T/N``, which is invertible and gives the same answer. It
+is also structurally full, so it costs ``O(N^2)`` to store and to factorise where the
+assembly itself costs ``O(N)``:
+
+```@example asm
+shifted = mass_operator(Matrix(Sp) .+ inv(N), bp)
+maximum(abs, φ - shifted \ rhs),
+Base.summarysize(Sp), Base.summarysize(Matrix(Sp) .+ inv(N))
+```
+
 ### Solving
 
 ```julia

scripts/mass_tolerance_margins.jl

@@ -0,0 +1,138 @@
+# What margin each of the three tolerances in `src/mass.jl` actually has, and therefore
+# whether the form they take is the right one.
+#
+# All three answer a structural question about an assembled matrix -- is this matrix
+# circulant, are the constants in its kernel, is the minor that drops one degree of freedom
+# invertible -- and the answer has to survive two things a written-down figure does not
+# survive: a rescaling of the assembly, and a change of element type. An absolute bound gets
+# both wrong. It calls an invertible matrix singular once the matrix is scaled below the
+# bound, and it calls every `Float32` assembly non-circulant because `Float32` rounding alone
+# exceeds a figure chosen for `Float64`.
+#
+# So each bound is relative and carries `eps(T)`. This script measures the two sides of each
+# one -- how close the accepted case comes to the bound, and how far the refused case sits
+# above it -- so that the form is chosen from evidence. The figures quoted in the comments in
+# `src/mass.jl` come from here.
+#
+# Run: julia --project=. scripts/mass_tolerance_margins.jl
+
+using SimpleSplines
+using LinearAlgebra
+using Printf
+using SparseArrays
+
+const MESHES = (UniformMesh, GradedMesh, RandomMesh)
+
+# ---------------------------------------------------------------------------------------
+# `_check_constant_kernel`: ‖M𝟙‖∞ against `n * eps(T) * ‖M‖∞`
+#
+# A stiffness matrix on a periodic basis has M𝟙 = 0 exactly in exact arithmetic -- the
+# integrand ∑ᵢ Bᵢ'(x) Bⱼ'(x) vanishes pointwise, since ∑ᵢ Bᵢ ≡ 1 -- so its residual is pure
+# rounding, and the question is only whether that rounding stays below the bound as n grows.
+function constant_kernel_margins(::Type{T}) where {T}
+    worst_singular, at_singular = zero(T), ""
+    best_invertible, at_invertible = T(Inf), ""
+
+    for meshtype in MESHES, p in 1:8, n in (2p + 2, 16, 33, 64, 256, 512, 1024)
+        n ≤ 2p + 1 && continue
+        q = SplineQuadrature(PeriodicBSplineBasis(meshtype(n, 2π), p))
+        for (A, singular) in ((T.(stiffness_matrix(q)), true), (T.(mass_matrix(q)), false))
+            ratio = norm(A * ones(T, n), Inf) / (n * eps(T) * norm(A, Inf))
+            if singular
+                ratio > worst_singular && ((worst_singular, at_singular) = (
+                    ratio, "$(meshtype) p=$p n=$n"))
+            else
+                ratio < best_invertible && ((best_invertible, at_invertible) = (
+                    ratio, "$(meshtype) p=$p n=$n"))
+            end
+        end
+    end
+    return (worst_singular, at_singular, best_invertible, at_invertible)
+end
+
+# ---------------------------------------------------------------------------------------
+# `_check_circulant`: max|M[i,j] - c[mod1(i-j+1,n)]| against `rtol * max|c|`
+#
+# The accepted side is a uniform periodic mesh, whose residual is assembly rounding. The
+# refused side is a graded mesh, which misses by a fraction of the entries themselves.
+function circulance_margins(::Type{T}) where {T}
+    worst_circulant, at_circulant = zero(T), ""
+    best_refused, at_refused = T(Inf), ""
+
+    for p in 1:8, n in (max(12, 2p + 2), 16, 33, 64, 256)
+
+        n ≤ 2p + 1 && continue
+        for (meshtype, circulant) in ((UniformMesh, true), (GradedMesh, false))
+            M = T.(Matrix(mass_matrix(SplineQuadrature(
+                PeriodicBSplineBasis(meshtype(n, 2π), p)))))
+            c = M[:, 1]
+            r = maximum(abs, [M[i, j] - c[mod1(i - j + 1, n)] for i in 1:n, j in 1:n]) /
+                maximum(abs, c)
+            if circulant
+                r > worst_circulant && ((worst_circulant, at_circulant) = (r, "p=$p n=$n"))
+            else
+                r < best_refused && ((best_refused, at_refused) = (r, "p=$p n=$n"))
+            end
+        end
+    end
+    return (worst_circulant, at_circulant, best_refused, at_refused)
+end
+
+# ---------------------------------------------------------------------------------------
+# `_check_definite_minor`: the pivot ratio of the Cholesky factor of the minor
+#
+# CHOLMOD reports success for a positive *semi*definite factorisation, so `issuccess` does not
+# separate a minor that is invertible from one that is not. The pivots do. `diag` of an LLᵀ
+# factor is the diagonal of L, so the pivot ratio is the square of the diagonal ratio; the
+# check compares the diagonal ratio against `sqrt(n * eps(T))`.
+#
+# The refused case is the projector onto the complement of span{𝟙, v} with v = (1,-1,1,-1):
+# it has the constants in its kernel, so it passes `_check_constant_kernel`, and v there as
+# well, so deflating the constants alone leaves a singular problem behind.
+function minor_pivot_margins()
+    worst_genuine, at_genuine = Inf, ""
+
+    for meshtype in MESHES, p in 1:4, n in (16, 33, 64, 256)
+        S = stiffness_matrix(SplineQuadrature(PeriodicBSplineBasis(meshtype(n, 2π), p)))
+        d = diag(cholesky(Symmetric(S[2:end, 2:end]); check = false))
+        lo, hi = extrema(abs, d)
+        ratio = (lo / hi)^2
+        ratio < worst_genuine && ((worst_genuine, at_genuine) = (
+            ratio, "$(meshtype) p=$p n=$n"))
+    end
+
+    v = [1.0, -1.0, 1.0, -1.0]
+    P = Matrix(I, 4, 4) .- ones(4, 4) ./ 4 .- (v * v') ./ 4
+    F = cholesky(Symmetric(sparse(P[2:end, 2:end])); check = false)
+    lo, hi = extrema(abs, diag(F))
+
+    return (worst_genuine, at_genuine, (lo / hi)^2, issuccess(F))
+end
+
+function main()
+    for T in (Float64, Float32)
+        ws, as, bi, ai = constant_kernel_margins(T)
+        println("\n", T, " -- the constants in the kernel, ‖M𝟙‖∞ / (n·eps·‖M‖∞)")
+        @printf("  accepted (stiffness) worst  %10.3e   at %s\n", ws, as)
+        @printf("  refused  (mass)      best   %10.3e   at %s\n", bi, ai)
+        @printf("  the gate is at 10; margins are %.1f× below and %.1e× above\n",
+            10 / ws, bi / 10)
+
+        wc, ac, br, ar = circulance_margins(T)
+        println("\n", T, " -- circulance, max deviation relative to max|c|")
+        @printf("  accepted (uniform)   worst  %10.3e = %7.1f·eps   at %s\n",
+            wc, wc / eps(T), ac)
+        @printf("  refused  (graded)    best   %10.3e                at %s\n", br, ar)
+        @printf("  the default rtol is sqrt(eps) = %.3e\n", sqrt(eps(T)))
+    end
+
+    wg, ag, bad, success = minor_pivot_margins()
+    println("\nthe minor's pivot ratio, min/max")
+    @printf("  genuine assembly     worst  %10.3e   at %s\n", wg, ag)
+    @printf("  kernel span{𝟙, v}           %10.3e   (issuccess reports %s)\n", bad, success)
+    @printf("  the bound is n·eps, i.e. %.3e at n = 3\n", 3 * eps())
+
+    return nothing
+end
+
+main()