Is there a simple way to add a sparse matrix to an LU decomposition of a dense matrix?

Question

I am solving a parabolic equation in the form:

$$ \left( {M \over\tau_j} + A \right) u^{j+1} = M f^j + {u^j \over \tau_j}, $$

where $A$ and $f$ are a dense stiffness matrix and the right hand side of the fractional diffusion problem from this paper:

[Acosta, Gabriel & Bersetche, Francisco & Borthagaray, Juan. (2017). A short FE implementation for a 2d homogeneous Dirichlet problem of a Fractional Laplacian. Computers & Mathematics with Applications. 784-816. 10.1016/j.camwa.2017.05.026. ]

$M$ is the matrix of mass, $u^j$ is the solution at a previous time step and $u^{j+1}$ is the vector of unknowns. $\tau_j$ is a number.

I have previously solved the non-parabolic problem $Au = B$ with LU decomposition through lapack's dgesv function.

$M$ is a sparse matrix with around 5 elements per row. It could also be reduced to a diagonal matrix (lumping of mass) by summing all coefficients to the diagonal element.

So my question is: is it possible to calculate $A = LU$ once and then only update the $LU$ factorization with the $M \over \sigma_j$ for each step?

I have a vague memory that I found a lecture notes or a textbook online about it a few months ago. Alas my google-fu is failing me right now. I tried googling for "parabolic equations solution with LU" and a few other variants about parabolic equations.

Answer 1

There is another factorization you could consider: the Hessenberg upper-triangular reduction. It's usually used as a preprocessing step in the QZ algorithm, but it has other uses as well.

Consider the reduction $(A,M) = Q^T(H,T)Z$. Then $\frac{M}{\mu_j} + A = Q^T(\frac{T}{\mu_j} + H)Z$, which is relatively cheap to solve as it only involves orthogonal matrices and a Hessenberg matrix. This is valid for any value of $\mu_j$, so long as the matrices $A$ and $M$ don't change.

However, the Hessenberg upper-triangular reduction is pretty expensive. This is only going to be worth it if you have a whole lot of systems to solve.

Answer 2

And just a few minutes after asking I found an answer. The procedure above is called "updating LU". This question has a nice generic answer with links to other more specific questions.

full rank update to cholesky decomposition

So the short answer is no, in my case. You cannot update LU decomposition with a full rank matrix in less than $O(n^3)$, which is the same as performing the decomposition again if the update is full rank. And the matrix of mass is full rank, I believe.

EDIT: Another option is going around the problem by replacing the LU decomposition with something else. For example:

1. One can use Hessenberg upper-triangular reduction as advised by Thijs Steel in another answer
2. With Hierarchical methods to compress the matrix $A$ and then do updates on the compressed matrix before factorization. Provided that $A$ is suitable for compression. I do use Hierarchical Semo-separable compression from STRUMPACK for that, but this question was to see if there is a way to "optimize" the direct gaussian solver.
3. Or other options that I am not aware of.