# Optimizing Cholesky factorization for multiple sparse matrices with same nonzero pattern

I'm using a Cholesky factorization to solve the linear step in a nonlinear system of equations (nonlinear finite element analysis). In the PETSc library, one can specify a parameter for SAME_NONZERO_PATTERN during successive solves.

The speedup is good, almost too good (100 times faster for the second solve). This makes me wonder, what sort of optimization one can do for Cholesky factorizations when solving for many sparse matrices with same nonzero pattern (but with different, possibly similar, values)?

-
Perhaps it is a related paper cs.ucdavis.edu/~bai/Winter09/nakatsukasabaigygi09.pdf –  Sunni Aug 14 '11 at 18:33

Almost all sparse direct solver packages split the factorization into two stages:

• symbolic factorization computes an ordering, often using nested dissection or a approximate minimum degree algorithm, such that the factors will be as sparse as possible and allocates space to hold the result. The value of the matrix entries are either not used or only used to make estimates about pivoting.

• numeric factorization computes the factorization given the ordering and sparsity computed by symbolic factorization.

In most circumstances, symbolic factorization is much less expensive than numeric factorization, so SAME_NONZERO_PATTERN offers limited benefit. This changes when running in parallel with many processes because the symbolic factorization does not scale as well (and some popular packages, including MUMPS, compute it in serial).

It is very unlikely that your factor of 100 comes from the symbolic factorization stage. Perhaps you are actually reusing the factors from the previous iteration? For performance questions like this, it usually helps to run with -log_summary and compare times and load balance for the Mat*FactorSym and Mat*FactorNum events. Send the output to petsc-users@mcs.anl.gov or petsc-maint@mcs.anl.gov if you would like help interpreting the output.

As for "updating" the factors after small changes in the matrix entries, attempts to do this have typically been unsuccessful. An alternative which I recommend is to use the old factorization as a preconditioner for the new matrix. For nonlinear problems using SNES, you can experiment using -snes_lag_preconditioner or -snes_lag_jacobian. If you use KSP directly, then pass SAME_PRECONDITIONER to KSPSetOperators().

-

I don't have much experience with sparse factorizations and am not sure if I understand your question correctly. But to my knowledge the sparse factorization methods first compute a rearrangement of the matrix elements so that the result of the factorization (the triangular matrix $\mathbf{L}$ for which $\mathbf{A}=\mathbf{LL}^T$, or $\mathbf{L}^T\mathbf{L}$) has a very sparse structure itself. So I assume the simple optimization is, that the rearrangement doesn't have to be recomputed each time and can simply be reused for the same sparsity pattern of $\mathbf{A}$, as it should be independent of the actual values.

-

A nice example is the case of smoothing signals using cubic splines. In that case, positive-definite systems of equations arise that have rich structure, which allows for the Cholesky factorization to be computed in linear complexity (instead of cubic complexity). In the example of cubic splines, the coefficient matrix $A$ is pentadiagonal, symmetric, persymmetric and Toeplitz. Then considering the equation of the Cholesky factorization $A = L L^T$, where say $L$ is lower triangular, it can be shown that $L$ has a very rich structure as well and that its elements can be computed by solving a coupled system of three second order difference equations. However, the amenability of a positive-definite sparse matrix to optimization of its Cholesky factorization, greatly depends on the particular structure of the sparsity.

-
Thank you for your answer, but my use case doesn't assume any particular structure, but rather the subsequent factorization of more matrices with the same nonzero pattern (where only the second factorization is speed up significantly). I suppose I'm looking for rearrangements operations to make the matrix banded, which could be applied to the same nonzero structure several times, or perhaps something else along these lines. –  Mikael Öhman Jul 18 '11 at 23:25