Performance of chol() on sparse matrices

[David Bateman]

Andreas Stahel wrote:
> Dear Octave users
> 
> when testing code for a class a timing result on Octave 3.0.3 puzzled me.
> generate a very sparse, symmetric, positive definite matrix Anxny (size 
> 62500x62500) and time a few commands
> 
> x=Anxny\b  -> 0.8 sec
> R=chol(Anxny) -> 7.3 sec
> x=R\(R'\b)  -> 2.3 sec
> [L,U,P]=splu(Anxny) -> 12 sec
> 
> I would expect the Cholesky back-substitution to be fastest and 
> cho(Anxny) to be comparable to Anxny\b !!
> 
> Would you happen to have hints on why this occurs

Try instead

[R, P, q] = chol (Anxny,'vector');
x(q) = R \ (R' \ b(q));

without the Q return value, chol can't use the sparsity preserving 
column transformations.

With the script below I see

SolveTime =  0.19601
CholTime =  0.66004
CholSolveTime =  0.30802
Chol2Time =  0.20801
Chol2SolveTime =  0.028001
LUTime =  1.3121

and norm (x4.'-x3) equal to  1.5781e-11


D.


nx=100; ny=100;
hx=1/(nx+1); hy=1/(ny+1);
Dxx=spdiags([-ones(nx,1) 2*ones(nx,1) -ones(nx,1)],[-1 0 1],nx,nx)/(hx);
Dyy=spdiags([-ones(ny,1) 2*ones(ny,1) -ones(ny,1)],[-1 0 1],ny,ny)/(hy);
Anxny=kron(speye(ny),Dxx) + kron(Dyy,speye(nx));

b=ones(nx*ny,1);
t0=cputime;
x2=Anxny\b;
SolveTime=cputime()-t0
t0=cputime;
R=chol(Anxny);
CholTime=cputime()-t0
t0=cputime;
x3=R\(R'\b);
CholSolveTime=cputime()-t0
t0=cputime;
[R,P,q] = chol(Anxny,'vector');
Chol2Time=cputime()-t0
t0=cputime;
x4(q) = R\(R'\b(q));
Chol2SolveTime=cputime()-t0
t0=cputime;
[L,U,P]=splu(Anxny);
LUTime=cputime()-t0
-- 
David Bateman                                dbateman at dbateman.org
35 rue Gambetta                              +33 1 46 04 02 18 (Home)
92100 Boulogne-Billancourt FRANCE            +33 6 72 01 06 33 (Mob)