Thread Subject: sparse matrix multiplication and cgs()

Hello,
   I solved Ax=b, A is large(of size 2e5), sparse(each row has no more than 26 nonzero), symmetric positive definite matrix. I did it in two ways,
1) L=chol(A,'lower'); tic, x=L\(L'\b),toc;
2) tic,[x,flag,resrel,iter]=cgs(A,b,[],n,[],[],x0),toc;

L is dense, iter=6. So theoretical 2) should be must faster than 1). But 1), 2) ends up with the same time.
Why?

GS <gshy2014@gamil.com> wrote in message <17247123.13121.1243869597327.JavaMail.jakarta@nitrogen.mathforum.org>...
> Hello,
> I solved Ax=b, A is large(of size 2e5), sparse(each row has no more than 26 nonzero), symmetric positive definite matrix. I did it in two ways,
> 1) L=chol(A,'lower'); tic, x=L\(L'\b),toc;
> 2) tic,[x,flag,resrel,iter]=cgs(A,b,[],n,[],[],x0),toc;
>
> L is dense, iter=6. So theoretical 2) should be must faster than 1).

Why not using PCG, which is design for spd matrix? Not sure about the "theoretical" claim above: CGS/PCG convergence depends on the condition number of A. If cond(A)=A it will converge in no time (1 iteration is enough in fact).

Iterative methods are not always faster than direct method.

Bruno

GS <gshy2014@gamil.com> wrote in message <17247123.13121.1243869597327.JavaMail.jakarta@nitrogen.mathforum.org>...
> Hello,
> I solved Ax=b, A is large(of size 2e5), sparse(each row has no more than 26 nonzero), symmetric positive definite matrix. I did it in two ways,
> 1) L=chol(A,'lower'); tic, x=L\(L'\b),toc;
> 2) tic,[x,flag,resrel,iter]=cgs(A,b,[],n,[],[],x0),toc;
>
> L is dense, iter=6. So theoretical 2) should be must faster than 1). But 1), 2) ends up with the same time.
> Why?

chol will be much faster if you allow it to permute the matrix (and L will not be dense). Try:

tic ; [L g q] = chol (A, 'lower') ; toc
tic ; x = q * (L' \ (L \ (q' * b))) ; toc

See also the FACTORIZE object, here:
http://www.mathworks.com/matlabcentral/fileexchange/24119

which does this for you. You can speed the forward/backsolve even further by using cs_lsolve for the L\(...) part, and cs_ltsolve for the L'\(...) part. See CSparse on my web page at http://www.cise.ufl.edu/research/sparse. The L'\(...) solve in MATLAB forms the explicit L transpose, which is slow. Alternatively, you could just compute LT = L' first, and then do the backslashes.

x=A\b does all this without forming L transpose explicitly.

THe other problem you have, which is more fundamental, is that your test is printing out the solution x (you have a comma instead of a semicolon). That dominates the run time, so it's not surprising that both methods take the same amount of time.

You can speed up chol even more by using a better ordering method that isn't in MATLAB (METIS). CHOLMOD has an interface to it (in SuiteSparse, at the same web page).

In any case, a problem of dimension 2e5 is probably small enough that a direct method will typically be faster ... at least if you allow it to do a fill-reducing ordering.

> Why not using PCG, which is design for spd matrix?
pcg is slower than cgs. pcg() needs about 16 iterations.

> Not sure about the "theoretical" claim above:
I mean the floating point operation count. Since there are no more than 26 nozeros each row/column, we only need O(n) (where n=size(A,1)) additions/multiplications for each iteratios. since we only need no more than 10 iterations in total. The total number of addition/multiplications is O(n). But for L\y, we need 2n^2 additions/multiplications.

>CGS/PCG convergence depends on the condition number of >A. If cond(A)=A it will converge in no time (1 >iteration is enough in fact).

cond(A) is not too bad.(about 1e3). cond(A) will affect the convergence, but it turns out that the it only needs no more than 6 iterations.

>
> Iterative methods are not always faster than direct
> method.
>
> Bruno