From 622e974ab61b2f805ad7ff47172f7a0c93347c9a Mon Sep 17 00:00:00 2001 From: Michael Kraus Date: Mon, 14 Sep 2026 15:03:28 +0200 Subject: [PATCH 1/4] Add a kernel keyword for singular assemblies The periodic stiffness matrix has the constants in its kernel, so it could not be handed to a MassOperator at all. `kernel = :project` deflates that kernel instead of rejecting the matrix: CirculantMass gives the constant Fourier mode a zero factor, which is the Moore-Penrose pseudoinverse, and FactorizedMass factorises the minor that drops the first degree of freedom and fixes the gauge afterwards. Both verify the assertion rather than trusting it. It is what a caller needs in place of the rank-one shift S + 11'/N, which removes the singularity by making the matrix structurally full. The circulant invertibility check also used an absolute tolerance, so it passed the ~1e-14 zero eigenvalue of a stiffness matrix and then divided by it. It is now relative, at N * eps times the largest eigenvalue. Co-Authored-By: Claude Opus 5 (1M context) --- CHANGELOG.md | 45 ++++++++++ docs/src/usage/quadrature.md | 50 +++++++++++ src/mass.jl | 155 ++++++++++++++++++++++++++++++----- test/mass_tests.jl | 73 +++++++++++++++++ 4 files changed, 302 insertions(+), 21 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 8247a86..751eab3 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -5,6 +5,51 @@ 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 from 11 264 stored entries +(188 584 B) to 1 048 576 (16 785 576 B), and construction from `O(N)` to `O(N²)`. The +deflation leaves the matrix untouched. Measured: the deflated solve reproduces the shifted one +to 1e-13, and returns a solution whose mean is zero to 1e-17. + +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. + +### 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. + ## [0.1.0] — 2026-09-07 The first release with any implementation in it. Before this the package was a diff --git a/docs/src/usage/quadrature.md b/docs/src/usage/quadrature.md index 1ffbf52..f57194f 100644 --- a/docs/src/usage/quadrature.md +++ b/docs/src/usage/quadrature.md @@ -259,6 +259,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 diff --git a/src/mass.jl b/src/mass.jl index ae81dd6..dc1c403 100644 --- a/src/mass.jl +++ b/src/mass.jl @@ -22,6 +22,18 @@ Three representations are provided, and which one is built is decided by the bas All three answer `\`, `ldiv!` and `Matrix`. +A mass matrix is positive definite, but the same three representations carry the other +assemblies of a [`SplineQuadrature`](@ref), and a stiffness matrix on a periodic basis is +*singular*: the constants the basis represents lie in its kernel. The `kernel` keyword says +what to do about that. `:reject` is the default and throws, since a singular assembly is +usually too coarse a quadrature rather than an intended one. `:project` states that the +kernel is the constants and asks for the solution that has no constant component — the +mean-free solution — which is what ``-\phi'' = \rho`` on a periodic domain asks for. Both +periodic representations implement it, each in the way its own structure allows, and both +verify the assertion rather than trusting it. Neither forms the rank-one shift +``\mathbb{M} + \mathbb{1}\mathbb{1}^T/N`` that would remove the singularity by filling the +matrix in completely. + !!! note "Where this does and does not matter" A mass solve is a small part of the cost of these discretisations — at ``N = 384`` it is around ``0.06`` ms against a ``2`` ms implicit step. The transform is used because it is @@ -34,7 +46,7 @@ abstract type MassOperator{T} end Base.eltype(::MassOperator{T}) where {T} = T @doc raw""" - FactorizedMass(M) + FactorizedMass(M; kernel = :reject) The sparse Cholesky factorisation of the mass matrix, for meshes on which it is not circulant. @@ -42,20 +54,55 @@ circulant. ``\mathbb{M}`` is symmetric positive definite and banded modulo ``N`` — a basis function overlaps only the ``2p+1`` others whose supports meet its own — so the factorisation is cheap and the solve is ``O(Np)``. + +With `kernel = :project` the matrix may be singular instead, with the constants in its +kernel. What is factorised is then the minor that drops the first degree of freedom, which is +positive definite exactly when the kernel is the constants and nothing more: a vector +supported away from the first index lies in the span of ``\mathbb{1}`` only if it is zero. +The dropped degree of freedom is the gauge, and [`mass_solve!`](@ref) fixes it afterwards by +taking the mean out. Deleting a row and a column preserves the sparsity; the rank-one shift +``\mathbb{M} + \mathbb{1}\mathbb{1}^T/N``, which also removes the singularity, makes the +matrix structurally full and the factorisation ``O(N^2)``. """ -struct FactorizedMass{T, MT, FT} <: MassOperator{T} +struct FactorizedMass{T, MT, FT, K} <: MassOperator{T} M::MT fact::FT - function FactorizedMass(M::MT) where {T, MT <: AbstractMatrix{T}} - F = cholesky(Symmetric(M); check = false) - issuccess(F) || throw(ArgumentError( - "the mass matrix is not positive definite; the quadrature is too coarse to " * - "resolve the basis")) - new{T, MT, typeof(F)}(M, F) + function FactorizedMass(M::MT; kernel = :reject) where {T, MT <: AbstractMatrix{T}} + if kernel === :reject + F = cholesky(Symmetric(M); check = false) + issuccess(F) || throw(ArgumentError( + "the mass matrix is not positive definite; the quadrature is too coarse to " * + "resolve the basis")) + return new{T, MT, typeof(F), :reject}(M, F) + elseif kernel === :project + _check_constant_kernel(M) + F = cholesky(Symmetric(M[2:end, 2:end]); check = false) + issuccess(F) || throw(ArgumentError( + "the matrix is singular beyond the constants: the minor that drops the " * + "first degree of freedom is not positive definite either")) + return new{T, MT, typeof(F), :project}(M, F) + end + throw(_not_a_kernel_mode(kernel)) end end +function _not_a_kernel_mode(kernel) + ArgumentError("kernel must be :reject or :project, not $(repr(kernel))") +end + +# `:project` deflates the constants specifically, so it is correct only where that is what the +# kernel holds. Verify rather than assume: a matrix singular for some other reason would take +# the same path, and both the positive-definite minor and the mean-free gauge would then be +# answering a question nobody asked. +function _check_constant_kernel(M::AbstractMatrix{T}) where {T} + residual = norm(M * ones(T, size(M, 2)), Inf) + residual ≤ sqrt(eps(T)) * max(one(T), norm(M, Inf)) || throw(ArgumentError( + "kernel = :project deflates the constants, but M𝟙 has norm $(residual); the kernel " * + "of this matrix is not what the projection assumes")) + return nothing +end + @doc raw""" BandedMass(M) @@ -108,7 +155,7 @@ function _bandwidth(M::AbstractMatrix) end @doc raw""" - CirculantMass(M, n) + CirculantMass(M, n; kernel = :reject) The mass matrix of a uniform periodic mesh, represented by the eigenvalues of its circulant structure and a pair of planned real transforms. @@ -131,16 +178,24 @@ does not. The first column is read from the assembled matrix rather than recomputed, and the construction checks that the matrix really is circulant — a silent mismatch here would give wrong answers on every mesh that is uniform by accident rather than by construction. + +`kernel = :project` accepts a matrix that is singular with the constants in its kernel, and +gives the constant mode a zero factor instead of an infinite one. That is the Moore–Penrose +pseudoinverse: the solve drops the constant component of the right-hand side and returns the +solution that has none of it. The transform does the whole of the work, so nothing is added +to the matrix and nothing is taken out of it. This is what makes a periodic stiffness matrix +usable here, and it is the rule an FFT Poisson solver applies when it sets the ``k = 0`` +factor to zero rather than dividing by it. """ struct CirculantMass{T, MT, PT, IT} <: MassOperator{T} M::MT - ĉ::Vector{Complex{T}} + ĉ⁻¹::Vector{Complex{T}} plan::PT iplan::IT buf::Vector{Complex{T}} n::Int - function CirculantMass(M::MT, n::Integer; atol = 1e-10) where { + function CirculantMass(M::MT, n::Integer; atol = 1e-10, kernel = :reject) where { T, MT <: AbstractMatrix{T}} size(M, 1) == n || throw(DimensionMismatch( "the mass matrix is $(size(M, 1))×$(size(M, 2)) but n = $(n)")) @@ -164,13 +219,43 @@ struct CirculantMass{T, MT, PT, IT} <: MassOperator{T} buf = Vector{Complex{T}}(undef, n ÷ 2 + 1) plan = plan_rfft(c; flags = FFTW.ESTIMATE | FFTW.UNALIGNED) iplan = plan_irfft(buf, n; flags = FFTW.ESTIMATE | FFTW.UNALIGNED) - ĉ = plan * c + ĉ⁻¹ = _reciprocal_eigenvalues(plan * c, Int(n), kernel) - all(x -> abs(x) > eps(T), ĉ) || throw(ArgumentError( - "the circulant mass matrix has a zero eigenvalue and is not invertible")) + new{T, MT, typeof(plan), typeof(iplan)}(M, ĉ⁻¹, plan, iplan, buf, Int(n)) + end +end - new{T, MT, typeof(plan), typeof(iplan)}(M, ĉ, plan, iplan, buf, Int(n)) +# The solve multiplies by these rather than dividing by the eigenvalues themselves. It moves +# n ÷ 2 + 1 divisions out of every solve and into the construction, and it is what lets the +# kernel carry a zero factor rather than a branch in the inner loop. +# +# One tolerance decides both modes, and it is *relative*: `n * eps` below the largest +# eigenvalue is where `pinv` puts the same cutoff, and an eigenvalue there is not +# distinguishable from zero by a transform that accumulated its own rounding on the way. An +# absolute bound instead calls a singular matrix invertible whenever its entries are large +# enough -- a periodic stiffness matrix has `|ĉ₀| ≈ 10⁻¹⁴` rather than zero, and the solve +# that follows divides by it. +function _reciprocal_eigenvalues(ĉ::Vector{Complex{T}}, n::Int, kernel) where {T} + tol = n * eps(T) * maximum(abs, ĉ) + if kernel === :reject + all(x -> abs(x) > tol, ĉ) || throw(ArgumentError( + "the circulant mass matrix has a zero eigenvalue and is not invertible; pass " * + "kernel = :project to solve in the complement of the kernel instead")) + return inv.(ĉ) + elseif kernel === :project + # The constant mode is ĉ[1] and nothing else, so this deflates the same subspace the + # factorised representation deflates, and it verifies that the subspace is the whole + # of the kernel rather than zeroing whichever mode happens to come out small. + abs(ĉ[1]) ≤ tol || throw(ArgumentError( + "kernel = :project deflates the constants, but the constant mode of this " * + "matrix has eigenvalue $(abs(ĉ[1])); its kernel is not what the projection " * + "assumes")) + all(x -> abs(x) > tol, @view ĉ[2:end]) || throw(ArgumentError( + "the matrix is singular beyond the constants: a second mode of the transform " * + "vanishes as well")) + return [i == 1 ? zero(z) : inv(z) for (i, z) in enumerate(ĉ)] end + throw(_not_a_kernel_mode(kernel)) end function _not_circulant(atol) @@ -239,11 +324,28 @@ function mass_solve! end # is the only representation here that does, and it is reached only by a periodic basis on a # non-uniform mesh -- the one case that is neither circulant nor banded. A bounded basis on # the same mesh goes through `BandedMass`, which allocates nothing. -function mass_solve!(y::AbstractVector, op::FactorizedMass, x::AbstractVector) +function mass_solve!(y::AbstractVector, op::FactorizedMass{T, MT, FT, :reject}, + x::AbstractVector) where {T, MT, FT} copyto!(y, op.fact \ x) return y end +# The constants are in the kernel, so the right-hand side determines a solution only once its +# own constant part is out of the way, and determines it only up to a constant. The first +# degree of freedom is the one the factorisation dropped: it is set to zero, which picks one +# solution, and the mean is taken out afterwards, which picks the one in the complement of +# the kernel. Both steps are O(N), and the temporary is the one the CHOLMOD solve allocates +# in any case. +function mass_solve!(y::AbstractVector, op::FactorizedMass{T, MT, FT, :project}, + x::AbstractVector) where {T, MT, FT} + n = length(x) + ŷ = op.fact \ @views(x[2:n] .- sum(x) / n) + y[1] = zero(eltype(y)) + copyto!(view(y, 2:n), ŷ) + y .-= sum(y) / n + return y +end + # The banded factor does implement an in-place `ldiv!`, so this path allocates nothing. function mass_solve!(y::AbstractVector, op::BandedMass, x::AbstractVector) y === x || copyto!(y, x) @@ -253,7 +355,7 @@ end function mass_solve!(y::AbstractVector, op::CirculantMass, x::AbstractVector) mul!(op.buf, op.plan, x) - op.buf ./= op.ĉ + op.buf .*= op.ĉ⁻¹ mul!(y, op.iplan, op.buf) return y end @@ -276,13 +378,18 @@ function Base.:\(op::MassOperator, X::AbstractMatrix) end @doc raw""" - mass_operator(M, basis) + mass_operator(M, basis; kernel = :reject) Build the [`MassOperator`](@ref) appropriate to `basis`: a [`CirculantMass`](@ref) for a periodic basis on a [`UniformMesh`](@ref), a [`FactorizedMass`](@ref) for a periodic basis on any other mesh, and a [`BandedMass`](@ref) for a bounded basis — clamped or recombined — whose mass matrix has no seam to wrap across. +`kernel = :project` passes the deflation on to whichever of the two periodic representations +is chosen, so that a caller with a singular assembly says what it wants rather than which +representation implements it. See [`MassOperator`](@ref) for what the deflation means and +each representation for how it is done. + !!! note "Circulance is a property of the basis, not of the mesh" A uniform mesh is necessary but not sufficient. The mass matrix is circulant only when every basis function is a translate of one cardinal spline, which needs the *periodic* @@ -292,8 +399,14 @@ whose mass matrix has no seam to wrap across. wrong answers everywhere except in [`CirculantMass`](@ref)'s own verification, which would reject it. """ -mass_operator(M::AbstractMatrix, ::AbstractBSplineBasis) = BandedMass(M) +function mass_operator(M::AbstractMatrix, ::AbstractBSplineBasis; kernel = :reject) + kernel === :reject || throw(ArgumentError( + "kernel = $(repr(kernel)) is implemented for a periodic basis only; a bounded basis " * + "takes the banded representation, which carries no deflation")) + return BandedMass(M) +end -function mass_operator(M::AbstractMatrix, b::PeriodicBSplineBasis) - mesh(b) isa UniformMesh ? CirculantMass(M, nbasis(b)) : FactorizedMass(M) +function mass_operator(M::AbstractMatrix, b::PeriodicBSplineBasis; kernel = :reject) + mesh(b) isa UniformMesh ? CirculantMass(M, nbasis(b); kernel) : + FactorizedMass(M; kernel) end diff --git a/test/mass_tests.jl b/test/mass_tests.jl index fca917e..f211493 100644 --- a/test/mass_tests.jl +++ b/test/mass_tests.jl @@ -57,6 +57,79 @@ using Test end end + @testset "$(rpad("a singular assembly is deflated rather than shifted",76))" begin + # the stiffness matrix of a periodic basis has the constants in its kernel, on every + # mesh. The textbook cure is the rank-one shift by the mean projector, which is what + # the deflation is checked against -- it must agree with it without ever forming it. + for meshtype in (UniformMesh, GradedMesh), p in 1:4, n in (16, 33) + b = PeriodicBSplineBasis(meshtype(n, 2π), p) + S = stiffness_matrix(SplineQuadrature(b)) + @test norm(S * ones(n), Inf) < 1e-10 + + op = mass_operator(S, b; kernel = :project) + @test op isa (meshtype === UniformMesh ? CirculantMass : FactorizedMass) + @test size(op) == (n, n) + + shifted = Matrix(S) .+ inv(n) + for _ in 1:3 + rhs = randn(n) + meanfree = rhs .- sum(rhs) / n + x = op \ rhs + + @test x ≈ shifted \ meanfree atol = 1e-10 + @test S * x ≈ meanfree atol = 1e-10 # it solves the system + @test sum(x) / n ≈ 0 atol = 1e-12 # in the complement of the kernel + + # the mean of the right-hand side is dropped by the deflation itself, so + # taking it out first changes nothing + @test op \ meanfree ≈ x atol = 1e-12 + + y = copy(rhs) # y and x may alias + @test mass_solve!(y, op, y) ≈ x atol = 1e-12 + end + end + end + + @testset "$(rpad("the deflation asserts a kernel, and both paths check it",76))" begin + # :project says what the kernel is rather than asking what it happens to be, so a + # mass matrix -- invertible, and with M𝟙 nowhere near zero -- is refused by it on + # either representation. That is what keeps the two meanings the same one. + for m in (UniformMesh(32, 2π), GradedMesh(32, 2π)), p in 1:4 + + b = PeriodicBSplineBasis(m, p) + M = mass_matrix(SplineQuadrature(b)) + @test_throws ArgumentError mass_operator(M, b; kernel = :project) + end + end + + @testset "$(rpad("what the deflation refuses",76))" begin + n = 16 + b = PeriodicBSplineBasis(UniformMesh(n, 2π), 3) + S = stiffness_matrix(SplineQuadrature(b)) + + # the default rejects a singular matrix, and the tolerance is relative: the zero + # eigenvalue of a stiffness matrix comes out of the transform at ~1e-14, not at zero, + # so an absolute bound of eps() would call it invertible and divide by it + @test_throws ArgumentError mass_operator(S, b) + @test_throws ArgumentError mass_operator(S, b; kernel = :reject) + + # a bounded basis has no deflation to reach + @test_throws ArgumentError mass_operator( + Matrix(mass_matrix(SplineQuadrature(BSplineBasis(UniformMesh(n, 2π), 3)))), + BSplineBasis(UniformMesh(n, 2π), 3); kernel = :project) + + # and neither mode accepts a name that is neither + @test_throws ArgumentError mass_operator(S, b; kernel = :ignore) + @test_throws ArgumentError mass_operator( + stiffness_matrix(SplineQuadrature(PeriodicBSplineBasis(GradedMesh(n, 2π), 3))), + PeriodicBSplineBasis(GradedMesh(n, 2π), 3); kernel = :ignore) + + # the deflation removes the constants specifically, so it verifies that they are + # what the kernel holds rather than assuming it + @test_throws ArgumentError SimpleSplines.FactorizedMass( + sparse(Diagonal([1.0, 0.0, 1.0, 1.0])); kernel = :project) + end + @testset "$(rpad("in-place solve, aliasing, and no allocation",76))" begin q = SplineQuadrature(PeriodicBSplineBasis(UniformMesh(32, 2π), 3)) op = mass_operator(q) From f49f39aa503f5a7c28ba0e6dea1d3de786ac1e05 Mon Sep 17 00:00:00 2001 From: Michael Kraus Date: Mon, 14 Sep 2026 15:24:02 +0200 Subject: [PATCH 2/4] Restore concrete inference, correct the prose `FactorizedMass`'s two-branch constructor returned a two-way Union even for a literal keyword, and it propagated into SplineQuadrature's own type parameter, which was concrete before this branch. `Base.@constprop :aggressive` restores it; a runtime Symbol still gives the Union, which is correct. The solve no longer divides by the eigenvalues, so the three places that described it as a pointwise division say multiplication instead. The changelog's stored-entry pair was measured at degree 5 rather than the degree 4 it names, and its two agreement figures were the best case rather than the range. Co-Authored-By: Claude Opus 5 (1M context) --- CHANGELOG.md | 9 +++++---- docs/src/usage/quadrature.md | 2 +- src/mass.jl | 12 +++++++----- 3 files changed, 13 insertions(+), 10 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 751eab3..f66e12b 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -23,10 +23,11 @@ Fourier mode is the one that vanishes and that no other does, and `FactorizedMas 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 from 11 264 stored entries -(188 584 B) to 1 048 576 (16 785 576 B), and construction from `O(N)` to `O(N²)`. The -deflation leaves the matrix untouched. Measured: the deflated solve reproduces the shifted one -to 1e-13, and returns a solution whose mean is zero to 1e-17. +on 1024 cells the shift takes a periodic stiffness matrix from 9 216 stored entries (155 816 B) +to 1 048 576 (16 785 576 B — it stays a `SparseMatrixCSC`, so every entry now carries a row +index beside it), 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. diff --git a/docs/src/usage/quadrature.md b/docs/src/usage/quadrature.md index f57194f..ff14daa 100644 --- a/docs/src/usage/quadrature.md +++ b/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 | diff --git a/src/mass.jl b/src/mass.jl index dc1c403..f079580 100644 --- a/src/mass.jl +++ b/src/mass.jl @@ -10,7 +10,7 @@ Three representations are provided, and which one is built is decided by the bas - [`CirculantMass`](@ref) for a periodic basis on a [`UniformMesh`](@ref), where the basis functions are translates of a single cardinal spline and ``\mathbb{M}`` is *circulant*, hence diagonalised by the discrete Fourier transform. A solve is two transforms and a - pointwise division, ``O(N \log N)``, with the transforms planned once. + pointwise multiplication, ``O(N \log N)``, with the transforms planned once. - [`BandedMass`](@ref) for a bounded basis. The overlaps are contiguous — there is no seam — so the matrix is banded outright and a banded Cholesky solves it in ``O(Np)`` with no allocation at all. @@ -68,7 +68,8 @@ struct FactorizedMass{T, MT, FT, K} <: MassOperator{T} M::MT fact::FT - function FactorizedMass(M::MT; kernel = :reject) where {T, MT <: AbstractMatrix{T}} + Base.@constprop :aggressive function FactorizedMass(M::MT; + kernel = :reject) where {T, MT <: AbstractMatrix{T}} if kernel === :reject F = cholesky(Symmetric(M); check = false) issuccess(F) || throw(ArgumentError( @@ -169,7 +170,7 @@ On a uniform mesh every basis function is a translate of one cardinal spline, so ``` with ``F`` the discrete Fourier transform. A solve is therefore a forward transform, a -pointwise division and an inverse transform. +pointwise multiplication by the reciprocals of ``\hat{c}``, and an inverse transform. The plans are created once, at construction, and the spectral buffer is preallocated, so a solve allocates nothing beyond its result. `\` allocates the result; [`mass_solve!`](@ref) @@ -334,8 +335,9 @@ end # own constant part is out of the way, and determines it only up to a constant. The first # degree of freedom is the one the factorisation dropped: it is set to zero, which picks one # solution, and the mean is taken out afterwards, which picks the one in the complement of -# the kernel. Both steps are O(N), and the temporary is the one the CHOLMOD solve allocates -# in any case. +# the kernel. Both steps are O(N). Two temporaries are allocated: the mean-free right-hand +# side, which is what the factor is given, and the one the CHOLMOD solve produces in any +# case. function mass_solve!(y::AbstractVector, op::FactorizedMass{T, MT, FT, :project}, x::AbstractVector) where {T, MT, FT} n = length(x) From 9329b6abf828821a5ac1a096799ee2c3ea296d2f Mon Sep 17 00:00:00 2001 From: Michael Kraus Date: Mon, 14 Sep 2026 16:26:42 +0200 Subject: [PATCH 3/4] Make the mass tolerances relative and type aware MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Three tolerances in the mass operators were absolute, and each accepted something it should have refused. The constant-kernel check of `kernel = :project` could not fall below `sqrt(eps(T))`, so it was absolute for any matrix smaller than that: an invertible mass matrix scaled below 1e-8 was accepted as having the constants in its kernel, and the solve returned a vector with relative residual of order one. The bound is now `N * eps(T)` times `‖M‖∞` with no floor, which is the standard the circulant path already used. The circulance check used `atol = 1e-10`, below `Float32` rounding, so no `Float32` assembly was circulant at all and `CirculantMass` could not be built in that element type. The keyword is now `rtol`, relative to the largest entry of the first column, defaulting to `sqrt(eps(T))`. `FactorizedMass` took CHOLMOD's `issuccess` as proof that the minor is positive definite. CHOLMOD calls a positive semidefinite factorisation a success, so a matrix whose kernel is one dimension larger than the constants was accepted and the solve returned amplified noise. The pivots are now read: the smallest must exceed `N * eps(T)` of the largest. An unknown `kernel` name on a bounded basis now reports that it is not a mode, rather than that it is implemented for a periodic basis. Co-Authored-By: Claude Opus 5 (1M context) --- CHANGELOG.md | 37 ++++++++-- docs/src/usage/quadrature.md | 8 ++- src/mass.jl | 98 +++++++++++++++++++------- test/mass_tests.jl | 130 +++++++++++++++++++++++++++++++++++ 4 files changed, 241 insertions(+), 32 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index f66e12b..246b13e 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -23,11 +23,12 @@ Fourier mode is the one that vanishes and that no other does, and `FactorizedMas 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 from 9 216 stored entries (155 816 B) -to 1 048 576 (16 785 576 B — it stays a `SparseMatrixCSC`, so every entry now carries a row -index beside it), 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. +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. @@ -41,6 +42,27 @@ what came back was amplified noise rather than an error. The tolerance is now re `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 `N * eps(T)` times `‖M‖∞`, with no floor, which is the + standard the circulant path already used. + - 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 @@ -51,6 +73,11 @@ A `CirculantMass` solve multiplies by precomputed reciprocal eigenvalues instead 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 diff --git a/docs/src/usage/quadrature.md b/docs/src/usage/quadrature.md index ff14daa..e2ec3cf 100644 --- a/docs/src/usage/quadrature.md +++ b/docs/src/usage/quadrature.md @@ -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())) diff --git a/src/mass.jl b/src/mass.jl index f079580..518d81e 100644 --- a/src/mass.jl +++ b/src/mass.jl @@ -79,9 +79,7 @@ struct FactorizedMass{T, MT, FT, K} <: MassOperator{T} elseif kernel === :project _check_constant_kernel(M) F = cholesky(Symmetric(M[2:end, 2:end]); check = false) - issuccess(F) || throw(ArgumentError( - "the matrix is singular beyond the constants: the minor that drops the " * - "first degree of freedom is not positive definite either")) + _check_definite_minor(F) return new{T, MT, typeof(F), :project}(M, F) end throw(_not_a_kernel_mode(kernel)) @@ -92,15 +90,49 @@ function _not_a_kernel_mode(kernel) ArgumentError("kernel must be :reject or :project, not $(repr(kernel))") end +function _singular_beyond_the_constants() + ArgumentError( + "the matrix is singular beyond the constants: the minor that drops the first degree " * + "of freedom is not positive definite either") +end + +# CHOLMOD calls the factorisation of a positive *semi*definite matrix a success -- a zero pivot +# is not a breakdown to it -- so `issuccess` alone does not say whether the minor is +# invertible, and a kernel one dimension larger than the constants passes it. The whole +# argument for `:project` rests on that minor, so the pivots are read rather than assumed. +# +# `diag` of an LLᵀ factor is the diagonal of L, so the pivots are its squares, and the minor is +# positive definite to working precision exactly when the smallest is above `n * eps` of the +# largest -- the same relative standard `_check_constant_kernel` and `_reciprocal_eigenvalues` +# hold their own assertions to. A matrix whose kernel is `span{𝟙, v}` puts that ratio at 2e-16, +# where a genuine assembly holds it above 0.19 over degrees 1 to 4 and 16 to 256 cells. +function _check_definite_minor(F) + issuccess(F) || throw(_singular_beyond_the_constants()) + d = diag(F) + lo, hi = extrema(abs, d) + lo > sqrt(length(d) * eps(eltype(d))) * hi || throw(_singular_beyond_the_constants()) + return nothing +end + # `:project` deflates the constants specifically, so it is correct only where that is what the # kernel holds. Verify rather than assume: a matrix singular for some other reason would take # the same path, and both the positive-definite minor and the mean-free gauge would then be # answering a question nobody asked. +# +# The bound is `n * eps(T)` relative to ‖M‖∞, which is the row-sum scale `M𝟙` is measured on, +# and it is the standard `_reciprocal_eigenvalues` holds the circulant path to. It has to be +# *relative* and it has to have no floor: a bound that cannot fall below one is absolute for +# every matrix smaller than that, and accepts any invertible matrix scaled down far enough. +# `M𝟙` of a genuine singular assembly is pure rounding and grows with n like the bound does -- +# over degrees 1 to 8 and 64 to 4096 cells it reaches 0.16 of it, where an invertible mass +# matrix sits 10¹³ above it in the same units. function _check_constant_kernel(M::AbstractMatrix{T}) where {T} - residual = norm(M * ones(T, size(M, 2)), Inf) - residual ≤ sqrt(eps(T)) * max(one(T), norm(M, Inf)) || throw(ArgumentError( - "kernel = :project deflates the constants, but M𝟙 has norm $(residual); the kernel " * - "of this matrix is not what the projection assumes")) + n = size(M, 2) + residual = norm(M * ones(T, n), Inf) + tol = n * eps(T) * norm(M, Inf) + residual ≤ tol || throw(ArgumentError( + "kernel = :project deflates the constants, but M𝟙 has norm $(residual) against a " * + "tolerance of $(tol); the kernel of this matrix is not what the projection assumes")) return nothing end @@ -156,7 +188,7 @@ function _bandwidth(M::AbstractMatrix) end @doc raw""" - CirculantMass(M, n; kernel = :reject) + CirculantMass(M, n; kernel = :reject, rtol = sqrt(eps(T))) The mass matrix of a uniform periodic mesh, represented by the eigenvalues of its circulant structure and a pair of planned real transforms. @@ -178,7 +210,9 @@ does not. The first column is read from the assembled matrix rather than recomputed, and the construction checks that the matrix really is circulant — a silent mismatch here would give -wrong answers on every mesh that is uniform by accident rather than by construction. +wrong answers on every mesh that is uniform by accident rather than by construction. `rtol` +is the tolerance of that check, *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. `kernel = :project` accepts a matrix that is singular with the constants in its kernel, and gives the constant mode a zero factor instead of an infinite one. That is the Moore–Penrose @@ -187,6 +221,11 @@ solution that has none of it. The transform does the whole of the work, so nothi to the matrix and nothing is taken out of it. This is what makes a periodic stiffness matrix usable here, and it is the rule an FFT Poisson solver applies when it sets the ``k = 0`` factor to zero rather than dividing by it. + +The deflation lives entirely in the stored reciprocal eigenvalues, so the type carries no +kernel-mode parameter and there is one [`mass_solve!`](@ref) for both modes. +[`FactorizedMass`](@ref) carries such a parameter because there the two modes are two +different solves, and the parameter is what selects between them. """ struct CirculantMass{T, MT, PT, IT} <: MassOperator{T} M::MT @@ -196,7 +235,8 @@ struct CirculantMass{T, MT, PT, IT} <: MassOperator{T} buf::Vector{Complex{T}} n::Int - function CirculantMass(M::MT, n::Integer; atol = 1e-10, kernel = :reject) where { + function CirculantMass(M::MT, n::Integer; rtol = sqrt(eps(T)), + kernel = :reject) where { T, MT <: AbstractMatrix{T}} size(M, 1) == n || throw(DimensionMismatch( "the mass matrix is $(size(M, 1))×$(size(M, 2)) but n = $(n)")) @@ -205,7 +245,7 @@ struct CirculantMass{T, MT, PT, IT} <: MassOperator{T} # Verify rather than assume. A matrix that is banded but not circulant would still # produce plausible numbers through the transform, and the failure would show up # much later as a wrong conservation law. - _check_circulant(M, c, Int(n), atol) + _check_circulant(M, c, Int(n), rtol) # The plans are made UNALIGNED so that they accept any strided argument -- a view # into a column of a matrix, in particular, whose alignment an aligned plan would @@ -259,19 +299,26 @@ function _reciprocal_eigenvalues(ĉ::Vector{Complex{T}}, n::Int, kernel) where { throw(_not_a_kernel_mode(kernel)) end -function _not_circulant(atol) +function _not_circulant(rtol) ArgumentError( - "the mass matrix is not circulant to within $(atol); a CirculantMass is only valid on " * - "a uniform mesh") + "the mass matrix is not circulant to a relative tolerance of $(rtol); a " * + "CirculantMass is only valid on a uniform mesh") end +# `rtol` is relative to the largest entry of the first column, so the verdict does not change +# when the assembly is scaled. It defaults to `sqrt(eps(T))`, which is what makes the question +# answerable in every element type: the residual of a genuinely circulant assembly is +# rounding, at most ~200·eps of that scale, while a graded mesh misses by 6e-2 of it. An +# absolute default instead fixes the answer to one element type -- at 1e-10 no `Float32` +# assembly is circulant at all, since `Float32` rounding alone exceeds it. +# # The generic check: every entry against the first column. Quadratic, which for a dense # matrix is what the question costs. -function _check_circulant(M::AbstractMatrix{T}, c::Vector{T}, n::Int, atol) where {T} - tol = atol * max(one(T), maximum(abs, c)) +function _check_circulant(M::AbstractMatrix{T}, c::Vector{T}, n::Int, rtol) where {T} + tol = rtol * maximum(abs, c) for j in 2:n, i in 1:n - abs(M[i, j] - c[mod1(i - j + 1, n)]) ≤ tol || throw(_not_circulant(atol)) + abs(M[i, j] - c[mod1(i - j + 1, n)]) ≤ tol || throw(_not_circulant(rtol)) end return nothing end @@ -282,8 +329,8 @@ end # entries in `c` covers the *unstored* positions as well: a column missing one of them cannot # reach the count. Both directions in O(nnz), which at n = 512 is 0.004 ms against the 9.1 ms # of probing all n^2 positions of a sparse matrix one `getindex` at a time. -function _check_circulant(M::SparseMatrixCSC{T}, c::Vector{T}, n::Int, atol) where {T} - tol = atol * max(one(T), maximum(abs, c)) +function _check_circulant(M::SparseMatrixCSC{T}, c::Vector{T}, n::Int, rtol) where {T} + tol = rtol * maximum(abs, c) rows = rowvals(M) vals = nonzeros(M) expected = count(x -> abs(x) > tol, c) @@ -292,10 +339,10 @@ function _check_circulant(M::SparseMatrixCSC{T}, c::Vector{T}, n::Int, atol) whe hits = 0 for t in nzrange(M, j) cᵣ = c[mod1(rows[t] - j + 1, n)] - abs(vals[t] - cᵣ) ≤ tol || throw(_not_circulant(atol)) + abs(vals[t] - cᵣ) ≤ tol || throw(_not_circulant(rtol)) abs(cᵣ) > tol && (hits += 1) end - hits == expected || throw(_not_circulant(atol)) + hits == expected || throw(_not_circulant(rtol)) end return nothing end @@ -402,9 +449,12 @@ each representation for how it is done. would reject it. """ function mass_operator(M::AbstractMatrix, ::AbstractBSplineBasis; kernel = :reject) - kernel === :reject || throw(ArgumentError( - "kernel = $(repr(kernel)) is implemented for a periodic basis only; a bounded basis " * - "takes the banded representation, which carries no deflation")) + # A name that is not a mode at all is reported as such. Reaching for the basis first would + # say that `kernel = :ignore` is implemented for a periodic basis, which it is not. + kernel === :project && throw(ArgumentError( + "kernel = :project is implemented for a periodic basis only; a bounded basis takes " * + "the banded representation, which carries no deflation")) + kernel === :reject || throw(_not_a_kernel_mode(kernel)) return BandedMass(M) end diff --git a/test/mass_tests.jl b/test/mass_tests.jl index f211493..5241a22 100644 --- a/test/mass_tests.jl +++ b/test/mass_tests.jl @@ -124,12 +124,142 @@ using Test stiffness_matrix(SplineQuadrature(PeriodicBSplineBasis(GradedMesh(n, 2π), 3))), PeriodicBSplineBasis(GradedMesh(n, 2π), 3); kernel = :ignore) + # a name that is not a mode is reported as one on a bounded basis too, rather than as + # a mode that some other basis implements + bb = BSplineBasis(UniformMesh(n, 2π), 3) + Mb = Matrix(mass_matrix(SplineQuadrature(bb))) + err = try + mass_operator(Mb, bb; kernel = :ignore) + catch e + e + end + @test err isa ArgumentError + @test occursin("must be :reject or :project", err.msg) + @test !occursin("periodic basis only", err.msg) + # the deflation removes the constants specifically, so it verifies that they are # what the kernel holds rather than assuming it @test_throws ArgumentError SimpleSplines.FactorizedMass( sparse(Diagonal([1.0, 0.0, 1.0, 1.0])); kernel = :project) end + @testset "$(rpad("a kernel larger than the constants is refused by both paths",76))" begin + # the projector onto the complement of span{𝟙, v} with v = (1,-1,1,-1): it has 𝟙 in + # its kernel, so it passes the constant-kernel check, and v there as well, so the + # deflation of the constants alone does not leave an invertible problem behind. Both + # v and 𝟙 are constant along diagonals, so the same matrix exercises both + # representations. + n = 4 + v = [1.0, -1.0, 1.0, -1.0] + P = Matrix(I, n, n) .- ones(n, n) ./ n .- (v * v') ./ n + + @test norm(P * ones(n), Inf) < 1e-14 # the constants are in the kernel + @test norm(P * v, Inf) < 1e-14 # and so is a second vector + + @test_throws ArgumentError CirculantMass(P, n; kernel = :project) + @test_throws ArgumentError SimpleSplines.FactorizedMass( + sparse(P); kernel = :project) + + # both say so, rather than reporting the kernel assertion they did accept + for e in (try + CirculantMass(P, n; kernel = :project) + catch e + e + end, + try + SimpleSplines.FactorizedMass(sparse(P); kernel = :project) + catch e + e + end) + @test occursin("singular beyond the constants", e.msg) + end + end + + @testset "$(rpad("the tolerances are relative, so scaling changes no verdict",76))" begin + # an absolute tolerance floored at one is absolute for every matrix below that scale, + # and it then accepts an invertible matrix as having the constants in its kernel + n = 32 + for meshtype in (UniformMesh, GradedMesh), p in 1:4, + s in (1e-12, 1e-8, 1.0, 1e8) + b = PeriodicBSplineBasis(meshtype(n, 2π), p) + q = SplineQuadrature(b) + + # an invertible assembly is refused at every scale ... + @test_throws ArgumentError mass_operator( + mass_matrix(q) .* s, b; kernel = :project) + + # ... and a singular one is accepted at every scale, and solves + S = stiffness_matrix(q) .* s + op = mass_operator(S, b; kernel = :project) + rhs = randn(n) + rhs .-= sum(rhs) / n + # `x` scales like 1/s, so the residual of `S * x` does not depend on s at all + x = op \ rhs + @test norm(S * x .- rhs, Inf) ≤ 1e-10 * norm(rhs, Inf) + @test abs(sum(x)) ≤ 1e-10 * norm(x, Inf) * n + end + + # circulance is judged relatively too: a uniform mesh stays circulant under scaling, + # and a graded one stays refused + M = Matrix(mass_matrix(SplineQuadrature( + PeriodicBSplineBasis(UniformMesh(16, 2π), 3)))) + G = Matrix(mass_matrix(SplineQuadrature( + PeriodicBSplineBasis(GradedMesh(16, 2π), 3)))) + for s in (1e-12, 1.0, 1e12) + @test CirculantMass(M .* s, 16) isa CirculantMass + @test_throws ArgumentError CirculantMass(G .* s, 16) + end + end + + @testset "$(rpad("the default tolerances follow the element type",76))" begin + # an absolute default fixes the answer to one element type: at 1e-10 no Float32 + # assembly is circulant at all, because Float32 rounding alone exceeds it + for p in 1:4, n in (16, 32, 64) + + q = SplineQuadrature(PeriodicBSplineBasis(UniformMesh(n, 2π), p)) + M32 = Float32.(mass_matrix(q)) + op = CirculantMass(M32, n) + @test op isa CirculantMass{Float32} + @test eltype(op) === Float32 + + b32 = randn(Float32, n) + @test norm(M32 * (op \ b32) .- b32, Inf) < 1.0f-3 + + # and the deflation works there too + S32 = Float32.(stiffness_matrix(q)) + opp = CirculantMass(S32, n; kernel = :project) + rhs = randn(Float32, n) + rhs .-= sum(rhs) / n + x = opp \ rhs + @test norm(S32 * x .- rhs, Inf) < 1.0f-3 + @test abs(sum(x)) < 1.0f-3 * n + end + end + + @testset "$(rpad("the deflated construct and solve infer and do not allocate",76))" begin + n = 32 + bu = PeriodicBSplineBasis(UniformMesh(n, 2π), 3) + bg = PeriodicBSplineBasis(GradedMesh(n, 2π), 3) + Su = stiffness_matrix(SplineQuadrature(bu)) + Sg = stiffness_matrix(SplineQuadrature(bg)) + + # the keyword is a literal at the call site, so it must constant-propagate into a + # concrete type -- a two-way Union here propagates into every consumer downstream. + # Written as a function rather than a closure: a closure over a local infers as Any + # on the 1.10 compat floor. + project(M, b) = mass_operator(M, b; kernel = :project) + @test @inferred(project(Su, bu)) isa CirculantMass + @test @inferred(project(Sg, bg)) isa FactorizedMass + + # the circulant deflation is a stored zero factor, not a branch, so the solve + # allocates no more than the undeflated one does: nothing + op = project(Su, bu) + x = randn(n) + y = similar(x) + mass_solve!(y, op, x) + @test (@allocated mass_solve!(y, op, x)) == 0 + end + @testset "$(rpad("in-place solve, aliasing, and no allocation",76))" begin q = SplineQuadrature(PeriodicBSplineBasis(UniformMesh(32, 2π), 3)) op = mass_operator(q) From 7789deb0902c010ff85e8ac7f40f90848e779769 Mon Sep 17 00:00:00 2001 From: Michael Kraus Date: Mon, 14 Sep 2026 16:35:03 +0200 Subject: [PATCH 4/4] Give the kernel gate a measured decade of headroom MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The archived sweep is wider than the one the bound was first set from: it reaches the degenerate corner, n barely above 2p+1, where every row of the stiffness matrix is nearly full. There the rounding in `M𝟙` saturates half of `N * eps(T) * ‖M‖∞` rather than a sixth of it, leaving a factor of two between a legitimate singular assembly and refusal. The gate now sits a decade above that scale. It costs nothing, because the two sides are not in competition: the smallest violation an invertible assembly produces is 1e4 (single) to 1e12 (double) above the scale, so the margins become 25x and 57x below, 6.6e11 and 1.2e3 above. `scripts/mass_tolerance_margins.jl` measures both sides of all three tolerances, and is where the figures in the comments come from. Co-Authored-By: Claude Opus 5 (1M context) --- CHANGELOG.md | 5 +- scripts/mass_tolerance_margins.jl | 138 ++++++++++++++++++++++++++++++ src/mass.jl | 26 +++--- 3 files changed, 156 insertions(+), 13 deletions(-) create mode 100644 scripts/mass_tolerance_margins.jl diff --git a/CHANGELOG.md b/CHANGELOG.md index 246b13e..5b9377c 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -49,8 +49,9 @@ of them were absolute, and each accepted something it should have refused: 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 `N * eps(T)` times `‖M‖∞`, with no floor, which is the - standard the circulant path already used. + 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 diff --git a/scripts/mass_tolerance_margins.jl b/scripts/mass_tolerance_margins.jl new file mode 100644 index 0000000..f8546de --- /dev/null +++ b/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() diff --git a/src/mass.jl b/src/mass.jl index 518d81e..05d0c79 100644 --- a/src/mass.jl +++ b/src/mass.jl @@ -105,7 +105,7 @@ end # positive definite to working precision exactly when the smallest is above `n * eps` of the # largest -- the same relative standard `_check_constant_kernel` and `_reciprocal_eigenvalues` # hold their own assertions to. A matrix whose kernel is `span{𝟙, v}` puts that ratio at 2e-16, -# where a genuine assembly holds it above 0.19 over degrees 1 to 4 and 16 to 256 cells. +# where a genuine assembly holds it above 0.18 -- `scripts/mass_tolerance_margins.jl`. function _check_definite_minor(F) issuccess(F) || throw(_singular_beyond_the_constants()) d = diag(F) @@ -119,17 +119,20 @@ end # the same path, and both the positive-definite minor and the mean-free gauge would then be # answering a question nobody asked. # -# The bound is `n * eps(T)` relative to ‖M‖∞, which is the row-sum scale `M𝟙` is measured on, -# and it is the standard `_reciprocal_eigenvalues` holds the circulant path to. It has to be -# *relative* and it has to have no floor: a bound that cannot fall below one is absolute for +# The scale is `n * eps(T)` relative to ‖M‖∞, the row-sum scale `M𝟙` is measured on. It has to +# be *relative* and it has to have no floor: a bound that cannot fall below one is absolute for # every matrix smaller than that, and accepts any invertible matrix scaled down far enough. -# `M𝟙` of a genuine singular assembly is pure rounding and grows with n like the bound does -- -# over degrees 1 to 8 and 64 to 4096 cells it reaches 0.16 of it, where an invertible mass -# matrix sits 10¹³ above it in the same units. +# +# `M𝟙` of a genuine singular assembly is rounding, and it saturates about half of that scale in +# the worst corner -- n barely above 2p+1, where every row is nearly full. The gate therefore +# sits a decade above the scale rather than on it. It can afford to: the two sides are not in +# competition here, since the smallest violation an invertible assembly produces is 10⁴ (single +# precision) to 10¹² (double) above the scale. `scripts/mass_tolerance_margins.jl` measures +# both sides. function _check_constant_kernel(M::AbstractMatrix{T}) where {T} n = size(M, 2) residual = norm(M * ones(T, n), Inf) - tol = n * eps(T) * norm(M, Inf) + tol = 10 * n * eps(T) * norm(M, Inf) residual ≤ tol || throw(ArgumentError( "kernel = :project deflates the constants, but M𝟙 has norm $(residual) against a " * "tolerance of $(tol); the kernel of this matrix is not what the projection assumes")) @@ -308,9 +311,10 @@ end # `rtol` is relative to the largest entry of the first column, so the verdict does not change # when the assembly is scaled. It defaults to `sqrt(eps(T))`, which is what makes the question # answerable in every element type: the residual of a genuinely circulant assembly is -# rounding, at most ~200·eps of that scale, while a graded mesh misses by 6e-2 of it. An -# absolute default instead fixes the answer to one element type -- at 1e-10 no `Float32` -# assembly is circulant at all, since `Float32` rounding alone exceeds it. +# rounding, at most 180·eps of that scale, while a graded mesh misses by 6e-2 of it +# (`scripts/mass_tolerance_margins.jl`). An absolute default instead fixes the answer to one +# element type -- at 1e-10 no `Float32` assembly is circulant at all, since `Float32` rounding +# alone exceeds it. # # The generic check: every entry against the first column. Quadratic, which for a dense # matrix is what the question costs.