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.


2 Answers 2


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.

  • $\begingroup$ 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. $\endgroup$ Mar 26, 2021 at 6:32
  • 1
    $\begingroup$ 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. $\endgroup$ Mar 26, 2021 at 7:26
  • $\begingroup$ 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. $\endgroup$ Mar 26, 2021 at 8:11
  • $\begingroup$ I have a model problem, so I can make it as big as I want. The sky is the limit. My work mail is [email protected]. I will not use it for my current work due to time constrains, but will discuss it with my advisor for a later work. $\endgroup$ Mar 26, 2021 at 9:13
  • $\begingroup$ Also thanks for your time and answer. $\endgroup$ Mar 26, 2021 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.
  • 3
    $\begingroup$ 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. $\endgroup$ Mar 24, 2021 at 12:50
  • $\begingroup$ 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. $\endgroup$
    – Charlie S
    Mar 24, 2021 at 14:38
  • $\begingroup$ 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. $\endgroup$ Mar 24, 2021 at 16:35
  • 1
    $\begingroup$ 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. $\endgroup$ Mar 26, 2021 at 6:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.