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

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.

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.

• That is an interesting idea. Are there effective parallel implementations of it? The part of the problem that I didn't mention is that I am using LU from LAPACK (MKL and PLASMA, etc). And it has very good parallel speed-up. – Dimitar Slavchev Mar 26 at 6:32
• Quick experiment on my laptop, a single LU of matrices of size 4000 takes about 0.39s. A HT reduction takes about 29s. That's almost 100 times slower. This is likely to get worse for larger matrices as the LU scales a bit better. – Thijs Steel Mar 26 at 7:26
• If you end up trying this idea, i have a pretty fast implementation of the HT reduction on GPU. The paper is not ready for publication, but i can provide you a copy via email if you want. – Thijs Steel Mar 26 at 8:11
• I have a model problem, so I can make it as big as I want. The sky is the limit. My work mail is dimitargslavchev@parallel.bas.bg. I will not use it for my current work due to time constrains, but will discuss it with my advisor for a later work. – Dimitar Slavchev Mar 26 at 9:13
• Also thanks for your time and answer. – Dimitar Slavchev Mar 26 at 9:38

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.
• You might be able to use the LU decomposition as a preconditioner in an iterative solution of the full system, allowing you to solve with just a few iterations. It might work for a few time steps and then you may need to recompute the LU factorization. – Amit Hochman Mar 24 at 12:50
• I agree with your assessment that the general full update of the LU decomposition is not worth the effort. However, if $M$ is diagonal, then it might be worth investigating a series of rank-1 updates to your decomposition. I suspect even that is not worth the headache involved, and may introduce undesirable numerical side effects since you intend to update your factorization many times. – Charlie S Mar 24 at 14:38
• Yes, this might be the most asked question here. There is no known way to make full-rank updates of LU decomposition, not even if you just want to add a multiple of the identity. – Federico Poloni Mar 24 at 16:35
• I am making a comparison between using Gaussian solver and HSS compression (from STRUMPACK). And with HSS compression, I can just use the compressed matrix $H$ over and over again. Well for the lumped matrix of mass at least, since the compression doesn't touch the diagonal. – Dimitar Slavchev Mar 26 at 6:39