Cholesky decomposition of sparse matrices using permutation matrices

Question

I am interested in the Cholesky decomposition of large sparse matrices. The problem I'm having is that the Cholesky factors are not necessarily sparse (just like the product of two sparse matrices is not necessarily sparse).

For example for a matrix with non-zeros only along the first row, first column, and diagonal the Cholesky factors have 100% fill-in (the lower and upper triangles are 100% dense). In the image below the gray is non zero and the white is zero.

One solution I'm aware is to find a permutation P matrix and do the Cholesky decomposition of P^TAP. For example with the same matrix by applying a permutation matrix which moves the first row to the last row and the first column to the last column the Cholesky factors are sparse.

My question is how to determine P in general?

To get an idea of the difference of the Cholesky decomposition of A and P^TAP from a more realistic matrix see the image below. I took all these images from http://www.seas.ucla.edu/~vandenbe/103/lectures/chol.pdf

According to the lecture notes

many heuristic methods (that we don’t cover) exist for selecting good permutation matrices P.

I would like to know what some of these methods are (code in C, C++, or even Java would be ideal).

Crucible

Nice people put whole books on the topic online: Y. Saad: Iterative methods for sparse linear systems, esp. chapter 3. Try "sparse reordering" on google, possibly with "hypergraph" as additional keyword. – Orthopedics

@LutzL, do you know of any libraries which do this? E.g. Eigen or MKL? Or one for Java? – Crucible

[Patrick R. Amestoy et al.: Linear algebra and sparse direct methods]( amestoy.perso.enseeiht.fr/COURS/unesco2004.pdf) chapter "Reordering sparse matrices", and other text on his page amestoy.perso.enseeiht.fr – Orthopedics

libsparse, which is now "SuiteSparse: UMFPACK, CHOLMOD, and many other". Also, there is netlib.org/sparse which is an ancestor or cousin,... – Orthopedics

Find many fascinating images (and very often the theory behind them) with google imagesearch google.de/search?tbm=isch&q=Reordering+sparse+matrices – Orthopedics

@LutzL, thanks for all the links! Reading through some of them I found en.wikipedia.org/wiki/Minimum_degree_algorithm. I think that's along the lines of what I'm looking for. Something like symamd. – Crucible

Is the matrix symmetric? Otherwise Chelesky decomposition does not apply. Also a triangular matrix is not dense, but 50% sparse. – Dilapidation

I'm voting to close this question as off-topic because it is primarily about Mathematics – Dilapidation

@ja72, of course the matrix is symmetric! It's also positive definitive. But I have even created a modified Cholesky decomposition using the LDL form for matrices which are ill-defined (which happens sometimes in my project). I'm quite familiar with Cholesky decomposition of dense matrices. Now I'm learning about sparse algorithms now. It's really a optimization issue (and memory). – Crucible

@ja72, your comment about triangular matrices being 50% sparse is a bit pedantic. I should have said 100% fill-in but I did not want to define fill-in. L+LT is 100% dense. – Crucible

@ja72, there is no right answer for P. It's NP-complete. There are a bunch of heuristic methods to find a P. It's a optimization question (both in time and memory). It's analogous to using e.g. a kd-tree in ray tracing to speed it up. Finding a suitable P is necessary for me to scale the research project I'm involved in up to a real world case. Do you have any useful comments on an algorithm to determine P? – Crucible

I think the ij-th element of L depends only on the values to the left and up in A. I mean there is a tree structure (or hierarchy) to the values and moving the more dense rows towards the bottom (and columns to the right) might better diagonalize the starting elements of L – Dilapidation

Maybe you can provide us with a simple 6×6 symmetric matrix you want triangulated (decomposed) so we can try some ideas. – Dilapidation

@ja72, it's easy to construct a symmetric positive definite matrix. Loop over one of the triangles (but not the diagonal) and write the same value to the transpose. Then set the diagonal values to n for an nxn matrix. See cholesky-decomposition-of-large-sparse-matrices-in-java – Crucible

That is the job of the OP. See stackoverflow.com/help/mcve – Dilapidation

Answer 1

The problem of finding an optimal permutation of rows and columns of a matrix for a minimum fill-in matrix-factorization is not a trivial trask (as pointed out in the comments). Thus, heuristic algorithms are used in praxis.

There are some libraries that implement heuristic renumbering/ordering-strategies, often based on graph-algorithms for the adjacency-graph of your matrix. One tries to reduce the bandwidth of the corresponding adjacency-matrix. An easy to implement algroithms is the Cuthill-McKee Algorithm or the Minimum-Degree Ordering algorithm. More about this problem can be found in the Book Yousef Saad: Iterative Methods for Sparse Linear Systems (2003), upon many others.

Many libraries implement heuristic algorithms, like

• Suitesparse A collection of libraries for direct solvers for largse sparse linear systems. Ordering methods implemented in the libraries AMD, CAMD, COLAMD, and CCOLAMD
• (Par-)Metis A library for Graph-partitioning, but provides Matrix reordering algorithms as well
• Boost.Graph Working on the adjacency graph directly and provides some ordering algorithms, like the mentioned Cuthill-McKee, and Minimum-Degree Ordering
• (PT-)Scotch for Graph-partitioning and sparse-matrix reordering

Some of these libraries provide also sparse Cholesky factorization methods and can be used directly.

Kimball

Sorry for the delayed acceptance. I have put this project on hold for a bit and will get back to later this year. – Crucible

In case you're interested I "solved" my problem in the end by using Eigen. – Crucible