A Laplacian is singular, so every solve with one has to remove its kernel first. There are three ways, they agree to fourteen digits, and the one everybody uses carries a free parameter that no account of the method mentions and that moves the condition number by nine hundred.
Every Laplacian is singular. That is not a defect of a particular graph or an artefact of a
discretisation; it is the definition. The rows sum to zero, so the vector of ones is in the kernel,
and a system Lx = b has either no solution or a whole line of them.
The line is not an inconvenience — it is the physics, and the kernel is known
exactly rather than computed. L is a network of conductances, b is a set of
injected currents, and x is a potential. Potentials are defined up to a constant, and a solver being
asked to return one has to be told which constant. Nothing in the matrix says.
Fig. 1 Thirty ways to remove the kernel of one system, by deleting one row and column, plotted against the
degree of the vertex removed. The best and worst are nine hundred times apart.
So a Laplacian solve has a preliminary step that a nonsingular solve does not, and there are three
standard ways to take it. This essay measures all three against an answer computed a fourth way, and
then measures the parameter that the most popular of them hides.
Ground a vertex. Pick one, delete its row and column — the same construction every effective
resistance in this field is computed through — solve the remaining (n − 1) × (n − 1)
system, and set the deleted potential to zero. This is the electrical reading — earthing a node — and
it is the standard move because it is free: the reduced matrix is symmetric positive definite for
any connected graph, so a Cholesky factorisation applies directly.
Project. Solve on the subspace orthogonal to the kernel. In practice this means subtracting the
mean from b before starting and from x afterwards, and running an iterative method that stays in the
subspace because the matrix maps it to itself. Conjugate gradients does this without being asked.
Regularise. Solve (L + τ·11ᵀ/n)x = b instead, with b orthogonal to the ones. That rank-one
addition moves the eigenvalue at zero to τ and leaves every other eigenvalue and eigenvector alone,
so for b in the range of L the answer is the same one — and the matrix is positive definite, so
again a Cholesky applies.
All three are exact in exact arithmetic and none is obviously better. The interesting question is
what they cost in the arithmetic there actually is.
Before comparing the three, it is worth asking what a solver does when handed the singular matrix
directly, because the answer is not what a reader of the paragraph above would predict.
The site’s own pivoted elimination, given L and a right-hand side that sums to zero:
graph smallest pivot flagged singular what solve returned
path 30 0 yes nothing — it refused
grid 6×7 6.7·10⁻¹⁶ no a vector
cycle 24 6.7·10⁻¹⁶ no a vector
star 30 8.7·10⁻¹⁵ no a vector
preferential 1.2·10⁻¹⁴ no a vector
complete 20 3.6·10⁻¹⁴ no a vector
Eight of the nine graphs tried do not trip the singularity test, and the one that does is the
path — the graph whose elimination happens to be exact, because every pivot along a path is a ratio
of small integers and every one of them is representable. Everywhere else the last pivot is the
accumulated rounding of a quantity that is mathematically zero, it comes out at 10⁻¹⁵ or so, it is
not zero, and nothing stops.
What comes back is not garbage. On the grid the returned vector has a relative residual of
7.4·10⁻¹⁶, and after subtracting its mean it agrees with the pseudoinverse solution to 2.7·10⁻¹⁵ —
it is a perfectly good point on the solution line. Which point is decided by what the last division
by 6.7·10⁻¹⁶ did, and here it shifted the whole answer by a constant of 2.54.
So the failure mode is the one this site keeps finding: not a crash, not a warning, but an answer
of the right shape whose arbitrary component was chosen by a rounding error. Two runs of the same
program on two machines that round that last pivot differently return potentials differing by a
constant, and every residual either of them reports is 10⁻¹⁶. Deciding that a zero has
arrived is the essay about the test that did not fire,
and this is that test failing on a matrix that is exactly singular by construction.
They agree, and the agreement is the uninteresting half #
Against a fourth route — the pseudoinverse assembled from the eigendecomposition, summing over the
n − 1 non-zero modes — the three routes give:
Every route is at the rounding level, and the three are within a factor of two of each other on
every graph. The forward errors track κ⁺(L) — the condition number of the matrix restricted to the
subspace it is inverted on — exactly as the error identity
says they should: the path has ten times the condition number and ten times the error.
That is a result and it is worth stating before the rest of the essay complicates it: the choice
of how to remove the kernel does not, on these graphs, cost accuracy. A reader who takes only that
away has taken away something true.
Grounding requires choosing a vertex, and no account of the method says which. In practice it is
whichever index the code drops — the last one, or the first, or whichever the data structure makes
convenient. It reads like a labelling convention.
It is a parameter, and on a graph with unequal degrees it is worth three orders of magnitude.
Grounding vertex v leaves the matrix L with row and column v deleted, and that matrix’s condition
number depends on v. Measured over every choice, on the same graph:
The star’s entry is the one to hold on to. Grounding the hub of a thirty-vertex star leaves the
identity matrix — every remaining vertex is joined to the ground and to nothing else, so the reduced
Laplacian is diagonal with every entry 1, and its condition number is exactly 1. Grounding a leaf
leaves a matrix with condition number 898. Two ways of solving the same system, one of which is
perfectly conditioned and one of which is not, differing only in which index was deleted.
The pattern across the table is that a high-degree vertex is the better ground, and it has a
reading rather than being a coincidence. Grounding fixes a potential, and fixing the potential of
something the rest of the graph is strongly attached to constrains more of the graph. Grounding a
leaf fixes the potential of a vertex whose only connection is one edge, which constrains almost
nothing, and the remaining system is left with a nearly-free direction.
The size of the effect is governed by how unequal the degrees are. On the barbell, where every
vertex has degree nine or ten, the spread is 1.29 and the choice is genuinely irrelevant. On the
star, where the degrees are 29 and 1, it is 898.
Fig. 2 The same measurement on a preferential-attachment graph, where the degrees run from two to
seventeen and the spread is a factor of ten.
Three orders of condition number sounds like three orders of accuracy, and it is not, which is the
part worth being careful about.
The forward error of a backward-stable solve is bounded by κ times the unit roundoff, so on the star
the badly grounded system has a bound 898 times looser. Both bounds are around 10⁻¹⁶ times something
small, so both solves return an answer good to thirteen or fourteen digits and nobody notices.
Where it does bite is an iterative solve. Conjugate gradients takes a number of steps that grows
like the square root of the condition number, so a factor of 898 in κ is a factor of thirty in the
iteration count — and unlike the accuracy, the cost is not sitting nineteen orders below anything
anybody cares about. The rate the condition number predicts
is the essay about that relationship; here the point is that a labelling convention has walked into
it.
And it bites in a second place, which is where this field goes next. A tree
preconditioner is scored by a combinatorial quantity, and
the scoring assumes a particular grounding; two implementations that ground differently are
comparing preconditioners for two different matrices.
The one place the three routes genuinely part company #
Everything above says the three are equivalent, and on a consistent right-hand side they are. The
right-hand side is not always consistent.
A Laplacian system is solvable only when b sums to zero — the injected currents have to balance. In
theory they do. In practice b is assembled from measurements, or from a discretisation, or from a
previous computation, and its sum is not zero but ε. Nothing in the setup guarantees otherwise, and
nothing reports it.
Adding a component ε along the vector of ones to an otherwise exact right-hand side, and measuring
each route’s answer against the pseudoinverse solution of the consistent part:
Two of the columns do not move at all and one of them tracks ε, amplified by ten.
The reason is what each route does with the inconsistent part. Projecting removes it explicitly —
the mean is subtracted before the iteration starts, so the component never enters. Regularising sends
it into the direction the rank-one term owns, where it contributes a multiple of the constant vector
that is removed when the answer is centred. Grounding does neither: the reduced system is
nonsingular, so it has a unique solution for every right-hand side, including inconsistent ones,
and it returns the exact answer to a problem in which the missing current was silently injected at
the grounded vertex.
That is a defensible thing to do and it is not what anybody thinks they asked for. It is also
invisible: the grounded solve’s residual on its own reduced system is 10⁻¹⁶, so every check the
solver can run reports success. The inconsistency is a property of the data that only the full
system knows about, and grounding is precisely the step that threw away the row that knew.
The practical rule is one line: if the right-hand side is not exactly balanced, balance it before
grounding. Subtracting the mean costs one pass and removes the entire effect. What it does not do
is tell anybody the data was inconsistent, which — given that the amplification is ten and the
inconsistency is often a bug rather than a rounding — is the thing actually worth reporting.
The projection route has a property the other two do not, and it comes straight from the
exactness of the kernel.
The kernel of the combinatorial Laplacian is spanned by the vector of ones, and L annihilates it
exactly: measured, ‖L·1‖∞ is 0, not 10⁻¹⁶. So the subspace the projection works on is known in
advance, exactly, without computing anything. There is no eigenvector to estimate and no basis to
orthogonalise, and the deflation cannot drift, because the direction being deflated is not an
approximation to anything.
Compare that with the ordinary situation. Deflating a converged eigenvector
from a Krylov space means deflating a vector known to 10⁻¹⁴, and the leak is at that level and
compounds over the iteration. Here there is no leak of that kind at all.
And the normalised Laplacian does not have this. Its kernel is spanned by D^{1/2}·1, which costs
n square roots to build and is not exact; measured on four graphs, ‖𝓛 D^{1/2}1‖∞ comes to between
2.2·10⁻¹⁶ and 3.8·10⁻¹⁵. Nothing in this essay’s measurements fails because of that — the numbers
above are at the rounding level for both matrices — but the kind of statement has changed. For one
matrix the deflation is structurally exact. For the other it holds to the rounding level, which is a
claim about the arithmetic rather than about the algebra, and a claim about the arithmetic is one an
iteration can wear away.
Fig. 3 The star: one vertex of degree twenty-nine and twenty-nine of degree one, which is the graph the
grounding choice matters most on.
The third route deserves a sentence, because it shares a name with a whole field on this site and is
not the same thing.
Regularisation proper changes the answer. A Tikhonov parameter trades
residual against solution norm, the answer moves as the parameter does, and choosing it is the
computation rather than a preliminary to it. Nothing like that is happening here: adding τ·11ᵀ/n
moves one eigenvalue from zero to τ and leaves the other n − 1 eigenpairs untouched, so for a
right-hand side in the range of L — which is any b that sums to zero, and any physical b does — the
solution is unchanged for every τ.
What τ does change is the conditioning. Setting it far below λ₂ makes the shifted matrix badly
conditioned in a direction the answer does not live in; setting it near λmax makes the matrix
well conditioned and does nothing else. The right choice is therefore somewhere in the middle of the
spectrum and is not delicate, which is exactly what a parameter that does not affect the answer
should look like.
That distinction — a parameter that moves the answer against one that moves only the arithmetic — is
worth keeping, because it is the difference between a modelling decision and an implementation
detail, and this field is about to spend three essays on a case where the two are confused.
On these graphs the three routes are equivalent to fourteen digits. The choice is not an
accuracy question and should not be presented as one.
Grounding has a hidden parameter and it is worth setting. Ground the highest-degree vertex. It
costs one pass over the degree array, it never hurts, and on a graph whose degrees are unequal it
buys between a factor of three and a factor of nine hundred in the conditioning — which becomes a
factor of thirty in the iteration count of anything iterative.
The projection route is the only one whose correctness is structural rather than numerical, and
only for the combinatorial Laplacian, whose kernel is exact. That is not a reason to prefer it; it
is a reason to know which of the three is being run when an iteration is long.
Fig. 4 A grid, whose degrees run from two to four and whose spread is a factor of three.Fig. 5 A barbell, where every degree is nine or ten and the choice does not matter.Fig. 6 A path, the worst-conditioned graph here, where the spread is small and the conditioning is not.Fig. 7 Two blocks, where the degrees are similar and the spread is correspondingly narrow.Fig. 8 The path, whose conditioning comes from its length rather than from its degrees.Fig. 9 The spectrum whose zero the whole essay is about removing.Fig. 10 The bounds on the smallest non-zero eigenvalue, which is what κ⁺ is one over.Fig. 11 An iteration whose step count is the square root of the condition number this essay moves.Fig. 12 Changing a condition number on purpose, in the field that does it for a living.Fig. 13 And what that buys in iterations, which is where the grounding choice is felt.Fig. 14 A solve measured against an answer that is known rather than estimated.Fig. 15 The two errors, which the table of three routes reports the second of.Fig. 16 The amplifier, which the grounding choice moves by nine hundred without touching the problem.Fig. 17 The identity the agreement of the three routes is read against.Fig. 18 A backward error measured against a structure rather than against a matrix.Fig. 19 What a singular system does to a solver that does not know it is singular.Fig. 20 A condition estimate, and the spread of what it might have returned.Fig. 21 A diagonal scaling recovering accuracy, which is what the normalised Laplacian is.Fig. 22 Condition numbers that move with the units, which a grounding choice is a discrete version of.Fig. 23 Orthogonality lost in a Krylov space, which is the leak an exact deflation does not have.Fig. 24 A residual that says nothing about the error, in the field about what a residual is worth.Fig. 25 A tolerance asked for against one obtained, which every iterative route here is stopped by.Fig. 26 The mantissa on a slider, on a solve rather than on a graph.Fig. 27 A factorisation with its residual printed, which is the site’s rule.Fig. 28 A threshold deciding an integer, which is the shape a component count also has.Fig. 29 Forming an inverse rather than solving, which the pseudoinverse route here does on purpose.Fig. 30 A bound never attained, in the field that measured it first.Fig. 31 Two orderings of one algebra, which the three routes here are a third instance of.