The basis nobody chose on purpose

The null-space method turns a constrained problem of size n + m into an unconstrained one of size n − m: write x = x_p + Zv with Ax_p = g and AZ = 0, and solve the reduced system ZᵀHZ. The essay that set it up measured its condition number as 21.13 and observed that it did not move while κ(A) crossed four decades.

Nothing in that derivation says which Z, and the set of legal answers is enormous. Any Z whose columns are independent and annihilated by A will do, every one gives the same x in exact arithmetic, and their condition numbers are decades apart.

The square that was quietly set to one #

The identity is

κ(ZᵀHZ) ≤ κ(H) · κ(Z)²

with the square arriving for the reason it arrives in the normal equations: ZᵀHZ is a Gram matrix in the H-inner product, and a Gram matrix’s condition number is the square of its factor’s.

The previous essay’s flat line was that identity with κ(Z) = 1, because the basis it used came from a QR of Aᵀ and had orthonormal columns by construction. So the flat line was true and it was also a choice, and the choice was made inside a routine rather than by anybody looking at the problem.

Set κ(Z) to 10⁸ and the bound says 10¹⁶. The measurement says 3.80·10¹⁶ against a κ(Z)² of 3.96·10¹⁶ — the square attained to two digits, on a matrix H whose own condition number is 100. The bound is not a bound here; it is the answer.

The two bases, and what each is for #

Orthonormal. Take a full QR of Aᵀ. The first m columns of Q span the row space of A and the last n − m span its null space, with κ(Z) = 1 exactly. That is the basis with no numerical question attached, and its price is density: Z is n × (n − m) with essentially every entry nonzero — measured at 100 per cent — even when A was sparse. On a problem where A is a few constraint rows across a large sparse stiffness matrix, Z is the largest object in the computation.

Fundamental. Split the columns of A into m basic and n − m nonbasic, so that A = [A_B A_N] after a permutation, and write

Z = [ −A_B⁻¹ A_N ] [ I ]

with the rows put back in order. The identity block is free, the other block is as sparse as A_B⁻¹A_N, and there is no orthogonalisation anywhere. Measured density: 50 per cent, half the orthonormal one’s.

The second is what a large-scale optimisation code uses, and it is what a simplex-based method produces for nothing, because the basic/nonbasic split is the simplex basis. Its condition number is a property of A_B, and A_B is chosen by whatever rule the code was using for other reasons.

Where the experiment had to be built rather than found #

The measurement needs a matrix on which the naive choice is bad and the matrix is fine, because otherwise the two explanations — a hard problem, and a badly made basis — cannot be separated.

A constraint matrix built as UΣVᵀ from random orthogonal factors cannot supply one. Its columns come from a random orthogonal V, so every subset of m of them is about as well conditioned as every other, and the naive rule is accidentally safe: measured κ(A_B) of 7.97·10⁵ for the first m columns against 1.02·10⁶ for the pivoted choice, which is no separation at all.

So the family is constructed: A = [B | N] with κ(B) = γ chosen and N well conditioned. The variable ordering is the problem’s rather than the solver’s, and a problem that lists its near-duplicate variables first is not exotic — it is what a model built up in blocks looks like, and it is what a mesh numbered by subdomain looks like.

That construction is the honest part of the page. The finding is not that the naive rule is always bad; it is that nothing stops it being arbitrarily bad, and that when it is, nothing in the problem’s own condition numbers says so.

The repair is pivoting, not orthogonalisation #

The obvious fix is to orthogonalise the fundamental basis, which costs a QR of an n × (n − m) matrix and gives back the density that was the reason for using it.

The cheaper fix is to choose the basic columns differently. Run a column-pivoted QR of A and take the m columns it selects: that is the greedy rule that maximises the residual norm at each step, it costs O(m²n), and it is what the rank-revealing factorisation already does for a different purpose.

Measured on the same problem at γ = 10⁸:

rule κ(Z) κ(ZᵀHZ) density error orthonormal 1.00 25.6 1.00 1.1·10⁻¹⁵ first m basic 1.99·10⁸ 3.80·10¹⁶ 0.50 5.2·10⁻² pivoted basic 2.06 31.9 0.50 6.7·10⁻¹⁶

The pivoted fundamental basis has the orthonormal one’s conditioning to within a factor of 1.25 and the naive one’s sparsity exactly. All of the sparsity and none of the loss, for one greedy pass over A.

What the basic set is, when nobody chose it #

The measurement above treats “which columns are basic” as a decision, and in most codes it is not made as one. It is inherited.

In a simplex-based method the basic set is the simplex basis, chosen by a pricing rule for reasons that are about the optimisation and not about conditioning at all. In a reduced-gradient method it is whichever variables are away from their bounds. In a finite-element code with multipoint constraints it is often the first variable each constraint mentions, because that is what the assembly loop reached first. None of those rules has any reason to produce a well conditioned A_B, and none of them reports κ(A_B).

The consequence is that the number this page is about — κ(Z), the quantity that decides how many digits the answer keeps — is set by a subroutine chosen for a different purpose. That is a recognisable shape on this site: the ordering that decides a factorisation’s memory is chosen by a graph heuristic that knows nothing about the values, and the tolerance that decides a rank is often a default nobody set.

What is different here is that a fix is cheap and local. Computing κ(A_B) costs an estimate; if it is large, one column-pivoted pass supplies a better basic set. The measurement’s value is that it says how large “large” has to be before it matters, and the answer is around 10⁸ for a five-per-cent error in double precision — which is √(1/u), and the next section is about why that is right for a different reason than the identity would give.

The square reaches the answer; the bound does not predict it #

One point is enough to observe that the square is attained in the condition number and not enough to say what it does to the answer. Swept over ten decades of γ, on the same construction:

The square reaches the answer. The fitted slope of the forward error against κ(Z) is 1.99 over the four clean decades, so the exponent that appears in the identity appears in the error too. That is the claim this page rests on and it holds exactly rather than approximately, which one point could not have established.

The bound does not. The last column runs from 1.7·10⁻² to 2.5·10⁻⁴, so κ(ZᵀHZ)·u overestimates the error by between 150 and four thousand times, and the overestimate itself varies by more than an order across the sweep. It is a bound and it is not a predictor, and the distinction matters exactly here, where somebody wants to know at what κ(A_B) to start worrying.

Which is why √(1/u) is right for a reason other than the one it looks like. Taking the identity at face value puts a five-per-cent error at κ(H)κ(Z)²u = 0.05 — that is, κ(Z) = 1.5·10⁷, thirteen times below the measured threshold of 2·10⁸. The √(1/u) heuristic, which involves neither H nor the identity, gives 6.7·10⁷ and lands within a factor of three. The rule of thumb beats the derivation, because the derivation is carrying a bound’s slack and the rule of thumb is carrying none.

There is a boundary at the far end too. Past γ = 10⁹ the error is 0.82 and then 0.32 — order one, and no longer monotone — so the sweep has run out of signal rather than continuing to worsen. The same shape a tensor that cannot be decomposed finds: an error that stops growing has not stopped mattering, it has finished, and past that point the conditioning is measuring how few starting configurations still work rather than how wrong the answer is.

The practical form of all this is a threshold a code can carry. The quantity available before anything is solved is κ(A_B), which a condition estimator supplies for the cost of a few solves with a factorisation the code already has, and κ(Z) is within a small factor of it on this family — 1.99 against 2 at every γ. So the rule is: estimate κ(A_B); if it is above about 10⁷, one column-pivoted pass over A is worth its O(m²n); if it is above 10⁹ the answer has no significant figures and the pass is not optional. Neither of those numbers comes out of the identity, and both come out of the sweep. What the identity supplies is the shape — that the threshold moves as the square root of the precision, so a code working in binary32 should read 10³ and 10⁴ instead, which is a range ordinary problems reach.

Every number in that paragraph is held to a test that fails if it moves: the slope of 1.99, the attainment of the square in κ(ZᵀHZ), the bound’s looseness and its variation, the position of the five-per-cent crossing against both √(1/u) and the identity, and the saturation past 10⁹.

What five per cent means here #

The naive basis’s forward error is 5.2·10⁻², measured against a solution computed in BigInt rationals from the same stored matrices. That is not a loss of precision in the usual sense — it is an answer wrong in its first significant figure.

And nothing about the computation announces it. The constraint is satisfied: Ax = g holds to the rounding level for every basis, because x_p was built to satisfy it and Z annihilates A exactly. The reduced system was solved by a Cholesky that completed without complaint, because ZᵀHZ really is positive definite. The residual of the reduced system is small. Every check a code would run passes.

The only quantity that says anything is κ(ZᵀHZ), and it is a number a code computes only if somebody asks it to. This is the same shape as the classical Gram–Schmidt failure, where the reconstruction ‖A − QR‖ stays at 10⁻¹⁶ while ‖QᵀQ − I‖ reaches one: the quantity that is wrong is not the quantity anybody looks at.

Why it is invisible to the constraint’s own conditioning #

κ(A) is 10⁶ at every point of the sweep. It does not move. The badness being measured is not in A; it is in a choice made about A, and no norm of A can see it.

That is worth separating from the second essay in the constraint field, where the range-space method’s error grew because κ(A) grew. There the problem got harder. Here the problem is fixed and the method’s internal decision is what moves.

The distinction is the site’s spine with a third author in it. Backward error is what the algorithm did; the condition number is what the problem did with it; and this is neither — it is what a reformulation did, and the reformulation was exact. Two essays here have now found the same author, and the other one is the elimination that squares.

What the theorem of the constraint preconditioner says about all this #

There is an apparent contradiction with the fourth essay in the constraint field, which states that the preconditioned spectrum is the generalised eigenvalues of the pencil (ZᵀHZ, ZᵀGZ) for any Z. If the basis matters as much as this page says, how can that be basis-free?

Because a generalised eigenvalue is invariant under congruence. Replacing Z by ZM sends both matrices of the pencil to Mᵀ(·)M, and det(MᵀXM − λMᵀYM) = det(M)²det(X − λY): the characteristic polynomial is scaled and its roots are not moved. Any two bases for one null space differ by such an M.

So the exact spectrum is basis-free and the computed one is not, and everything this page measures lives in the gap between those two sentences. That gap is the site’s whole subject, and it is unusually explicit here: an invariant that is invariant in the algebra and moves by eight orders in the arithmetic.

The sparsity that is being defended #

It is worth putting a number on what the fundamental basis is for, since the page has spent most of its length on its failure mode.

The orthonormal Z from a QR of Aᵀ has n(n − m) entries and essentially all of them are nonzero. The fundamental basis has an (n − m) × (n − m) identity plus an m × (n − m) block, so its nonzero count is (n − m)(m + 1) at worst and much less when A_B⁻¹A_N is itself sparse — which it is whenever the constraints are local, as they are for a mesh, a network or a set of multipoint ties.

On the figure’s small problem that is 0.50 against 1.00 in density, which is a factor of two. On a problem with n = 10⁵ and m = 10³ it is the difference between an object of 10⁸ entries and one of 10⁶, and the second is the only one anybody can hold. The reduced Hessian inherits the same distinction: ZᵀHZ formed from a dense Z is dense, and formed from a local Z keeps a banded structure that a sparse Cholesky can exploit.

So the naive rule is not a mistake to be replaced by the orthonormal basis. It is the right family of basis, chosen badly, and the repair is inside the family.

Two refusals, and what each guards #

The claim this page has to close is the one that reads as a definition rather than an assumption: any set of n − m independent directions annihilated by A is as good a basis as any other, since they all describe the same subspace and give the same answer. Every clause of that is true except the conclusion. The assertion is fed the naive basis’s reduced Hessian and required to reject the claim that it is within a factor of a million of the orthonormal one’s — and 3.8·10¹⁶ against 25.6 is a factor of 1.5·10¹⁵.

The second refusal guards the other direction. The fundamental basis has an identity block in it, and an identity block is orthonormal, so it is easy to conclude that the basis is. It is not: the other block is A_B⁻¹A_N and has no reason to be anything in particular. The assertion is fed κ(Z) = 1.99·10⁸ and required to refuse.

And the third covers the sparsity claim from the other side: a QR of Aᵀ is asserted to be dense, because the entire case for the fundamental basis rests on the orthonormal one not being sparse, and a reader is entitled to have that measured rather than assumed.

Two stops give two points and a direction, and a direction is not a rate — the difference between inheriting γ and inheriting its square is exactly what two points cannot distinguish. The slider can.

Two of the three columns have not moved at all across four decades of γ, which is the control the third column is read against: whatever is happening to the variable-reduction basis is not happening to the problem.

Two of the three bases are unaffected by γ entirely. κ(Z) is 1 for the orthogonal basis at all six stops and 2.06 for the third at five of the six, and both errors stay between 2.9·10⁻¹⁶ and 2.3·10⁻¹⁵ across ten decades of γ. Whatever the constraint matrix is doing, those two bases do not inherit it.

The variable-reduction basis inherits it exactly. κ(Z) reads 4.44, 206, 2.00·10⁴, 1.99·10⁶, 1.99·10⁸ and 1.99·10¹⁰ — dividing by γ gives 4.44, 2.06, 2.00, 1.99, 1.99, 1.99, so κ(Z) = 2γ to three figures from the second stop on. The basis nobody chose is exactly twice as badly conditioned as the block somebody happened to put first.

And the error is quadratic in that, not linear. Dividing the variable-reduction error by κ(Z)²·u gives 0.019, 0.012 and 0.012 at γ = 10⁴, 10⁶ and 10⁸ — constant to a quarter across four decades. So the cost is κ(Z)² times the unit roundoff, which is the squaring this collection keeps finding wherever a Gram matrix or a non-orthonormal basis appears, and it is why the difference between κ(Z) = 1 and κ(Z) = 2γ is not a factor of 2γ in the answer but a factor of 4γ².

The last row is at the ceiling rather than on the curve: an error of 0.317 is a third of the answer, and the run from 0.0518 to 0.317 is a factor of six where the previous steps were four orders. The quadratic law should be read from the middle of the table, and its own endpoint is where the answer has stopped having digits to lose.

A third basis this page does not build #

There is a middle option that appears in the literature and is worth naming so that the comparison does not read as exhaustive: an orthogonal fundamental basis, obtained by orthogonalising only the m × (n − m) block rather than the whole of Z.

It does not exist in general. The columns of the fundamental basis are not orthogonal to each other, and making them so mixes the identity block into itself — after which the identity block is gone and so is the sparsity. What can be done instead is to keep the fundamental form and precondition the reduced problem, which is exactly what the constraint preconditioner of the constraint field does and which leaves the conditioning of Z where it was while removing its effect on the iteration count.

That is a genuine third answer and it changes the question rather than answering it: the reduced Hessian is still κ(Z)² and is never formed, and what governs the work is the preconditioned spectrum instead. It is available only to an iterative method, and the two bases compared above are what a direct one has.