3
$\begingroup$

I asked the same question on mo.se and it was suggested that scicomp would be a better forum for it. So here it is:

I am writing code for solving linear equations of the form

$$A_{n\times n}\cdot x=1_n$$

where $n$ is on the order of $10^6$ and $A$ is a symmetric matrix with approx $10^3$ nonzero entries in each row. This makes its size barely manageable, but inverting it is infeasible, and I'm not sure which decomposition suitable for solving linear equations would lead to sparse matrices.

Thus the questions:

1) Any hope that $LU$ decomposition of a symmetric sparse matrix would be sparse?

2) Is it possible to take advantage that r.h.s. is a scalar to simplify solution of the above equation?

3) What would be the best numerically stable algorithm to handle linear equations of that size?

Improve this question
$\endgroup$
4
  • 1
    $\begingroup$ Is there any structure to $A$ (such as banded)? Do you have some idea how the eigenvalues are distributed? Have you tried a Krylov method such as conjugate gradients, or is there a reason you need a direct method? $\endgroup$
    – Christian Clason
    Commented May 9, 2014 at 23:33
  • $\begingroup$ Unless your matrix has a very specific structure, you're almost certainly better off with an iterative method. $\endgroup$
    – Doug Lipinski
    Commented May 10, 2014 at 0:52
  • $\begingroup$ @ChristianClason, $A$ consists of something like $e^{-d_{ij}^2}$, where $d_{ij}$ are distances between 2-D objects. If these were 1-D objects $A$ would be banded; in 2-D the rows can be arranged in one of the 2-D direction to make the matrix almost multi-banded, with tiny values that we consider zeros between bands. It's not perfectly multi-banded because the 2-D items are not uniformly distributed, so the multiple bands won't be perfectly aligned. $\endgroup$
    – Michael
    Commented May 12, 2014 at 17:17
  • $\begingroup$ A direct solver would definitely appreciate a rearrangement to reduce banding (something like Cuthill-McKee), but if $A$ itself is barely manageable, the factors are very unlikely to be. Your best bet is still a Krylov solver (CG if $A$ is positive definite, MINRES if not). $\endgroup$
    – Christian Clason
    Commented May 12, 2014 at 18:18

1 Answer 1

Reset to default
5
$\begingroup$
  1. In your case, almost no hope. Only very specific types of sparse matrices have a sparse $LU$ decomposition. For example, (trivially) diagonal and tridiagonal matrices.
  2. Not really. What you're basically saying is you want the sum of all the columns of $A^{-1}$. If instead the RHS was sparse, you could use selective inversion (SelInv) techniques, potentially.
  3. Your matrix is real symmetric I presume, so you can use an iterative method like conjugate gradients (CG). All it requires is that you can apply the matrix-vector product onto an arbitrary matrix. You may need a preconditioner, or to apply something like algebraic/geometric multigrid techniques. If you know something about the low rank off diagonal structure of $A$, you could consider hierarchical matrix factorization methods (this is pretty advanced and there is very little code out there for it). Finally, your matrices are not that sparse. You could consider treating them densely and using an out-of-core sparse LU. This is really more of a last resort if nothing else works.

If you know that $A$ is positive definite, then replace LU above with Cholesky, which is slightly more efficient.

Improve this answer
$\endgroup$
1
  • $\begingroup$ Sparse LU factorization is a widely used technology. It's not perfect, but it can work well in many situations. This particular problem is quite large for sparse direct methods though. Iterative methods such as GMRES should be applicable. However, unless A is symmetric and positive definite, CG isn't applicable. $\endgroup$
    – Brian Borchers
    Commented May 10, 2014 at 3:41

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.