Where the multipliers go

A factorisation with nothing to pivot for found Cholesky needing no interchanges at all on a positive definite matrix, and when symmetry is not enough built the symmetric indefinite factorisation $PAP^{\mathsf T} = LDL^{\mathsf T}$, with $D$ block diagonal in one-by-one and two-by-two pieces, and measured Bunch and Kaufman’s rule for choosing between them. The rule bounds the growth of the entries in $D$. It does not bound the entries of $L$, and the essay drew the consequence in one sentence: “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.” It named the repair — the bounded variant, which searches along rows and columns the way rook pivoting does for unsymmetric elimination — and left the sentence as a warning.

The sentence makes three predictions, and only one of them turns out to be true.

A matrix that makes the multipliers large #

The classic way to make Bunch–Kaufman’s multipliers unbounded fits in three rows. Take $\begin{pmatrix} 0 & \varepsilon & 0 \\ \varepsilon & 0 & 1 \\ 0 & 1 & 1\end{pmatrix}$. The first diagonal entry is zero and the largest entry in its column is ε, so the rule looks at the second column, finds its largest entry is 1, sees the second diagonal is zero too, and takes the two-by-two pivot on the first two rows. That block is $\begin{pmatrix}0&\varepsilon\\ \varepsilon&0\end{pmatrix}$, its inverse has entries $1/\varepsilon$, and the third row’s multipliers are its entries times that inverse: $1/\varepsilon$.

Embedded in something larger, the pattern keeps its shape. The family here is 24 by 24: the three rows above, with the second and third also coupled to a seeded random symmetric block of 21 rows. The coupling ε runs from $10^{-1}$ to $10^{-10}$. Beside Bunch–Kaufman runs the bounded rule of Ashcraft, Grimes and Lewis: from the first column, walk to the row holding its largest off-diagonal entry, then to the largest entry in that row’s column, and so on, until either a diagonal entry is large enough to pivot on alone or an entry is found that is largest in both its row and its column, and take the two-by-two block around it. Every multiplier it produces is at most $1/(1-\alpha) = 2.78$ with Bunch–Kaufman’s constant α, and on the thirty random symmetric matrices used below as the control its worst is 1.56.

Bunch–Kaufman’s largest multiplier is 12 at ε = 0.1 and 1.2·10¹⁰ at 10⁻¹⁰, exactly $1.2/\varepsilon$ throughout; the bounded rule’s is 1.6 at every ε. The growth in D, the quantity Bunch–Kaufman was designed to control, is 2.51 for one rule and 2.57 for the other. So the family does what it was built to do: two backward-stable-looking factorisations of the same matrix, identical in growth, differing by ten orders of magnitude in the size of one factor.

What the solve depends on is a product #

The dashed lines in the figure are the quantity the backward error of a solve with these factors is actually written in. The standard analysis of a symmetric indefinite solve bounds the backward error by a modest multiple of the unit roundoff times $\||L||D||L^{\mathsf T}|\|$, not of ‖L‖ alone. For Bunch–Kaufman that product is 6.72 times ‖A‖ at every ε from $10^{-1}$ to $10^{-10}$; for the bounded rule it is 9.47. The multipliers grow by a factor of a billion and the product does not move, because in the columns where L is large, D is small in exactly the proportion that cancels it: the two-by-two block holds ε, and $1/\varepsilon$ times ε times $1/\varepsilon$ is paired in the product with the entries that make it come out at one.

The measurement agrees with the analysis rather than with the warning. Bunch–Kaufman’s backward error is 1.6·10⁻¹⁶ at every coupling, and after one step of refinement using the same factors it is 6.3·10⁻¹⁷ or 1.1·10⁻¹⁶. The bounded rule’s is slightly worse, between 4.8 and 6.1·10⁻¹⁶, and refinement brings it to the same place. Refinement solves with the factors and so, by the warning, should have inherited their size; it inherits the product, which is small, and does what refinement does.

The forward errors follow the matrix’s condition number, which grows as $1/\varepsilon^2$ — from 720 to 6.7·10²⁰ over the sweep — and both rules land within a small factor of it at every coupling, sometimes one ahead and sometimes the other. At $10^{-8}$ Bunch–Kaufman’s forward error is 4.9·10⁻⁷ and the bounded rule’s 2.6·10⁻⁸; at $10^{-10}$ they are 3.8·10⁻⁵ and 1.9·10⁻⁵. Neither is a property of the factors. Both are what a backward-stable solve of a matrix with that condition number is entitled to.

The same holds for a condition estimate of the kind libraries compute from the factorisation. It estimates $\|A^{-1}\|$ by solving with the factors a few times, and each of those solves is the backward-stable solve above. An estimate that can be fooled found the estimator’s failures in the vectors it happens to try, not in the factors it solves with, and nothing here gives it a new one.

Why the product stays small, in three rows #

The cancellation can be seen exactly on the three-row pattern, before the random block is added. Bunch–Kaufman’s factors there are a two-by-two block $E = \begin{pmatrix}0&\varepsilon\\ \varepsilon&0\end{pmatrix}$, a one-by-one pivot of 1, and a third row of $L$ equal to $(1/\varepsilon, 0)$ in the block’s columns. The product of absolute values has as its third row, in the block’s columns, $(1/\varepsilon, 0)\,|E| = (0, 1)$: the large multiplier meets the zero on E’s diagonal and the ε off it, and what comes out is one. Carried on through $|L|^{\mathsf T}$, the corner entry of $|L||D||L^{\mathsf T}|$ is $0 \cdot 1/\varepsilon + 1 \cdot 0 + 1 \cdot 1 = 1$, exactly the corresponding entry of $|A|$. Every entry of the product can be checked the same way, and none contains ε at all. The multiplier is $1/\varepsilon$ because the pivot is ε, and the analysis of the solve never sees one without the other.

That is the difference from unsymmetric elimination, where the swap that is not optional showed a large multiplier doing real damage. There, a multiplier of $1/\varepsilon$ multiplies a row of ordinary size and writes a large number into the Schur complement — growth — and the growth is what the backward error pays for. Here the large multipliers are paired with a pivot block whose inverse produced them, and the Schur complement they update is, in this family, left with ordinary entries. Partial pivoting bounds multipliers by one and lets growth run; Bunch–Kaufman bounds growth and lets multipliers run; and on this matrix the second is the one a solve can afford.

The conditioning has to go somewhere #

If the size of L is not an error, it is a measurement of something, and the condition numbers of the two factors say what.

At ε = $10^{-8}$ the matrix’s condition number is 6.7·10¹⁶. Bunch–Kaufman’s L has condition number 1.1·10¹⁷ and its D 7.3·10⁸. The bounded rule’s L has condition number 37.5 — at that coupling and at every other — and its D has 5.9·10¹⁶. The ill-conditioning of A is a fact about A and every exact factorisation has to hold it somewhere; the two rules put it in different factors. Bunch–Kaufman divides by the small ε in the pivot and so writes the large reciprocal into the multipliers; the bounded rule walks away from the ε coupling to an entry largest in its row and column, eliminates with bounded multipliers, and leaves the small quantity to surface in D as a small pivot.

That reading also dissolves the warning’s premise. A large multiplier is not a sign that something has gone wrong in the elimination, the way a large growth factor is for unsymmetric Gaussian elimination. It is a sign that the matrix is ill conditioned and that the rule chose to express that in L. The exact answer to a nearby problem is the reason that distinction matters: the backward error is about the algorithm, the condition number is about the problem, and a factor that stores the condition number is doing the second job, not failing the first.

The use that reads L by itself #

The third prediction is the one that holds, and it holds sharply. Some uses of the factorisation do not solve with all of it; they take one factor and interpret it. The most common in optimisation is extracting a direction of negative curvature: when the factored matrix is a Hessian that is not positive definite, the most negative eigenvalue of D, with its eigenvector carried back through $L^{\mathsf T}$, gives a direction along which the quadratic model decreases, and a trust-region or line-search method steps along it. A minimum the Hessian cannot see and a loop that asks the null space why both read the inertia of this factorisation; this reads its vectors.

The matrix’s smallest eigenvalue is −8.40 at every coupling — the random block decides it, and ε does not touch it. The bounded rule’s direction has a Rayleigh quotient of 13% of that at every coupling, the same number to six digits. Bunch–Kaufman’s has 5% at ε = 0.1, 0.07% at $10^{-2}$, and then falls by exactly $\varepsilon^2$: 6.9·10⁻¹⁶ of the curvature at $10^{-8}$, 6.9·10⁻²⁰ at $10^{-10}$. The direction is dominated by the $1/\varepsilon$ multipliers, which point it into the corner of the matrix where the ε coupling lives and where the curvature is of order $\varepsilon^2$. It is a direction of negative curvature in the exact sense — its quotient is negative — and a useless one.

This is the ‖L‖ the earlier essay meant. The direction is $L^{-\mathsf T}$ applied to one vector with nothing afterwards to undo it, so the multipliers go straight into the answer, and a factor holding $10^{17}$ worth of conditioning turns into a direction that is a thousand-million times too long in the wrong components. The same goes for a null-space basis built from the columns of L, for any preconditioner that uses L alone, and for a modified factorisation that changes D and keeps L — each of them applies one factor without its partner.

The multipliers live where the pivot divided by ε #

Column by column, the damage is local. Bunch–Kaufman’s large multipliers sit in the first two columns — the two-by-two pivot on the ε coupling — and every later column of its L is ordinary, below two. At $\varepsilon = 10^{-6}$ one column holds entries above ten and its largest is 1.2·10⁶; turn the dial and that column’s entries scale as $1/\varepsilon$ while the others stay put. The bounded rule’s columns are all below two at every setting, because its search never chooses the ε coupling as a pivot at all.

That locality is what makes the curvature failure so complete and the solve so unaffected. The solve uses the large column and the small pivot together, and they cancel. The curvature direction does not come from that pivot block: the ε block’s eigenvalues are $\pm\varepsilon$, near zero, and D’s most negative eigenvalue sits in a later two-by-two block of ordinary size, $-3.9$ or so. But the direction is carried back through $L^{\mathsf T}$ by back substitution, and back substitution adds to each component the multipliers below it times the components already computed. The first two components receive the $1/\varepsilon$ multipliers of the first two columns, so a direction that started in the middle of the matrix ends up a million times longer in its first two coordinates than anywhere else, pointing into the corner where the ε coupling lives and the curvature is of order $\varepsilon^2$. The large column and the small pivot cancel in a solve because the solve uses both; the curvature direction uses the column and never meets the pivot.

The control #

None of this happens by itself. On thirty random symmetric matrices with standard normal entries, Bunch–Kaufman’s largest multiplier is between 1.26 and 3.91, and the bounded rule’s between 1.19 and 1.56, inside its bound of 2.78. On the regularised saddle-point matrices of the earlier essay the bounded rule’s multipliers are exactly one. Large multipliers need a structure that asks for them — a small coupling between a row with no diagonal and the rest — and saddle-point systems with a weakly coupled constraint are exactly where optimisation codes meet one. A pivot that searches one row and one column made the same observation about rook pivoting for unsymmetric matrices: the search costs almost nothing on typical matrices and earns its keep on the structured ones.

What the warning should have said #

The earlier sentence is right about one of its three consumers and wrong about two, and the dividing line is clean. A use that applies the whole factorisation — a solve, a refinement, a condition estimate built from solves — inherits $|L||D||L^{\mathsf T}|$, and Bunch–Kaufman keeps that bounded. A use that applies one factor alone — a curvature direction, a basis read from L, a factorisation whose D is modified — inherits ‖L‖, and Bunch–Kaufman does not keep that bounded. Which rule to use is therefore not a question about accuracy in general. It is a question about what the caller does with the factors after the solve, which a library cannot know and the caller always does.

None of this touches what the factorisation is most often used for besides solving. The inertia — how many positive, negative and zero eigenvalues the matrix has — is read from the signs of D’s blocks, and both rules report the same counts at every coupling, because Sylvester’s law of inertia holds for any congruence however badly conditioned its factor. The shift that stops at the first right count relies on exactly that, and it is safe with either rule; it is the vectors, not the counts, that the multipliers reach.

The cost of choosing the bounded rule is small where it has been measured. It takes the same number of two-by-two blocks on the random matrices within one or two, its growth is 2.57 against 2.51 here, and the extra comparisons of its search are, in the unsymmetric measurements, a few per cent of the elimination. What it buys is a factorisation in which L is well conditioned by construction, so that the conditioning of the problem appears where a reader of D can see it — as a pivot of size ε — rather than spread through a factor that one consumer in three reads alone.

What one planted structure does not show #

One family, at one size, with one random block drawn and two more checked for the same behaviour; the ε coupling is placed in the first row, where Bunch–Kaufman meets it first. The same coupling deeper in the matrix is met after other pivots have changed it, and the multipliers then depend on what those pivots did. The curvature direction measured is the simplest one — the most negative eigenvalue of D, carried back once — and the methods that use this factorisation refine it; whether a refinement step recovers what the large multipliers lose is not measured. And the condition estimate is argued rather than run.

Still open: a curvature direction that is refined, and the ε coupling met late #

Refining the direction. A Newton-type method rarely steps along the raw direction; it takes a few iterations of a Lanczos or inverse iteration from it. The prediction with a sign is that inverse iteration with the same factors recovers the curvature within two or three steps, because each step is a solve and the solve is stable — which would make the bounded rule’s advantage a matter of the first step only.

The coupling met late. Put the ε coupling in the last rows instead of the first, where Bunch–Kaufman reaches it after twenty pivots that have mixed it into the Schur complement. Whether the multipliers still grow as $1/\varepsilon$ or whether the earlier eliminations have diluted the coupling into something the rule handles with bounded multipliers is the measurement that would say how often a real matrix triggers this.

A modified factorisation. Optimisation codes often replace D by a positive definite $\tilde D$ and use $L\tilde DL^{\mathsf T}$ as a model Hessian. The difference from A is $L(\tilde D - D)L^{\mathsf T}$, and with Bunch–Kaufman’s L it could be ten orders of magnitude larger than the change made to D. That is the second use of one factor alone, and the prediction is that it fails the same way.