Efficient solution of a linear system with a diagonal Laplace matrix +
Dear 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 - the sum of the Laplace matrix L and the diagonal matrix D
- Laplace matrix L is sparse, about 25 non-zeros per row
- the system is large, with as many unknowns as there are pixels in the input image (usually> 1 million).
The Laplace matrix L does not change between successive runs of the algorithm; I can plot this matrix on preprocessing and possibly calculate its factorization. The diagonal matrix D and the right vector b change with each run of the algorithm.
I am trying to figure out what would be the fastest method to solve the system at runtime; I don't mind wasting time on preprocessing (e.g. calculating the factorization L).
My initial idea was to pre-compute the Cholesky L factorization and then update the runtime factorization with the values from D (update rank 1 with cholupdate) and quickly solve the back substitution problem. Unfortunately, the Cholesky factorization is not as rare as the original matrix L, but simply loading it from disk takes 5.48 s; as a comparison, it takes 8.30s to directly solve the backslash system.
Given the shape of my matrices, is there any other method you would recommend to speed up the solution at runtime, no matter how long it takes for preprocessing time?
Thanks in advance for your help!
source to share
Assuming you are working on a grid (since you mention images, although this is not guaranteed), you are more interested in speed than accuracy (since 5 seconds seems too slow for 1 million unknowns), I see several options.
First, forget about exact methods like Cholesky (+ reordering). Even if they allow you to store factorization and reuse it across multiple rhs, you will most likely have to store giant matrices that seem unsolvable in your case (hopefully you are reordering rows / columns with reverse Cuthill McKee or something else though - this greatly simplifies factorization).
Depending on your boundary conditions, I would first try Matlab poisolv
, which solves the Poisson problem using FFT and possible reprogramming if you need Dirichlet boundary conditions instead of periodic ones. This is very fast, but may not be suitable for your problem (you mentioned having 25 nnz for a Laplacian + identity matrix: why? Is this a high order Laplace matrix, in which case you might be more interested in accuracy than I assume? Or is - is this actually a different problem than the one you are describing?).
Then you can try multigrid solvers which are very fast for imaging and smooth problems. You can use a simple relaxation technique for each iteration and each level of a multigrid program, or use more convenient techniques (such as a preset conjugate gradient level). Alternatively, you can do a simpler preset conjugate gradient (or even SSOR) without the multiseries, and if you are only interested in an approximate solution, you can stop iterations until convergence is complete.
My arguments for iterative solvers:
- you can stop before convergence if you want an approximate problem
- you can reuse other results to initialize your solution (for example, if your different runs correspond to different video frames, then using the solution of the previous frame as the initialization of the next one makes sense).
Of course, a forward solver, for which you can precompile, store, and store the factorization, also makes sense (although I don't understand your argument for updating rank-1 if your matrix is constant), since only back substitution is left to be done at runtime. But this ignores the structure of the problem (regular grid, possible interest in limited precision results, etc.), I would choose methods that have been developed for cases like Fourier-like methods or multirid. Both methods can be implemented on the GPU for faster results (remember, GPUs are adapted for image / texture processing!).
Finally, you can get some interesting answers from scicomp.stackexchange , which is more numerical oriented.
source to share