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 — including everything a constraint written as a weight at infinity covers.

The 2×2 pivot #

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

E = [ aₖₖ  aₖᵣ ]
    [ aᵣₖ  aᵣᵣ ]

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ₖₖ| ≥ α·λ. 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.

The bound for a run of 1×1 pivots is (1 + 1/α)ⁿ⁻¹, and α = (1 + √17)/8 ≈ 0.6404 is the constant every implementation uses. Saying that α minimises that bound is the natural next sentence and it cannot be right: (1 + 1/α)ⁿ⁻¹ falls monotonically in α, so it is minimised by taking α as large as the rule allows and has no stationary point anywhere.

What α balances is the two paths against each other. Two 1×1 steps grow by at most (1 + 1/α)²; one 2×2 step, covering the same two rows, by at most (3 − α)/(1 − α) — which increases in α where the first decreases. Setting them equal:

(1 + 1/α)²  =  (3 − α)/(1 − α)      at  α = (1 + √17)/8,

both sides 6.5616, agreeing to 9·10⁻¹⁶. So the constant minimises the larger of two bounds rather than one bound, which is why it sits in the interior and why it is an algebraic number rather than a round one.

And the choice barely matters #

Sweeping α over two hundred and fifty random 8×8 symmetric indefinite matrices plus three saddle points:

   α        mean growth    worst |L|    2×2 blocks taken
  0.2          3.689         21.27            325
  0.5          1.737          6.32            517
  0.6404       1.556          6.32            610
  0.8          1.502          6.96            694
  0.9          1.516         19.85            720

The measured growth is flat from α = 0.5 upward — 1.50 to 1.74, a sixth — while the quoted bound moves over that range by a factor of twelve. The classical constant is not the measured optimum; 0.8 is, by 3.6%, which is far inside the scatter. This is the same relationship the essay’s other numbers have: the bound is exponential, the measurement is between 1.00 and 1.33 on the saddle-point family, and the constant chosen to optimise the bound is optimising something that does not happen.

The last column is what α is actually for — it is monotone in how often a 2×2 block is taken, 325 to 720 across the sweep — and the one before it is the price of pushing it too far. The worst entry of L is 6.32 up to α = 0.8 and 19.85 at 0.9, which is the weakness Bunch–Kaufman is known for and which the rules below are built to fix, located on the parameter that controls it.

That column is worth one more look, because it makes the case for the alternative rules concretely rather than by citation. At α = 0.9 the growth in D is 1.52 — better than the classical constant’s 1.56 — and the largest entry of L is three times larger. A code choosing α by the number the essay’s figures report would choose 0.9 and would be choosing the worse rule, because the quantity it improved is not the quantity that binds.

Growth in D is what the bound is about and the size of L is what the factors are used for. A condition estimate reads L, an iterative refinement solves with it, a null-space basis is built from it, and every one of those inherits ‖L‖ rather than the growth factor. So the essay’s own headline comparison — Bunch–Kaufman’s 1.25 against the diagonal rule’s 5·10⁷ — is a comparison on the axis where the classical rule wins outright, and there is a second axis on which it is merely adequate and rook pivoting is the repair.

That is the honest reason there are four rules rather than one, and it is sharper than “no single choice dominates”: the rules optimise different quantities, the quantity a bound is stated in is not always the quantity a caller depends on, and α is a single knob that trades between them monotonically.

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⁷, 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.

The constant is a half, and it is worth writing down rather than hiding in a c. Growth 1/(2τ) predicts 50, 5,000, 5·10⁵ and 5·10⁷ against a measured 50, 5,019, 5·10⁵ and 5·10⁷ — the worst disagreement is four parts in a thousand, at the stop where τ is largest and the asymptotic form has least room. A reader who carries away 1/τ will over-estimate the damage by a factor of two at every τ, which is a small error and the wrong kind: it is in the constant of the thing the essay is about.

Two more decades of regularisation, which is the direction a penalty method moves in as it converges:

Bunch–Kaufman’s readout is identical at τ = 0, 10⁻⁸, 10⁻⁶, 10⁻⁴ and 10⁻² — four blocks, growth 1.254, residual 5.8·10⁻¹⁷, to every digit printed. It is not that the block rule degrades more slowly than the diagonal one; it does not degrade at all, because the quantity that is going to zero is never in a denominator.

And the two rules’ residuals are the same statement about growth read twice. Divide each residual by its own growth factor and the slider returns 4.4·10⁻¹⁷, 6.2·10⁻¹⁷, 2.6·10⁻¹⁷, 3.2·10⁻¹⁷ and 4.5·10⁻¹⁷ — constant within a factor of 2.4 across seven orders of magnitude of growth, and about a fifth of machine epsilon. The backward-error argument is not being illustrated here, it is being measured: the residual is what the growth makes it, and the pivot rule’s whole job is the growth.

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 σₘᵢₙ(B)²/‖H‖, which come from the constraint block being coupled through H⁻¹. So

κ  ≈  ‖H‖ · max(1, ‖H‖/σₘᵢₙ(B)²)

and a nearly redundant constraint — two rows of B nearly parallel — takes σₘᵢₙ(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/(2τ) 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.