solving sparse linear system
Peter
Given a sparse invertible matrix A and a vector b. As far as I know
there are two types of efficient ways to solve Ax=b: gradient-based
iterative methods (such as conjugate gradient, Gauss-Seidel, Jacobi
method, etc) and direct sparse inversion (such as sparse LU, SVD,
Gauss elimination, etc). Of course there are also other types of
methods, but these two are most favored in practical applications.
I know it depends very much on the specific problem, but is there some
simple rule of thumb which determines which one is computationally
more effective?
Any help, information, keyword, idea, book title etc. would be
appreciated
Peter
bdb
http://en.wikipedia.org/wiki/Positive-definite_matrix
Unfortunately, there is no rule, it depends on your system but for very
larges systems iterative method are recommended.
My opinion : it is always better to build the system which facilitate the
system solving. For that you need a good understanding of the physics
behind the matrix.
Peter
> The gradient based method requires a positive definite matrixhttp://en.wikipedia.org/wiki/Positive-definite_matrix
> Unfortunately, there is no rule, it depends on your system but for very
> larges systems iterative method are recommended.
> My opinion : it is always better to build the system which facilitate the
> system solving. For that you need a good understanding of the physics
> behind the matrix.
Thanks for the answer.
I give two examples below. So why is conjugate gradient inefficient in
the #1st case and why does it work well in the #2nd?
Example#1:
Laplacian U(x,y) + k^2 n(x,y)^2 U(x,y) = S(x,y)
where U is the wave field, Laplacian is the Laplacian operator, n(x,y)
is the index of refraction, k is the constant wavenumber and S(x,y) is
the source term. This is a sparse linear system. By some finite
difference of finite element scheme it can be written as:
L U = S
where L is a sparse self adjoint matrix, U is the unknown vector we
would like to solve for, S are the known source terms. In this case
without preconditioning conjugate gradient converges very slowly,
althoug it is sparse and large scale and positive definite. And sparse
LU or some other direct method is much more efficient.
Example #2:
non-uniform Fourier transform of a function can be approximated as a
uniform fft on grid and an interpolation in the non-uniform points:
Up FFT f = F
where f is the function, FFT is the 2d fft operator on grid, Up is the
interplation operator and F is the result. Both operators are large
scale and their matrices are sparse.
In this case conjugate gradient on the normal equations is very
effective to compute the inverse and any direct method is impossible
to implement.
Best regards
Peter