Symbolic factorization of a large sparse matrix

I'm solving a series of problems Ax=b, where A changes, but still has the same sparsity pattern. To exploit that, in the initialization, I only perform symbolic factoriazation once.

Currently, I use pardiso:

https://www.mathworks.com/matlabcentral/answers/497401-pardiso-reuse-symbolic-factorization

I was wondering how come matlab doesn't offer something like this, considering that the mldivide (based on Tim Davis's SuiteSparse) is supposed to be the fastest option?

Answer 1

[Zohar]

1. When a symmetric matrix (Hessian) is involved, MA57 is used (ldl), which is faster than suitesparse.
2. The symbolic factorization is the reordering that produces sparser factors. It's done by default with amd (metis is another option).

More details:

https://www.youtube.com/watch?v=1dGRTOwBkQs&ab_channel=TimDavis

From [Boyd09, app. C.3]:

"For the sparse Cholesky factorization, the re-ordering permutation P is often determined using only sparsity pattern of the matrix A, and not the particular numerical values of the nonzero elements of A. Once P is chosen, we can also determine the sparsity pattern of L without knowing the numerical values of the nonzero entries of A. These two steps combined are called the symbolic factorization of A, and form the first step in a sparse Cholesky factorization. In contrast, the permutation matrices in a sparse LU factorization do depend on the numerical values in A, in addition to its sparsity pattern.

The symbolic factorization is then followed by the numerical factorization, i.e., the calculation of the nonzero elements of L. Software packages for sparse Cholesky factorization often include separate routines for the symbolic and the numerical factorization. This is useful in many applications, because the cost of the symbolic factorization is significant, and often comparable to the numerical factorization."

Watch for the following pitfall: when using analyzePattern (style Eigen Pardiso wrapper), the nonzero pattern needs to be preserved throughout the iterations. Matlab by default (in certain common operations) prunes nonzeros that are numerically zeros thus upsetting the pattern. I used to have unexplained crashes with Pardiso matlab wrapper.