Deferring Numerical Pivoting with the Woodbury Matrix Identity

I propose here a way to factorize a sparse symmetric and indefinite matrix \(A\) with what I’ll call “Woodbury deferred pivots”, which allow us to retain a static execution plan even in the presence of numerical pivoting. A lot of the sparse linear solver research I do focuses on methods where we can obtain the execution graph statically purely from the results of the symbolic phase. The Cholesky factorization for symmetric positive definite systems has this property because it requires no numerical pivoting for stability, and the QR factorization for general nonsymmetric systems also has this property, but at considerable additional cost for storing the factors. Symmetric indefinite systems however in general require pivoting. This mixes numerical computation with irregular operations like atomically updating index sets or growing factor storage dynamically. Numerical pivoting may not seriously damage our static execution graph if it remains within a supernode’s diagonal block (this would have no impact on the sparsity pattern of our factors), but if it selects a 2x2 pivot outside of this block that entails a dynamic change to the factor’s sparsity pattern. If this happens repeatedly then it may break the static load balancing built up in the execution plan after the symbolic analysis phase, degrading performance.

I propose an application of the Woodbury matrix identity to help mitigate this. If a pivot would cause cross-supernode dynamic updates I instead threshold that value, with the threshold chosen to promote (but not necessarily guarantee) a stable factorization. Every such threshold is a rank-1 modification to our factor, in fact the factorization simply equals a diagonally shifted \(A\) where each diagonal entry is the threshold minus actual pivot

\[ A + diag(\delta) = LDL^T \]

Now we could simply accept this as-is and use it as a preconditioner, but usually it will require as many iterations as the rank change. This means if you required e.g. 20 thresholds during the factorization then solving \(Ax=b\) using \(LDL^T\) as a preconditioner theoretically has an upper bound of 20 iterations for e.g. GMRES, but in practice this upper bound is often reached exactly thus requiring 20 full scans of \(L\). For systems with large and expensive factors this is a very heavy cost. We can however completely recover the factor equation by observing

We can omit the zero columns of the diagonal shift through low rank form. Take \( E \) to be the column vectors associated with the modified pivots, and \( C \) to be the signed modifications. Then \( diag(\delta) = ECE^T \) and \( Y = L^{-1} E \). This gives

where solving with \(L,L^T\) proceeds as usual but the inner term requires a Woodbury update:

The factor \(Y\) involves solving a sparse triangular system against columns of the identity, we can omit storage of many nonzeros since we only need to store the nonzeros that appear along the elimination tree ancestor path extending from the thresholded value, so this is not a purely dense storage. We may also alternatively not precompute the Woodbury factors at the cost of needing to scan over the factors twice as many times per Krylov iteration, in the case where we are memory constrained and the number of thresholds is large (in the 100s) this could be a beneficial strategy.

This computation structure allows us to proceed with the factorization on its originally defined static schedule, and then have a final “clean up” phase which recovers the information lost from thresholding in the form of a low rank update. Solving systems involving our factorization now will require solving \(D-YCY^T\) which is where the Woodbury matrix identity will show up. I suspect this may be equivalent to or very similar to another form of pivot deferral where pivots do not get resolved immediately but rather get moved up the elimination tree to be handled at another stage in the factorization. The primary difference in structure is that with a Woodbury-deferred pivot no immediate pivot logic need be executed during the standard numeric phase, all of the Woodbury-deferred pivots instead get handled in a single batch operation post-factorization.

I will now look into some results approaching it this way

Results

For all results I do a supernodal sparse \( LDL^T \) factorization. For full pivoting I use the deferred pivoting approach, which moves a pivot further up the elimination tree instead of immediately resolving the dynamic pivoting decision (this still requires some pivot logic to execute during the numeric phase, in particular a future supernode must grow in size to accommodate the deferral)

I test three matrices

1. 3D Laplacian stencil with numerical entries sampled uniformly at random from [-1,1] and then symmetrized
2. 3D Helmholtz equation
3. An augmented system beginning with a random tridiagonal matrix A (not symmetric), but symmetrized into [[0,A^T],[A,0]].

And I compare four approaches to solving these systems, all of which used the same supernodal \(LDL^T\) factorization core with relaxed supernodes:

• Deferral - Classic deferred pivots
• Uncorrected - Thresholded pivots with no Woodbury correction
• WB - Woodbury-deferred pivots with precomputed \(Y\)
• WB on demand - Woodbury-deferred pivots without precomputed \( Y \)

Excess numerical storage / symbolic baseline

For each method and matrix I report the excess storage, that is the additional storage required to accomodate the pivoting logic as a ratio against a hypothetically un-pivoted system

and for absolute numbers (MiB).

GMRES iterations (effective L passes)

I ran unrestarted GMRES with all factorization approaches until convergence. I report the number of GMRES steps required for final convergence and in parentheses I also report the number of effective passes over \( L \).

Note that I’m not making a statement here that GMRES is the ideal solver for this situation, but it gives a better apples-to-apples comparison than running some with classic refinement and others with Krylov machinery.

and this resulted in the following backward errors post-convergence

Quantizing \( Y \)

While not really the focus of this blog post I did experiment with quantizing the \(Y\) factor to see if we could get the best of both worlds in terms of precomputing \(Y \), minimizing our extra storage, and minimizing the extra passes over \( L \) and I think there is potentially some merit to the idea. I just focused on the “augmented tridiagonal” system for these experiments as this one was specifically constructed to require a lot of pivoting

Relative excess storage:

Absolute storage:

GMRES convergence (effective L passes in parentheses)

Discussion

The key finding here is that we can recover direct-solver-like backward error while preserving the original static sparsity structure of the triangular factor. Woodbury-deferred pivots allow us to defer pivot logic until after full completion of the numeric phase and then follow up with calculating a Woodbury correction factor that represents the deferred pivots. This allows us to avoid mixing numeric logic with dynamic factor pattern updates, and the complete Woodbury correction may in principle be computed in a single batched kernel requiring only a single additional scan of the factors \(L,D\). This feels very similar to ordinary deferred pivoting (which is why I compared against it) with the exception that ordinary deferred pivoting still requires some symbolic update logic to execute during the numeric phase. Possibly “defer until root supernode” would be approximately equivalent to a Woodbury-deferral, but it still would require updating index sets because other columns along the elimination tree ancestor path may have nonzero intersection with the deferred pivot.

I did not time these results yet because the implementations were not optimized, I will have another post to share those results as they come up. I am also interested to see if the static nature of this approach may have benefits for GPU implementation, but this is speculative and I have not yet tried to do this.

Related work

I don’t exhaustively cover all possible literature but rather what seemed to be some key moments in related work. The earliest reference I could find that explicitly uses the Woodbury matrix formula to correct for thresholded pivots appeared in [1], but it was not implemented. The same thresholded pivot scheme is called “static pivoting” in [5], but there is no mention to the Woodbury updates and it is promoted rather as a preconditioner for GMRES and FGMRES. A more recent paper [3] implements a very similar scheme to what I proposed here but is instead a blocked algorithm intended for dense matrices. [2] operated on sparse matrices and allowed for supernode-local pivoting like I did here, and it mentions the possibility of a woodbury update but chooses instead to rely purely on Krylov refinement for convergence.

I don’t claim total novelty here as the approach has been referenced in passing in papers or applied in very related if not exact circumstances and considering that, it’s possible or even likely that this exact approach has been published and I simply did not find it in my literature survey.

That said what appears to be somewhat new here is the combination of static pivoting with supernodal-local pivots, and then cleaning them up with a subsequent Woodbury update. I also investigated different ways to mitigate the extra storage from lazy evaluation of the woodbury update (required 2x scans over \(L \)) to quantization.

1. Xiaoye S. Li and James W. Demmel (1998), “Making Sparse Gaussian Elimination Scalable by Static Pivoting,” SC ’98. Full paper

2. Olaf Schenk and Klaus Gärtner (2006), “On Fast Factorization Pivoting Methods for Sparse Symmetric Indefinite Systems,” ETNA 23, 158–179. Full paper

3. Neil Lindquist, Piotr Luszczek, and Jack Dongarra (2023), “Using Additive Modifications in LU Factorization Instead of Pivoting,” ICS ’23, 14–24. Full paper

4. Iain S. Duff and Stéphane Pralet (2007), “Towards Stable Mixed Pivoting Strategies for the Sequential and Parallel Solution of Sparse Symmetric Indefinite Systems,” SIAM J. Matrix Anal. Appl. 29(3), 1007–1024. Publisher page

5. Duff, Erisman, and Reid, Direct Methods for Sparse Matrices, Second edition, §13.4, “Static Pivoting”. Publisher link

6. S. Duff and J. K. Reid (1983), The Multifrontal Solution of Indefinite Sparse Symmetric Linear Equations ACM Transactions on Mathematical Software 9(3), 302–325 DOI