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), for the same reason it always is.

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.

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.

At other settings #