Pivoting and sparsity

9/17/2026; edited 9/27/2026

The first note ended with two different reasons to permute a matrix. A fill-reducing ordering chooses an elimination order from the nonzero pattern before the numerical factorization begins. Numerical pivoting chooses acceptable pivots from the values that appear during elimination. The first tries to make the factors cheaper. The second tries to keep the elimination numerically stable.

For a matrix known to be symmetric positive definite, this conflict disappears in the standard Cholesky setting. Cholesky needs no stability pivoting, so the permutation can target fill alone. The conflict returns for a general symmetric indefinite matrix. A symbolic analysis can predict the work for its chosen order, but it cannot determine whether the corresponding pivots will be safe. This is the central difficulty in a sparse indefinite LDL factorization.

We consider nonsingular symmetric systems; singular systems require additional rank-deficiency handling. The target factorization has the form

$$PAP^T=LDL^T,$$

where $P$ applies one permutation to both rows and columns and $D$ contains $1\times1$ and $2\times2$ diagonal blocks. The ordering chooses a prospective elimination order. Numerical pivoting chooses the pivot blocks and may revise that order.

What symbolic factorization promises

Let $A=A^T\in\mathbb{R}^{n\times n}$, and associate an undirected graph with its structural nonzero pattern. Vertex $i$ represents row and column $i$. For $i\neq j$, an edge $(i,j)$ is present when $a_{ij}$ is structurally nonzero. Diagonal entries do not create graph edges, and the numerical values of the off-diagonal entries do not matter at this stage.

Let $P_s$ be a permutation matrix, form $\widetilde A=P_sAP_s^T$, and perform scalar elimination on $\widetilde A$ in the order $1,\ldots,n$, with every scalar pivot nonzero. Eliminating a vertex joins all of its remaining neighbors into a clique. An edge added this way is fill. Repeating the operation produces a filled graph, and the later neighbors of vertex $k$ give the structural nonzeros below the diagonal in column $k$ of $L$.

This graph calculation is the symbolic factorization. It predicts the factor pattern, storage, and much of the arithmetic without computing a Schur complement. These predictions are structural. An exact numerical cancellation may make a predicted entry zero, but an implementation cannot rely on that accident when allocating storage.

The symbolic phase can also compare candidate orderings. A good ordering tries to keep the neighbor sets small, which reduces both fill and work. Finding the ordering with minimum fill is hard in general, so practical solvers use heuristics such as minimum degree and nested dissection. Davis (2006) gives a careful account of this separation between ordering, symbolic analysis, and numerical factorization.

The promise is especially useful for a sequence of matrices with the same pattern. If the ordering and elimination decisions remain fixed, the ordering and symbolic factorization can be computed once and reused. Only the numerical factors need to be recomputed when the values change. Multiple right-hand sides require even less work: the same numerical factors can be reused without refactorization.

There is an important condition hidden in this promise. The predicted pattern belongs to a fixed elimination order and fixed pivot blocks. If numerical stability changes either one, the symbolic plan is no longer guaranteed to describe the factorization exactly.

The same pattern can demand a different pivot

A four-vertex path gives a small example. Consider

$$A(d)= \begin{bmatrix} d&1&0&0\\ 1&2&1&0\\ 0&1&2&1\\ 0&0&1&2 \end{bmatrix}.$$

For every nonzero $d$, the structural graph is the path $1-2-3-4$. Symbolic elimination in the order $1,2,3,4$ predicts no fill because each eliminated vertex has at most one remaining neighbor. The symbolic prediction is therefore the same for $d=1$ and $d=10^{-2}$, even though $A(1)$ is positive definite and $A(10^{-2})$ is indefinite.

The numerical decisions are different. At a given elimination step, let $a_{ik}$ denote an entry of the current Schur complement, not necessarily the corresponding input entry. A common threshold test accepts a nonzero scalar pivot $a_{kk}$ only if

$$|a_{kk}| \geq \tau\max_{i\neq k}|a_{ik}|,$$

where $0<\tau<1$ is chosen by the solver and the maximum runs over the remaining indices. If the test passes, the multipliers in that pivot column have magnitude at most $1/\tau$. With $\tau=10^{-1}$, the planned first pivot is accepted for $A(1)$ and rejected for $A(10^{-2})$.

A block pivot needs its own test. If the current matrix is partitioned as

$$\begin{bmatrix}B&C^T\\C&S\end{bmatrix},$$

then eliminating the $1\times1$ or $2\times2$ block $B$ produces the block multipliers $CB^{-1}$. Practical block tests bound these multipliers or the associated update growth. Merely checking that $B$ is nonsingular is not enough.

Suppose the solver responds by exchanging vertices 1 and 2 and eliminating the second vertex of the path first. Its neighbors are 1 and 3, so the graph acquires the fill edge $(1,3)$. The exact Schur complement on vertices $1,3,4$ is

$$\begin{bmatrix} -49/100&-1/2&0\\ -1/2&3/2&1\\ 0&1&2 \end{bmatrix}.$$

The new entry $-1/2$ is precisely the fill predicted by eliminating vertex 2 first. Nothing about the input pattern changed. The numerical value of one diagonal entry changed the pivot order, and the new order changed the factor pattern.

Exchanging vertices is not the only response. The leading $2\times2$ block can instead serve as a block pivot. For $d=10^{-2}$, the largest entry of $CB^{-1}$ has magnitude $50/49$, well below $1/\tau=10$, so the block passes the analogous entrywise multiplier test. This choice creates no new off-diagonal entry in the trailing Schur complement. It nevertheless creates factor fill within the two-column pivot block: multiplication by $B^{-1}$ mixes the pivot columns, giving $L_{31}=50/49$ even though $A_{31}=0$. It also replaces two planned scalar pivots with one block pivot.

The example says something more basic: the nonzero pattern alone cannot determine which numerical path a sparse indefinite factorization will take. A dynamically pivoted numerical phase must be allowed to depart from the original scalar order. If a planned pivot fails the acceptance test, preserving both the scalar order and the scalar pivot structure requires weakening the test or perturbing the matrix.

Pivoting inside a sparse factorization

Searching the whole active matrix after every rejected pivot would discard much of the benefit of symbolic analysis. Sparse solvers instead combine a global structural ordering with local numerical freedom. The multifrontal method makes this compromise clear.

The symbolic ordering organizes the factorization into fronts. A front is a dense work array containing some variables ready for elimination together with variables that couple them to later work. Schur-complement updates from earlier fronts are assembled there. A variable is fully summed once every contribution to its row and column has arrived. Only then can its value be tested reliably as a pivot.

Within the fully summed part of a front, an indefinite solver searches for acceptable $1\times1$ and $2\times2$ pivots. The candidate search is restricted to the fully summed variables. In the classical threshold test, the multiplier test accounts for their coupling to all remaining variables represented in the front, including those in the update part. This is much cheaper and more predictable than complete pivoting over the whole reduced matrix. It also means that an acceptable pivot elsewhere in the matrix may not yet be available. Some restricted and parallel variants inspect only a subset of the update variables. They reduce search or communication at the cost of weakening the multiplier safeguard.

If a fully summed variable has no acceptable pivot in the current front, the solver can delay it. The variable and its assembled entries are passed to the parent front, where more candidates are available. This preserves the numerical safeguard, but the parent is now larger than symbolic analysis predicted. Several delayed pivots can meet in one ancestor, enlarge a dense frontal matrix, and increase work and storage far more than their count alone suggests. Duff and Reid (1983) introduced this strategy in the original multifrontal method.

Delaying a pivot postpones the decision rather than guaranteeing acceptance. At the root front, which has no parent, the solver must resolve any remaining unacceptable pivots or report failure.

The threshold $\tau$ exposes the local tradeoff. For a fixed Schur complement, increasing $\tau$ tightens the acceptance test and the multiplier bound. Over a complete factorization, this can cause more searches and delays. Threshold pivoting therefore supplies a tunable numerical test, not a simultaneous global bound on fill and stability.

The symmetric rook strategies of Ashcraft, Grimes, and Lewis (1998) illustrate why stronger searches are difficult to import unchanged. A front-restricted implementation must update and inspect additional rows and columns while choosing $1\times1$ or $2\times2$ diagonal blocks. If a search points outside the fully summed block, the solver cannot simply pivot there; it must continue searching, weaken the test, or delay a variable. These searches and data movements can add overhead and interfere with large blocked updates.

Preprocessing moves the conflict

Scaling and matching try to make the symbolic prediction numerically plausible before elimination starts. Diagonal scaling reduces disparities among rows and columns. Weighted matching seeks large entries that can serve as prospective pivots. For a symmetric problem, the preprocessing must preserve symmetry, and a large off-diagonal match naturally suggests a $2\times2$ block rather than an unsymmetric row permutation.

The prospective $1\times1$ and $2\times2$ pivots can then constrain the fill-reducing ordering. This mixes numerical information into the symbolic phase instead of waiting for a failure during factorization. Duff and Pralet (2005) develop symmetric scaling, matching, and constrained orderings for this purpose.

Preprocessing does not remove value dependence. Schur complements contain updated values that were not present in the input matrix, so a promising original entry can still become an unacceptable pivot later. Scaling and matching can reduce the number of surprises, but they do not guarantee that every pivot will be accepted. Scaling can improve the conditioning of the transformed system, but it does not guarantee good conditioning or eliminate sensitivity inherent in the original problem.

Dynamic and static responses

A dynamic strategy keeps the numerical threshold and changes the factorization when the symbolic choice fails. It may exchange candidates within a front, use a $2\times2$ block, or delay a variable. The computed factors then belong to the original matrix, up to rounding error, but their structure and cost may exceed the symbolic estimate.

Static pivoting makes the opposite choice. It follows the planned structure and perturbs a small or unacceptable pivot so elimination can continue. This makes the memory schedule predictable and communication easier, which is attractive in parallel solvers. The factors now represent a perturbed matrix by design. Iterative refinement may recover an accurate solution to the original system, but convergence is not automatic, and a perturbation can corrupt inertia information needed by an optimization algorithm. Mixed strategies first accept pivots under numerical tests, then weaken the test or perturb only after too many delays. Duff and Pralet (2007) analyze such combinations.

The main devices can be summarized as follows.

No row of the table dominates the others. A robust sparse solver is a collection of compromises among them. MA57, for example, separates analysis from numerical factorization, supports symmetric scaling and frontal $1\times1$ or $2\times2$ pivoting, can restart if its storage estimate is exceeded, and provides refinement and error analysis.

What can actually be reused

It is easy to say that the symbolic factorization is reusable when the pattern stays fixed. For an indefinite problem, it is better to say exactly which object is being reused.

• The ordering can usually be reused. It remains a useful structural starting point even if later pivoting departs from it.

• The assembly tree and symbolic fronts can often be reused as a scaffold. Delayed pivots may enlarge the fronts or require extra workspace during each numerical factorization.

• A preallocated structural envelope can be reused if the numerical factorization stays within it. Reusing the exact factor pattern and storage offsets additionally requires compatible pivot blocks, permutations, and storage conventions. A code can reserve extra space or restart when the envelope is exceeded.

• The numerical factors are reusable for new right-hand sides, but not for a matrix whose numerical entries have changed.

This distinction matters in parameter studies and optimization. A sequence of KKT matrices may keep the same block pattern while its diagonal values and inertia change at every iteration. Reusing the ordering can still save substantial work. Reusing an exact pivot sequence without tests is a stronger decision. It means accepting unchecked pivots, perturbing them to enforce a static plan, or assuming numerical behavior that the pattern alone does not guarantee.

For SPD matrices, the hierarchy collapses to the simpler case. A fixed pattern and permutation determine the symbolic Cholesky structure. Numerical cancellation may make some predicted entries zero without invalidating the storage plan. The symbolic analysis, storage plan, and elimination dependencies can all be reused across new positive definite values. This is one reason the symbolic and numerical phases separate so cleanly in sparse Cholesky.

Diagnostics belong to the result

A successful return code is not enough to evaluate the compromise. Useful structural diagnostics are the number of delayed pivots, the number and magnitudes of perturbations, the nonzeros in the factors, the largest front, peak memory, and any restart caused by underestimated storage. Comparing these quantities with the symbolic estimates shows how far numerical pivoting moved the computation from its plan.

Numerical diagnostics answer a different question. Solvers should report scaling information, estimates of element growth or conditioning when available, and the behavior of iterative refinement. For a computed solution $\hat{x}$ with residual $r=b-A\hat{x}$, the componentwise relative backward error

$$\eta = \max_i \frac{|r_i|}{\left(|A||\hat{x}|+|b|\right)_i}$$

is the smallest common bound for which perturbations satisfying

$$|\Delta A|\leq\eta|A|,\qquad |\Delta b|\leq\eta|b|, \qquad (A+\Delta A)\hat{x}=b+\Delta b$$

can explain the residual, with the usual conventions for a zero denominator. Here $\Delta A$ is not required to be symmetric; this is an unstructured componentwise backward error. The bound preserves structural zeros in $A$ and $b$. It is more informative than a normwise residual when rows and variables have very different scales. Arioli, Demmel, and Duff (1989) develop this measure for sparse systems and connect it to iterative refinement.

A small backward error says that the computed solution solves a nearby problem. It does not say that the forward error is small when the original problem is ill-conditioned. Pivoting tries to control errors introduced by the factorization; it cannot remove sensitivity already present in $A$.

The structural machinery comes next

The symbolic phase sees a graph while the numerical phase sees Schur-complement values. Sparse indefinite solvers reconcile the two with threshold tests, block pivots, delayed variables, preprocessing, perturbations, and refinement. Each device protects one part of the computation by giving up certainty somewhere else.

This note used fronts only to explain the numerical consequences of pivoting. The next note will develop the structural machinery itself in the cleaner SPD setting: elimination trees, approximate minimum degree, supernodes, and multifrontal factorization. It will then return to indefinite matrices and QDLDL to show how numerical pivoting and quasi-definite structure interact with the same machinery.