I have a very large sparse matrix, 'A', in Ax = b. What work in advance of getting 'b' can be done to reduce solving time?

Question

This question borders between a programming and math question (more math). I have a little matrix knowledge but this is past my ability, so any help is very much appreciated.

Question

I have a very large sparse matrix (say one million rows square). I'd like to solve this for many different sets of constant values (to clarify, by constant values, I mean the 1xn matrix 'b' in Ax = b). Performance is absolutely crucial - so it would make sense to do whatever simplification I can in advance of introducing each set of constants.

Goal

I'd like to get the time complexity for solving each constant set down to O(n)... I may be dreaming however.

Other Info

• The matrix really is a million (maybe even 100 million) rows
• It's very sparse
• It has a bandwidth ranging between 2 and as much as a few hundred
• Even for larger bandwidths, each row has max. 7 or 8 non zero entries
• It takes an (extremely rough) diagonal form
• Each matrix 'b' in Ax = b is actually very sparse itself
• Due to the nature of the problem, a solution should always be available
• All nonzero values are real, floating point numbers. They can be negative however.

Answer 1

Instead of reinventing the wheel, look at the major linear algebra packages. Many have procedures for large sparse matrices. One method is "sparse $LU$ decomposition" (google that phrase for tutorials on it). Once the matrix is factored, solving multiple equations with different right sides becomes fast. The difficult part is factoring into sparse factors, which is tricky, and that's why it would be best to use existing code rather than write your own. It is theoretically possible for a sparse matrix to be impossible to factor into sparse factors (the additional non-zero entries are called "fill-in") but there are heuristics that attempt to keep the factors sparse and in practice they are often successful. The only sure way to know if it will work is to try it.

Note that if your matrix is symmetric there will be special versions for that case.

Answer 2

A system of linear equations with million or more independent variables, even with a sparse matrix, justifies attempting an iterative approach. If the matrix is symmetric and positive definite, the Conjugate Gradients Method is the way to go. LU-factorization, being attractive, brings on the problem of matrix fill-in: any linear solver based on Gaussian approach (including LU-factorization) will have to store the matrix on the disk. Performance, as a result, will be significantly slowed down. However, since the modern processors are multi-core, the iterative methods can be applied simultaneously to several right-hand sides, thus improving the performance significantly; and the whole matrix will be stored completely in the Random-Access Memory. Here is one such program, allowing solution of systems with multiple right-hand sides using the iterative approach: http://members.ozemail.com.au/~comecau/CMA_Sparse.htm

Answer 3

There is a huge amount of research on this problem --- it turns out that solving huge linear systems is pretty ubiquitous in applied maths, because linear algebra is one of the few things that we can do efficiently for large problems.

Sparse LU decompositions have already been mentioned; the other main approach is Krylov subspace methods (such as GMRES) with a vast choice of preconditioners --- there is a whole lot of strategies, and many of them work only for specific problems.

I would suggest contacting someone in person; see who works in numerical analysis at your local institution. $10^6$ to $10^8$ is already quite a large problem, and it is likely that some tricks of the trade are needed to get things working. All this especially if you say that "performance is crucial".

Also, in any case try to use an existing program or library (Matlab, for some quick experimentation if the problem fits in memory for it, or Scipy, PetSC, Trilinos, ...). If you are not an expert, rewriting things is unlikely to match state-of-the-art code in terms of performance.