[michakraus]

A MassOperator could not be built from a singular matrix at all, so a caller with a periodic
stiffness matrix had to make it invertible first. The only way to do that without changing the
answer is the rank-one shift by the mean projector, S + 𝟙𝟙ᵀ/N — and that is a scalar added to
every entry, so it turns an O(N) banded assembly into a structurally full one.

This branch lets the operator deflate the kernel instead.

What it adds

kernel = :project on mass_operator, CirculantMass and FactorizedMass. It states that the
kernel is the constants, and asks for the solution that has no constant component — the mean-free
solution, which is the one -φ'' = ρ determines on a periodic domain. kernel = :reject is the
default and is the previous behaviour exactly.

representation	how it deflates
CirculantMass	the constant Fourier mode gets a zero factor instead of an infinite one. That is the Moore–Penrose pseudoinverse of a circulant operator, and it is the rule PoissonSolverFFT already applies when it sets the k = 0 factor to zero.
FactorizedMass	the Cholesky factor is taken of the minor that drops the first degree of freedom, which is positive definite exactly when the kernel is the constants and nothing more. mass_solve! sets that degree of freedom to zero and takes the mean out afterwards.

Neither touches the matrix. A bounded basis raises rather than accepting the keyword; the banded
representation carries no deflation.

Both verify the assertion rather than trusting it. The circulant path checks that the constant
mode is the one that vanishes and that no other mode does; the factorised path checks M𝟙 ≈ 0.
A matrix singular for some other reason takes neither path — it is rejected with a message saying
so, rather than getting a plausible answer to a question nobody asked.

What the shift costs, measured

The periodic stiffness matrix at degree 4, on a UniformMesh:

N	stored entries, S	stored entries, S .+ inv(N)
32	288	1 024
128	1 152	16 384
512	4 608	262 144
2048	18 432	4 194 304

At N = 2048 the shifted matrix alone measures Base.summarysize 67 125 416 B. It stays a
SparseMatrixCSC, so it carries a row index beside each of the N² entries.

And S + 𝟙𝟙ᵀ/N is circulant like S itself, so its whole effect is visible in the transform.
It changes one Fourier coefficient, from zero to one:

n= 32  ‖ℱ(c + 𝟙/N) - ℱ(c) - e₀‖∞ = 2.78e-17
n=128  ‖ℱ(c + 𝟙/N) - ℱ(c) - e₀‖∞ = 5.20e-18
n=512  ‖ℱ(c + 𝟙/N) - ℱ(c) - e₀‖∞ = 3.47e-18

The O(N²) fill-in buys exactly that one number. The deflation sets the same coefficient without
forming anything. Over 32 to 1024 cells the two solves agree to 1e-15 relative, on both
representations, and the mean of the deflated result is zero to 1e-15.

The bug fix that came with it

The circulant invertibility check was abs(λ) > eps(T) — absolute, and unscaled. The zero
eigenvalue of a periodic stiffness matrix comes out of the transform at about 1e-14 rather than
at zero, so the check passed it and the solve then divided by it. What came back was amplified
noise rather than an error. The tolerance is now relative, at N * eps(T) times the largest
eigenvalue, where pinv puts the same cutoff.

This is what made kernel = :reject meaningless without the fix: :reject has to reject what
:project projects.

Breaking, and why it is worth it

• FactorizedMass gained a fourth type parameter, which carries the kernel mode. A signature
  spelled FactorizedMass{T, MT, FT} needs the extra parameter. The constructor is unchanged.
• A CirculantMass solve multiplies by precomputed reciprocal eigenvalues instead of dividing
  by the eigenvalues. It moves N ÷ 2 + 1 divisions out of every solve and into the construction,
  and mirrors PoissonSolverFFT's k⁻². Results can differ in the last ULP; no doctest moved.

Its relation to #10 and #11

The case #10 is about is the rank-one shift in the periodic B-spline Poisson solve, and after this
branch that caller has no shift to build: PoissonSolvers hands over the sparse assembly itself
(JuliaPlasma/PoissonSolvers.jl, branch kernel-projection). The whole solver at degree 4 on 2048
cells goes from 72 551 536 B to 5 426 160 B, and nothing is O(N²) anywhere along the path.

That leaves #10's literal request — a Circulant type and a first-column entry point, which is
#11 — as a smaller convenience for a caller that knows a circulant matrix in O(N) and is not
coming from a singular assembly. It is no longer what the reported case needs. Worth deciding on
its own merits, or closing.

Verification

• Pkg.test() on this branch: 64 394 passed, 0 failed, Aqua.test_all included
  (ambiguities, piracy, stale deps, unbound parameters, undefined exports all green).
• Docs build (julia --project=docs -e 'include("docs/make.jl")'): exit 0, doctests ran, no
  missing-docstring and no broken-@ref warnings. No doctest moved despite the reciprocal change.
• JuliaFormatter --check, fatou lint and the NFC check: clean on every changed file.

julia-branch-verifier ran the seven-point pass against the diff. What it found:

Fixed

• src/mass.jl:71 — the two-branch kernel constructor defeated constant propagation.
  FactorizedMass(M) and FactorizedMass(M; kernel = :project) both inferred as a two-way
  Union, and it propagated into SplineQuadrature's MO type parameter, which was concrete
  before this branch. Base.@constprop :aggressive restores concrete inference for a literal
  keyword; a runtime Symbol still gives the Union, which is correct. Verified on the 1.10
  [compat] floor as well.
• src/mass.jl:12, src/mass.jl:173 and docs/src/usage/quadrature.md:235 — three places
  described the circulant solve as "a pointwise division". The solve multiplies by stored
  reciprocals now, and under :project it multiplies the constant mode by zero.
• src/mass.jl:339 — a comment claimed the :project solve allocates one temporary. Measured in
  a fresh process at --check-bounds=auto: 1456 B against 1296 B for :reject, the difference
  being the mean-free right-hand side.
• CHANGELOG.md — the stored-entry pair it quoted, 11 264 entries and 188 584 B, is degree 5;
  at the degree 4 it names the figures are 9 216 and 155 816 B. Its two agreement figures were the
  best case rather than the range, and are now stated relative and over the range measured.

Raised and not accepted

• It reported 16 785 576 B for the shifted matrix as unreproducible, measuring 8 388 656 B.
  It had measured Matrix(S) .+ inv(N), a dense Matrix{Float64}. The expression the changelog
  names is S .+ inv(N), which stays a SparseMatrixCSC and so carries an Int64 row index
  beside every Float64: 1 048 576 stored entries at exactly 16 B each. The figure stands.

Checked and clean

• Correctness of the :project derivation: with M symmetric and M𝟙 = 0, solving the minor
  and setting y[1] = 0 reproduces the dropped row automatically, since 𝟙ᵀMy = (M𝟙)ᵀy = 0.
• The y === x aliasing contract, verified rather than reasoned: 5 trials × both representations
  × mass_solve!, ldiv! and the column-view path of \(op, ::AbstractMatrix), all agreeing
  with the non-aliased result to 1e-12, residual ‖Sx − b‖∞ ≈ 2e-15.
• Both new tolerance gates at N up to 1024 and degree up to 4: margins of 10⁵–10⁷ between
  the kernel mode and the cutoff, and between the cutoff and the next mode. Float32 builds on
  both paths.
• Type piracy: every added method is on an owned type or an owned name. No renamed field has a
  reader outside src/mass.jl, and the new fourth type parameter is spelled out nowhere.
• Duplicated dependency functionality: pinv and nullspace are dense and O(N³), CHOLMOD has
  no rank-deficient path, and FFTW has no deflation concept.
• Allocations: fatou lint clean, and the CirculantMass solve holds its asserted 0 bytes.
• Comments in the present tense throughout, with no history outside CHANGELOG.md.

Pre-existing, left alone

• The rejected rank-one shift is explained three times in src/mass.jl. Collapsing that to one
  statement and two cross-references is an editorial call rather than part of this change.

Not checked

• The sandbox blocks ps, so the machine was assumed idle for the allocation figures.
• Full-package inference on the 1.10 min job — only the construct was verified there. No test
  asserts inference, so a difference could not turn CI red.

🤖 Generated with Claude Code

[Copilot]

🟡 Changes recommended

Unresolved moderate issues affect scale-relative validation and offset-indexed inputs; additional null-mode tests are also needed.

Get a fresh assessment by requesting another Copilot review.

Pull request overview

Adds opt-in kernel projection for singular periodic assemblies while preserving rejection by default.

Changes:

• Implements circulant and factorized deflation.
• Adds singularity validation and reciprocal eigenvalue caching.
• Expands tests, documentation, and changelog entries.

File summaries

File	Summary
src/mass.jl	Implements projection and updated circulant solves. Findings: moderate issues with relative tolerance, offset-indexed matrices, and offset-axis vectors; missing nonconstant-null-mode coverage is a nit.
test/mass_tests.jl	Adds regression and validation tests; lacks cases for additional null modes in both representations (nit).
docs/src/usage/quadrature.md	Documents singular assembly usage.
CHANGELOG.md	Records the feature, fixes, and breaking changes.

Review details

Suppressed comments (4)

src/mass.jl:101

• This check is not scale-relative for matrices with norm below one: an invertible matrix such as αI with α < √eps(T) has M𝟙 below the absolute floor, so its minor factorizes and kernel = :project silently treats a nonsingular matrix as having a constant kernel. It also makes max invalid when T is complex. Scale the tolerance by norm(M, Inf) itself so the kernel assertion is relative.

    residual ≤ sqrt(eps(T)) * max(one(T), norm(M, Inf)) || throw(ArgumentError(

src/mass.jl:81

• M is declared as an AbstractMatrix, but M[2:end, 2:end] assumes one-based axes and does not necessarily remove exactly one row and column. With an offset-indexed matrix this selects the wrong minor (or changes its size), so kernel = :project can reject or factor the wrong system; derive the dropped ranges from the matrix axes instead.

            F = cholesky(Symmetric(M[2:end, 2:end]); check = false)

src/mass.jl:256

• The new guard for a second, nonconstant null mode is not exercised by the added mass tests: the failure cases cover only M𝟙 ≠ 0, while the successful circulant cases have all other modes nonzero. Add a small circulant matrix with a zero nonconstant Fourier mode (for example [[1, 1], [1, 1]]) and assert that kernel = :project throws, so this safety check is protected against regression.

        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"))

test/mass_tests.jl:130

• These assertions only exercise the M𝟙 ≠ 0 guard. They do not cover the other promised rejection: a matrix with M𝟙 = 0 but an additional null mode, which should hit the circulant second-mode check and the factorized minor failure. Add at least one such case for each representation so a future change cannot silently accept a larger kernel.

        # 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)

• Files reviewed: 4/4 changed files
• Comments generated: 1
• Review effort level: Lite

---

💡 Add a code-review agent skill or configure MCP servers for context-aware, tailored reviews. Learn more in the docs.

[michakraus]

Independent review — kernel-projection

Reviewed from a detached worktree at f49f39a, not from the PR body. Where I repeat one of the
julia-branch-verifier claims below it is because I re-measured it myself; where I disagree I say
so. Every number here comes from a probe I ran in a fresh process on this branch.

Verdict: request changes. One line — the :project gate in _check_constant_kernel — is not
scale invariant, and under a rescaling it accepts a matrix whose kernel is not the constants and
then returns a silently wrong answer. That is the same defect class as the absolute-tolerance bug
this branch sets out to fix, on the other representation. Everything else below is advisory.

CI: green, and genuinely so. All eleven check runs on f49f39a report success, including
Julia nightly - ubuntu-latest and all three Julia min (1.10) platforms — so nightly is not
being waved through by a continue-on-error. I also ran Pkg.test() locally on the branch: passed,
Aqua.test_all included (piracy, ambiguities, stale deps, compat bounds, persistent tasks).

Findings

#	severity	file:line	claim	evidence
1	correctness-risk	src/mass.jl:101	the constant-kernel gate is not scale invariant. max(one(T), norm(M, Inf)) floors the bound at 1, so for ‖M‖∞ < 1 the test is absolute at sqrt(eps(T)) ≈ 1.5e-8. An invertible, well-conditioned matrix scaled below that is accepted by kernel = :project and the solve returns a wrong vector with no error.	below
2	quality	src/mass.jl:101 vs src/mass.jl:240	the two paths verify the same assertion to standards ~8 orders apart — sqrt(eps) against n·eps·max|λ̂|. The same assembly therefore gets a strict check on a UniformMesh and a lax one on a GradedMesh.	below
3	quality	src/mass.jl:274, src/mass.jl:295	pre-existing and adjacent: _check_circulant's atol = 1e-10 is the other absolute tolerance in the same constructor, and it makes CirculantMass unreachable in Float32 at every degree. The PR body's "Float32 builds on both paths" holds only for a hand-converted matrix.	below
4	quality	test/mass_tests.jl:96-133	neither "singular beyond the constants" branch is exercised by the suite — not _reciprocal_eigenvalues' second-vanishing-mode error, not FactorizedMass' minor-Cholesky error. Both are correct; they are just untested.	below
5	nit	src/mass.jl:191	CirculantMass carries no kernel-mode type parameter while FactorizedMass gained one, so mass_operator(…; kernel = :project) records the choice in the type on one mesh and not on the other. A consumer cannot dispatch on it uniformly.	read
6	nit	src/mass.jl:405	an unrecognised kernel name on a bounded basis reports kernel = :ignore is implemented for a periodic basis only, which says :ignore is valid somewhere. The periodic path gets the right message.	below
7	nit	CHANGELOG.md:31	16 785 576 B is right, but it is the size of S .+ inv(N), an expression that appears nowhere in the repository. The docs and the new test both write Matrix(S) .+ inv(N), which is 8 388 656 B.	below
8	nit	test/mass_tests.jl	no @inferred and no @allocated assertion on the new :project paths, so neither the @constprop fix nor the 0-byte circulant :project solve has a regression guard.	read

Finding 1 — the gate fails open under rescaling

mass_matrix of PeriodicBSplineBasis(GradedMesh(32, 2π), 3), scaled by s. This matrix is SPD
and invertible; M𝟙 is nowhere near zero; kernel = :project must refuse it at every s.

s=1e+00     ||Ms*1||=2.193e-01  gate=1.490e-08  raises
s=1e-04     ||Ms*1||=2.193e-05  gate=1.490e-08  raises
s=1e-08     ||Ms*1||=2.193e-09  gate=1.490e-08  ACCEPTED  ||Ms*y - x̄||inf=7.74e-01  rel=4.57e-01
s=1e-10     ||Ms*1||=2.193e-11  gate=1.490e-08  ACCEPTED  ||Ms*y - x̄||inf=3.61e+00  rel=1.76e+00
s=1e-12     ||Ms*1||=2.193e-13  gate=1.490e-08  ACCEPTED  ||Ms*y - x̄||inf=1.25e+00  rel=6.36e-01

The relative residual of the returned solution is 0.46, 1.8, 0.64 — not a degraded answer, a
wrong one, and no exception. The circulant path refuses the same matrix at every scale, because its
gate is purely relative:

circulant, s = 1e+00, 1e-08, 1e-12  ->  all three raise
  "kernel = :project deflates the constants, but …"

A mass matrix on a domain of length L with N cells has entries ~L/N, so ‖M‖∞ < 1 is already
the normal case — at N = 32, L = 2π above, ‖M‖∞ = 0.105. The floor is live today; only the
1e-8 margin keeps it from biting. A micrometre-scale domain, or any non-dimensionalisation that
shrinks the assembly, reaches it.

The fix is the same one the branch already applies on the circulant side: make the bound relative
and drop the max(one(T), …) floor, e.g. norm(M*𝟙, Inf) ≤ tol * norm(M, Inf). A zero matrix is
the only case the floor protects, and that one is better handled explicitly.

Finding 2 — the two gates disagree by eight orders

At n = 64, degree 3, max|λ̂| ≈ 0.098: the circulant gate is n·eps·max|λ̂| = 1.4e-15, the
factorised gate is sqrt(eps)·max(1, ‖M‖∞) = 1.0e-7. Perturbing a stiffness matrix symmetrically
so that M𝟙 grows:

||P*1||inf=1.359e-11  gate=1.012e-07  ACCEPTED  ||P y - rhs||inf = 2.009e-11
||P*1||inf=1.358e-09  gate=1.012e-07  ACCEPTED  ||P y - rhs||inf = 2.739e-09
||P*1||inf=1.358e-08  gate=1.012e-07  ACCEPTED  ||P y - rhs||inf = 4.750e-09
||P*1||inf=1.358e-07  gate=1.012e-07  raises

So the residual the :project solve delivers degrades in proportion to whatever the gate lets past,
and the gate lets 1e8 times more past than its sibling does. If the loose bound is deliberate —
to absorb the quadrature's own error in the assembly — that deserves a sentence at
src/mass.jl:95-98, because right now the two checks read as one policy and behave as two.

Finding 3 — the other absolute tolerance, in the same constructor

periodic uniform Float32 p=1..4, n=16,32,64
  -> all 12 RAISE: "the mass matrix is not circulant to within 1.0e-10"

Every Float32 periodic-uniform SplineQuadrature fails at construction, for every degree, on the
pre-existing atol = 1e-10 default. Pre-existing, yes — but it is eight lines from the tolerance
this branch fixes, it is the same mistake, and it means the Float32 evidence in the PR body was
necessarily gathered from a hand-converted matrix rather than from any assembly this package can
produce. Scaling it as atol * max(one(T), maximum(abs, c)) would need atol to default to
something T-dependent; that is a small change and belongs here rather than in a follow-up nobody
will open.

Finding 4 — the two rejection branches nobody tests

Both fire correctly when I build the inputs by hand, so this is coverage, not a defect:

E1  circulant with the k=0 and k=n/2 modes both zero
      kernel=:reject  -> "…has a zero eigenvalue and is not invertible…"
      kernel=:project -> "…singular beyond the constants: a second mode of the transform vanishes…"
E2  PSD with kernel = span{𝟙, v}, v ⟂ 𝟙
      FactorizedMass  -> "…singular beyond the constants: the minor that drops the first…"

sparse(Diagonal([1.0, 0.0, 1.0, 1.0])) at test/mass_tests.jl:131 exercises
_check_constant_kernel only; it never reaches the minor. Two more @test_throws would close this.

Finding 6 — the message

bounded basis,  kernel = :ignore -> "kernel = :ignore is implemented for a periodic basis only;
                                     a bounded basis takes the banded representation…"
periodic basis, kernel = :ignore -> "kernel must be :reject or :project, not :ignore"

Finding 7 — the changelog figure

Reproduced exactly, degree 4, N = 1024:

nnz(S)=9216   summarysize(S)=155816
S .+ inv(N)           SparseMatrixCSC   nnz=1048576  summarysize=16785576
Matrix(S) .+ inv(N)   Matrix                         summarysize= 8388656

The changelog's number stands, and its parenthetical about the row index is correct — the
branch-verifier's "unreproducible" was the dense expression, and the PR body's rebuttal is right.
But the only spelling a reader will find in the repository is Matrix(S) .+ inv(N), in
docs/src/usage/quadrature.md and in the new test, which is half that. Naming the expression in the
changelog removes the trap.

Verified and fine

• The FactorizedMass derivation is correct, and I checked the dropped row directly. With M
  symmetric and M𝟙 = 0, 𝟙ᵀMy = (M𝟙)ᵀy = 0, so row 1 of the residual is minus the sum of the
  others, which the minor has already matched. Numerically, on a GradedMesh, degrees 1–4,
  N ∈ {16, 33, 64}: full residual ‖Sy − x̄‖∞ ≤ 6.5e-14 and the row-1 residual alone
  ≤ 6.5e-14, i.e. at roundoff, never worse than the rest. The reconstruction is real, not assumed.
• Both :project paths are the Moore–Penrose pseudoinverse, which is a stronger statement than
  the branch's own test against the rank-one shift. ‖op \ x − pinv(Matrix(S))·x‖∞, over
  {Uniform, Graded} × p ∈ {2,4} × N ∈ {16,32}: 1.7e-15 to 7.2e-14.
• The circulant kernel check cannot be fooled about which mode vanishes. ĉ[1] is the
  constant mode; rfft returns indices 0 … ⌊n/2⌋ and the omitted half is the conjugate of the
  included one, so all(abs > tol, ĉ[2:end]) really does cover the whole spectrum, for odd n too.
• The factorised "nothing else vanishes" argument is sound, and it is stronger than the docstring
  claims. A PD minor plus M𝟙 = 0 forces M PSD with kernel exactly span{𝟙}: any
  w = ṽ + α𝟙 has wᵀMw = ṽᵀMṽ > 0 unless ṽ = 0. So the branch needs no separate PSD check. I
  also fed it a non-symmetric matrix with zero row sums but Symmetric(A)𝟙 ≠ 0 — it is rejected by
  the minor, not accepted.
• The n · eps(T) · max|λ̂| tolerance is correctly scaled and has enormous headroom in
  Float64. |λ̂min| / tol, uniform periodic, over p = 1…8 and n ∈ {32, 256, 2048}: worst
  case 1.30e+09 at p = 8, n = 2048. It is also pinv's own rtol = min(size)·eps, as the
  comment says. Two caveats, neither blocking: in Float32 the worst margin falls to 2.4 at
  p = 8, n = 2048, so around n ≈ 5000 a legitimate matrix would be called singular; and n·eps
  grows a factor n / log n faster than the FFT's actual accumulated error, so the bound is loose
  in the direction that costs, not the one that is unsafe.
• Aliasing holds on every path, measured rather than argued. y === x for
  mass_solve!, ldiv!(op, x) and the column-view path of \(op, ::AbstractMatrix), on both
  :project representations: difference from the non-aliased result exactly 0.00e+00 in all
  six. The :project FactorizedMass is safe because x[2:n] .- sum(x)/n materialises and the
  CHOLMOD solve completes before y[1] is written.
• Allocations, fresh process, --check-bounds=auto, n = 64, warmed, minimum of 20:
  circulant :reject 0 B, circulant :project 0 B, factorized :reject 3984 B,
  factorized :project 4560 B. The :project circulant solve really does stay allocation-free, and
  the 576 B difference on the factorised path matches the one extra length-63 temporary the
  comment at src/mass.jl:338 describes. (My absolute figures differ from the body's 1456/1296
  because I measured a different n; the difference is the claim, and it holds.)
• Inference is concrete on both 1.13 and the 1.10 [compat] floor. With a literal keyword,
  FactorizedMass(M; kernel = :project) infers
  FactorizedMass{Float64, SparseMatrixCSC{Float64,Int64}, CHOLMOD.Factor{Float64,Int64}, :project},
  and SplineQuadrature(b) carries a fully concrete MO parameter on both meshes. A runtime
  Symbol gives the two-way Union, which is the correct answer. The Base.@constprop :aggressive
  fix does what the body says it does, on 1.10 as well — I re-ran it there.
• Bounded basis: mass_operator(M, ::AbstractBSplineBasis; kernel = :project) raises, and
  PeriodicBSplineBasis <: AbstractBSplineBasis so the periodic method is the more specific one —
  the dispatch is not accidental. BandedMass takes no kernel keyword at all, which is the right
  way to carry "no deflation here".
• Pkg.test() on the branch: passed, Aqua.test_all green (piracy, ambiguities, stale deps,
  unbound parameters, undefined exports, compat bounds, persistent tasks). fatou lint src/mass.jl:
  clean. No new dependency, so [compat] needs nothing. Changelog entry present, breaking changes
  named, Unreleased — targeting 0.2.0 is the right bump for a 0.x breaking change. Comments are
  present tense throughout the new code with no history outside CHANGELOG.md.

Not checked

• The docs build. CI's Documentation and Doctests jobs are both green on f49f39a, so I took
  those rather than rebuilding.
• The downstream PoissonSolvers.jl kernel-projection branch and the 72 551 536 → 5 426 160 B
  figure, which lives in another repository.
• ExplicitImports. The new code adds no import; norm comes from the module's pre-existing
  using LinearAlgebra.
• Behaviour for n = 1, and for y/x of different lengths in the :project solve, which reads
  n = length(x). Both degenerate.