Efficiently solving linear system with Laplacian + diagonal matrix

Dear fellow Matlab users,

In my implementation of an image processing algorithm, I have to solve a large linear system of the form A * x = b, where:

• matrix A = L+D is the sum of a Laplacian matrix L and a diagonal matrix D
• Laplacian matrix L is sparse, with about 25 non-zeros per row
• the system is large, with as many unknowns as there are pixels in the input image (typically > 1 million).

The Laplacian matrix L does not change between successive runs of the algorithm; I can construct this matrix in preprocessing, and possibly compute its factorization. The diagonal matrix D and right-side vector b change at each run of the algorithm.

I am trying to find out what would be the fastest method to solve the system at runtime; I do not mind spending time on preprocessing (for computing a factorization of L, for example).

My initial idea was to pre-compute a Cholesky factorization of L, then update the factorization at runtime with values from D (rank-1 update with cholupdate), and solve quickly the problem with back-substitution. Unfortunately, the Cholesky factorization is not as sparse as the original L matrix, and just loading it from disk already takes 5.48s; as a comparison, it takes 8.30s to directly solve the system with backslash.

Given the shape of my matrices, is there any other method that you would recommend to speedup the solving at runtime, no matter how long it takes at preprocessing time?

Thanks in advance for your help!