When symmetry is not enough

A factorisation with nothing to pivot for ended on the one condition Cholesky needs: the matrix must be positive definite. This essay is about what happens on a symmetric matrix that is not, and the answer is not a small adjustment.

Consider

[ 0  1 ]
[ 1  0 ]

It is symmetric. Its eigenvalues are +1 and −1, so it is nonsingular and its condition number is exactly 1 — as well conditioned as a matrix can be. And it has no nonzero diagonal entry.

A factorisation PAPᵀ = LDLᵀ with D diagonal has to divide by a diagonal entry of the permuted matrix at its first step. A symmetric permutation moves diagonal entries to other diagonal positions; it cannot move an off-diagonal entry onto the diagonal, because that is what makes it symmetric. So every diagonal entry available to be chosen is zero, and the factorisation does not exist.

Not “is inaccurate”. Does not exist. No amount of searching helps, and no amount of precision helps, and the matrix is perfectly conditioned.

Why a symmetric permutation is the constraint #

The obvious response is to permute rows only, as LU does, and take the 1 into the pivot position. That works and it is not a symmetric factorisation any more: PA is not symmetric, L and U are unrelated, and the storage and arithmetic halved by symmetry are back.

Which is the whole trade. Symmetry buys half of everything and costs the freedom to permute rows without permuting columns. For a positive definite matrix that costs nothing, because there was never a reason to permute. For an indefinite one it costs the existence of the factorisation.

The saddle point, where the zero is not an accident #

The 2×2 above is easy to dismiss as a puzzle. The structure it has is the structure of one of the most common systems in computational science.

[ H   Bᵀ ] [x]   [f]
[ B   0  ] [λ] = [g]

The bottom-right block is zero by construction. It is the block of Lagrange multipliers in a constrained optimisation, the pressure block in an incompressible flow discretisation, the dual block in an interior-point step, and the multiplier block in a constrained least-squares problem. Nobody put a zero there by accident; the zero is what says the constraint has no term of its own.

So a factorisation that needs a nonzero diagonal entry there has no matrix to work on, and the problem class in which this arises is enormous.

The 2×2 pivot #

Bunch and Kaufman’s answer is to allow a pivot block of size two: take

E = [ a_kk  a_kr ]
    [ a_rk  a_rr ]

whole, invert it — it is 2×2, so the inverse is four divisions — and eliminate two variables at once. Symmetry survives, because the block is symmetric and the elimination is A ← A − LELᵀ. The factorisation exists for every symmetric nonsingular matrix, and D is block diagonal with 1×1 and 2×2 blocks rather than diagonal.

On [[0, 1], [1, 0]] the whole matrix is one 2×2 pivot, D is the matrix itself, L is the identity, and the residual is exactly zero — asserted as an equality rather than a tolerance, because there is no arithmetic in the factorisation of a block that is already its own D.

When to take one, which is where the strange constant comes from #

The rule has to decide, at each step, between a 1×1 pivot and a 2×2. Taking 2×2 always would waste the definite case; taking 1×1 whenever the diagonal is nonzero would divide by numbers that are nearly zero.

Bunch and Kaufman’s rule compares the largest off-diagonal magnitude in the current column, λ, against the diagonal entry, and takes a 1×1 pivot when |a_kk| ≥ α·λ. Below that it looks one step further — at the largest entry in the row of the offending candidate — and either takes a different 1×1 pivot or takes the 2×2 block.

α = (1 + √17)/8 ≈ 0.6404 is the value minimising the resulting growth bound, and the bound is (1 + 1/α)ⁿ⁻¹. That is exponential and looks alarming beside Cholesky’s exactly-1, and in practice the measured growth is small: the assertions here find 1.00 to 1.33 across every matrix and every τ drawn, against a bound in the thousands.

What the obvious repair does instead #

The clean failure at τ = 0 is easy to dismiss as a boundary case. What makes the block rule necessary rather than tidy is what happens near it.

Regularise the primal block by τI — which is what a penalty method, an augmented Lagrangian, or a proximal term does — and the diagonal is no longer zero. Both 1×1 rules now succeed. They succeed by dividing by τ.

The measurement: at τ = 10⁻², 10⁻⁴, 10⁻⁶ and 10⁻⁸ the diagonal rule’s growth factor is 50, 5.0·10³, 5.0·10⁵ and 5.0·10⁷ — exactly c/τ — and its residual follows, reaching 2.2·10⁻⁹ at the smallest τ. Bunch–Kaufman’s growth is 1.25 at all four and its residual is 5.8·10⁻¹⁷ at all four.

Seven orders of magnitude of accuracy, on a matrix where the 1×1 rule did not fail and reported nothing. That is the case that ships.

Where the blocks are taken, and where they are not #

The other half of the claim is that the generality costs nothing when it is not wanted, and it is the half this site has been wrong about before, so it is measured.

On a positive definite matrix, Bunch–Kaufman takes zero 2×2 pivots and makes zero interchanges. The rule reduces to no pivoting at all, Cholesky’s argument survives intact, and the residual is at rounding. On the saddle-point matrix at τ = 0 it takes n/2 blocks — one for every pair — and at τ = 1, where the diagonal is as large as the coupling, it takes almost none.

The other rules, and why there are several #

Bunch–Kaufman is what LAPACK’s xSYTRF runs by default, and it is not the only rule, for a reason worth stating.

Bunch–Parlett searches the whole remaining submatrix, like complete pivoting, and has a polynomial growth bound rather than an exponential one. It costs O(n³) comparisons, which is the same objection complete pivoting faces in the unsymmetric case.

Bounded Bunch–Kaufman, or rook pivoting, alternates between a column search and a row search until an entry is largest in both. It bounds the entries of L, which Bunch–Kaufman does not: the classical rule can produce an L with arbitrarily large entries even while the growth in D is modest, and a large L is a problem for anything that later uses the factors — a condition estimate, an iterative refinement, or a null-space basis.

Aasen’s method reduces to tridiagonal LTLᵀ instead of block diagonal, which bounds L by construction and moves the difficulty into a tridiagonal solve.

That there are four rules rather than one is the signature of a problem where no single choice dominates, and it is the same signature the unsymmetric case has: partial, scaled partial, rook and complete, each buying a different thing.

The alternatives that avoid the problem, and what they cost #

A saddle-point system can be solved without a symmetric indefinite factorisation at all, and the two standard ways are worth naming because both trade the difficulty for a different one.

The Schur complement, or range-space method. Eliminate x first: from the top block x = H⁻¹(f − Bᵀλ), substitute, and solve (B H⁻¹ Bᵀ) λ = B H⁻¹ f − g. The Schur complement S = B H⁻¹ Bᵀ is symmetric positive definite when H is, so Cholesky applies and the whole of the previous essay comes back.

What it costs: H⁻¹ appears explicitly, so H must be cheaply invertible — diagonal, or block diagonal, or a mass matrix — and if it is not, forming S is more expensive than the original system. And S is dense even when B is sparse, because B H⁻¹ Bᵀ couples every pair of constraints that share a variable. The conditioning is worse too: κ(S) is roughly κ(H)·κ(B)², so a route chosen to get back to Cholesky can square a condition number to do it — which is the road that squares the problem reappearing in a different field.

The null-space method. Find Z with BZ = 0, solve the reduced problem (ZᵀHZ) y = Zᵀ(f − H x₀) in the null space of the constraints, which is again SPD when H is definite on that subspace.

What it costs: a null-space basis. Computing one that is well conditioned needs a QR of Bᵀ; computing one cheaply — by partitioning B into a square block and the rest — gives a Z whose conditioning is the conditioning of that block, which nothing chose on numerical grounds. And if H is only definite on the null space, which is the usual assumption in optimisation, the reduced problem is fine and the full one is indefinite — so the method needs a Z before it can say anything at all.

Against both of those, a symmetric indefinite factorisation of the whole system needs no assumption on H beyond nonsingularity of the pair, no explicit inverse, and no basis. That is why it is the default in a general-purpose code, and why the pivot rule inside it had to be invented.

The conditioning, which the factorisation does not fix #

A saddle-point matrix’s condition number does not behave like H’s or B’s, and the block structure is what makes it legible.

The eigenvalues separate into two groups: those of order ‖H‖, which come from the primal block, and those of order σ_min(B)²/‖H‖, which come from the constraint block being coupled through H⁻¹. So

κ  ≈  ‖H‖ · max(1, ‖H‖/σ_min(B)²)

and a nearly redundant constraint — two rows of B nearly parallel — takes σ_min(B) towards zero and κ up quadratically. That is a property of the problem and not of the factorisation, and no pivot rule in this essay improves it.

What the factorisation controls is the growth, which is the part the algorithm is responsible for, and the separation is exactly the site’s own identity: forward error ⪅ κ × backward error, with κ belonging to the constraints and the backward error belonging to the pivot rule. The measurement above — 1.25 against 5·10⁷ — is entirely about the second factor.

Why not just use LU #

Because half of everything is a lot, and because the alternative loses something structural as well as something quantitative.

An unsymmetric LU of a symmetric matrix computes both triangles, stores both, and performs twice the arithmetic. On a sparse saddle-point system it also computes a different elimination tree, because the row permutation and the column permutation are no longer tied, and the fill it produces has no reason to respect the symmetry of the sparsity pattern.

And it discards the inertia — the counts of positive, negative and zero eigenvalues — which LDLᵀ hands over for free. Sylvester’s law of inertia says the signs in D are the signs of the eigenvalues, so a symmetric indefinite factorisation answers “how many negative eigenvalues does this Hessian have” as a by-product, at no cost, and that question is the one an optimisation code is actually asking. A 2×2 block contributes one positive and one negative, since its determinant is negative — which is why the rule takes them.

What the residual is measured against #

A detail about the badge, because a permuted factorisation makes the obvious choice wrong.

The quantity printed is ‖PAPᵀ − LDLᵀ‖_F / ‖A‖_F, with the permutation applied to A rather than undone from the factors. Both are defensible and they differ in what they blame: the first asks whether the factorisation of the permuted matrix is accurate, the second whether the factors reconstruct the original.

They are the same number here, because a permutation is orthogonal and the Frobenius norm is orthogonally invariant. That is worth checking rather than assuming — it fails for a norm that is not unitarily invariant, and the ∞-norm this site uses elsewhere is one of those.

Two rules that are the same rule #

ldl’s three pivot options are written through one body, and the reason is the site’s standing practice rather than economy: what is being compared has to be the rule and not two implementations of an elimination.

The none and diagonal branches differ by a single symmetric interchange chosen before the step; everything after that — the divisions, the rank-one update, the ordering of the arithmetic — is character-for-character the same code. So when the diagonal rule’s growth is 5·10⁷ and the block rule’s is 1.25, the difference is the pivot and cannot be anything else.

That is the same discipline pivotStrategies applies to the unsymmetric case in the pivot that reads the units, and it is why both essays can attribute a factor of 10⁷ to a rule rather than to a routine.

What is worth carrying #

A symmetric factorisation with 1×1 pivots does not always exist, and the matrix that shows it is 2×2, perfectly conditioned and has a zero diagonal by construction rather than by accident.

The failure at zero is the easy case. The one that ships is small-but-nonzero, where the 1×1 rule succeeds, divides by τ, and returns a factorisation whose growth is 1/τ with nothing reporting it.

A 2×2 pivot costs nothing where it is not wanted. Zero blocks and zero interchanges on a positive definite matrix, measured, so the general routine is not a tax on the common case.

And the block structure is information. The inertia falls out of the signs in D, which is the question a constrained optimisation is asking anyway, and no unsymmetric factorisation offers it.