Thread Subject: Solving Ax=b for large sparse matrix(1e+6 X 1e+5) and cholinc failure

I have to solve Ax=b where A is a large sparse matrix (1e+6
x 1e+5) with about

7e+6 non-zero elements.

 

I tried to solve A'Ax=A'b using pcg, but the new matrix
(A'A) is very bad conditioned. I tried to use cholinc to
precondition the system, but A’A has about 1e+7 non-zero
elements, and Matlab gives me out of memory error.

 

I tried also to solve the system using:

x=(A’A)\(A’b)

 

but here it is what I get:

 

sp\: bandwidth = 5864+1+5864.

sp\: is A diagonal? no.

sp\: is band density (0.01) > bandden (0.50) to try banded
solver? no.

sp\: is A triangular? no.

sp\: is A morally triangular? no.

sp\: is A a candidate for Cholesky (symmetric, real positive
diagonal)? yes.

sp\: is CHOLMOD's symbolic Cholesky factorization (with
automatic reordering) successful? yes.

sp\: is CHOLMOD's numeric Cholesky factorization successful? no.

 

CHOLMOD version 1.1.0, May 5, 2006: : status: ERROR, out of
memory

  Architecture: Linux

    sizeof(int): 4

    sizeof(UF_long): 4

    sizeof(void *): 4

    sizeof(double): 8

    sizeof(Int): 4 (CHOLMOD's basic integer)

    sizeof(BLAS_INT): 4 (integer used in the BLAS)

  Results from most recent analysis:

    Cholesky flop count: 8.0627e+11

    Nonzeros in L: 1.9098e+08

  memory blocks in use: 11

  memory in use (MB): 8.8

  peak memory usage (MB): 1961.0

  maxrank: update/downdate rank: 8

  supernodal control: 1 40 (supernodal if flops/lnz >= 40)

  nmethods=0: default strategy: Try user permutation if
given. Try AMD.

    Select best ordering tried.

    method 0: user permutation (if given)

    method 1: AMD (or COLAMD if factorizing AA')

        prune_dense: for pruning dense nodes: 10

        a dense node has degree >= max(16,(10)*sqrt(n))

        flop count: 8.0627e+11

        nnz(L): 1.9098e+08

  final_asis: TRUE, leave as is

  dbound: LDL' diagonal threshold: 0

    Entries with abs. value less than dbound are replaced
with +/- dbound.

  grow0: memory reallocation: 1.2

  grow1: memory reallocation: 1.2

  grow2: memory reallocation: 5

  nrelax, zrelax: supernodal amalgamation rule:

    s = # columns in two adjacent supernodes

    z = % of zeros in new supernode if they are merged.

    Two supernodes are merged if (s <= 4) or (no new zero
entries) or

    (s <= 16 and z < 80%) or (s <= 48 and z < 10%) or (z < 5%)

  OK

sp\: But A is real symmetric; try MA57.

MA57 Control Parameters

    Ordering used: AMD ordering using MC47.

    Block Size for Level 3 BLAS: 16

    Merge assembly tree nodes if number of eliminations

    is less than: 16

    Matrix is scaled using symmetrized version of HSL code MC64.

    Maximum iterative refinement steps: 10

    Pivot Threshold: 1.000000e-02

    Zero Pivot Tolerance: 1.000000e-20

sp\: is MA57's symbolic analysis successful? yes.

    Forecast number of Real entries in factors: 191653054

    Forecast number of Integer entries in factors: 1466442

    Forecast maximum frontal size: 8071

    Number of nodes in Assembly tree: 5001

    Forecast length of FACT array (real)

    without numerical pivot: 288897514

    Forecast length of ifact array (integer)

    without numerical pivot 15335967

    Length of fact required for a successful

    completion of the numerical phase allowing

    data compression (without numerical pivoting) 284102031

    Length of ifact required for a successful

    completion of the numerical phase allowing

    data compression (without numerical pivoting) 15335967

    Number of data compresses 0

    Forecast number of floating point additions

    for the assembly processes 6.036817e+08

    Forecast number of floating point operations

    to perform the elimination operations

    counting multiply-add pair as two operations 7.745417e+11

    *** Return code from MA57AD: 1

??? Error using ==> mldivide

Out of memory. Type HELP MEMORY for your options.

 

Can anyone suggest me a solution or a different strategy to
solve the original problem Ax=b?


Thank you
 

"Alessio Rucci" <alessio.rucci@libero.it> wrote in message
<fguq3f$ig1$1@fred.mathworks.com>...
> I have to solve Ax=b where A is a large sparse matrix (1e+6
> x 1e+5) with about
>
> 7e+6 non-zero elements.
>
>
>
> I tried to solve A'Ax=A'b using pcg, but the new matrix
> (A'A) is very bad conditioned. I tried to use cholinc to
> precondition the system, but A’A has about 1e+7 non-zero
> elements, and Matlab gives me out of memory error.
>
>
>
> I tried also to solve the system using:
>
> x=(A’A)\(A’b)
>
>
>
> but here it is what I get:
>
>
>
> sp\: bandwidth = 5864+1+5864.
>
> sp\: is A diagonal? no.
>
> sp\: is band density (0.01) > bandden (0.50) to try banded
> solver? no.
>
> sp\: is A triangular? no.
>
> sp\: is A morally triangular? no.
>
> sp\: is A a candidate for Cholesky (symmetric, real positive
> diagonal)? yes.
>
> sp\: is CHOLMOD's symbolic Cholesky factorization (with
> automatic reordering) successful? yes.
>
> sp\: is CHOLMOD's numeric Cholesky factorization
successful? no.
>
>
>
> CHOLMOD version 1.1.0, May 5, 2006: : status: ERROR, out of
> memory
>
> Architecture: Linux
>
> sizeof(int): 4
>
> sizeof(UF_long): 4
>
> sizeof(void *): 4
>
> sizeof(double): 8
>
> sizeof(Int): 4 (CHOLMOD's basic integer)
>
> sizeof(BLAS_INT): 4 (integer used in the BLAS)
>
> Results from most recent analysis:
>
> Cholesky flop count: 8.0627e+11
>
> Nonzeros in L: 1.9098e+08
>
> memory blocks in use: 11
>
> memory in use (MB): 8.8
>
> peak memory usage (MB): 1961.0
>
> maxrank: update/downdate rank: 8
>
> supernodal control: 1 40 (supernodal if flops/lnz >= 40)
>
> nmethods=0: default strategy: Try user permutation if
> given. Try AMD.
>
> Select best ordering tried.
>
> method 0: user permutation (if given)
>
> method 1: AMD (or COLAMD if factorizing AA')
>
> prune_dense: for pruning dense nodes: 10
>
> a dense node has degree >= max(16,(10)*sqrt(n))
>
> flop count: 8.0627e+11
>
> nnz(L): 1.9098e+08
>
> final_asis: TRUE, leave as is
>
> dbound: LDL' diagonal threshold: 0
>
> Entries with abs. value less than dbound are replaced
> with +/- dbound.
>
> grow0: memory reallocation: 1.2
>
> grow1: memory reallocation: 1.2
>
> grow2: memory reallocation: 5
>
> nrelax, zrelax: supernodal amalgamation rule:
>
> s = # columns in two adjacent supernodes
>
> z = % of zeros in new supernode if they are merged.
>
> Two supernodes are merged if (s <= 4) or (no new zero
> entries) or
>
> (s <= 16 and z < 80%) or (s <= 48 and z < 10%) or (z < 5%)
>
> OK
>
> sp\: But A is real symmetric; try MA57.
>
> MA57 Control Parameters
>
> Ordering used: AMD ordering using MC47.
>
> Block Size for Level 3 BLAS: 16
>
> Merge assembly tree nodes if number of eliminations
>
> is less than: 16
>
> Matrix is scaled using symmetrized version of HSL code
MC64.
>
> Maximum iterative refinement steps: 10
>
> Pivot Threshold:
1.000000e-02
>
> Zero Pivot Tolerance:
1.000000e-20
>
> sp\: is MA57's symbolic analysis successful? yes.
>
> Forecast number of Real entries in factors: 191653054
>
> Forecast number of Integer entries in factors: 1466442
>
> Forecast maximum frontal size: 8071
>
> Number of nodes in Assembly tree: 5001
>
> Forecast length of FACT array (real)
>
> without numerical pivot: 288897514
>
> Forecast length of ifact array (integer)
>
> without numerical pivot 15335967
>
> Length of fact required for a successful
>
> completion of the numerical phase allowing
>
> data compression (without numerical pivoting) 284102031
>
> Length of ifact required for a successful
>
> completion of the numerical phase allowing
>
> data compression (without numerical pivoting) 15335967
>
> Number of data compresses 0
>
> Forecast number of floating point additions
>
> for the assembly processes
6.036817e+08
>
> Forecast number of floating point operations
>
> to perform the elimination operations
>
> counting multiply-add pair as two operations
7.745417e+11
>
> *** Return code from MA57AD: 1
>
> ??? Error using ==> mldivide
>
> Out of memory. Type HELP MEMORY for your options.
>
>
>
> Can anyone suggest me a solution or a different strategy to
> solve the original problem Ax=b?
>
>
> Thank you
>

My last post seems to have gotten lost. Or not.
I clicked "submit" but I got sent to the "update your
profile" page. So this might be a duplicate.

Upload the matrix to http://www.cise.ufl.edu/~web-gfs
and let me take a crack at it. I'd like to add it to
my collection (http://www.cise.ufl.edu/research/sparse).
Send A itself, not A'*A, as a mat file.

x=(A'*A)\(A'*b) is not a good idea, since it is very
ill conditioned. The failure here is CHOLMOD running
out of memory, however. It's puzzling that backslash
then tries MA57. CHOLMOD and MA57 both use AMD, and
the latter tends to use a tad bit more memory. So if
CHOLMOD+AMD runs out of memory, MA57+AMD almost certainly
will.

You should use sparse QR, although that doesn't use
the BLAS and is thus slower than sparse Cholesky.

Try using METIS. There is an interface to the METIS
ordering method in SuiteSparse on the File Exchange.
METIS can give better orderings than AMD, although it
can segfault sometimes :-(

"
> >
>
> My last post seems to have gotten lost. Or not.
> I clicked "submit" but I got sent to the "update your
> profile" page. So this might be a duplicate.
>
> Upload the matrix to http://www.cise.ufl.edu/~web-gfs
> and let me take a crack at it. I'd like to add it to
> my collection (http://www.cise.ufl.edu/research/sparse).
> Send A itself, not A'*A, as a mat file.
>
> x=(A'*A)\(A'*b) is not a good idea, since it is very
> ill conditioned. The failure here is CHOLMOD running
> out of memory, however. It's puzzling that backslash
> then tries MA57. CHOLMOD and MA57 both use AMD, and
> the latter tends to use a tad bit more memory. So if
> CHOLMOD+AMD runs out of memory, MA57+AMD almost certainly
> will.
>
> You should use sparse QR, although that doesn't use
> the BLAS and is thus slower than sparse Cholesky.
>
> Try using METIS. There is an interface to the METIS
> ordering method in SuiteSparse on the File Exchange.
> METIS can give better orderings than AMD, although it
> can segfault sometimes :-(
>
>
>
>
>
>

Ok,I'll send you the matrix

Hi;

> > I have to solve Ax=b where A is a large sparse matrix
(1e+6 x 1e+5) with about 7e+6 non-zero elements.

> > I tried also to solve the system using:
> >
> > x=(A’A)\(A’b)
> >

Why don't you solve for x by simply putting

x = A \ b

I tried it using a sparse matrix and it works fine. If b is
a row matrix and A is transposed you can also try

x = b / A

There are 2 functions that you can use which do exactly the
same thing

x = mldivide(A,b); %Equivalent to A \ b
x = mrdivide(b,A); %Equivalent to b / A

Why don't you give it a try.

Sterren

>
> Upload the matrix to http://www.cise.ufl.edu/~web-gfs
> and let me take a crack at it. I'd like to add it to
> my collection (http://www.cise.ufl.edu/research/sparse).
> Send A itself, not A'*A, as a mat file.
>


Can I send you the matrix by a jumbo mail (your mail is that
in your profile?)?It makes my some problem to upload the
matrix on your website...

> Can I send you the matrix by a jumbo mail (your mail is that
> in your profile?)?It makes my some problem to upload the
> matrix on your website...

Honestly, it's not difficult at all. Let's say you have a
set of simultaneous equations

3a+b+7c = 7
11a+4b+6c = 3
9a+2b = -20

and you want to solve for a,b,c

Get this in the form Ax = b:

>>A = [3 1 7; 11 4 6; 9 2 0]

A =

     3 1 7
    11 4 6
     9 2 0

>> b = [7; 3; -20]

b =

     7
     3
   -20

>> x = A \ b

x =

   -4.9750
   12.3875
    1.3625

a = -4.9750, b = 12.3875, c = 1.3625

This can work with sparse matrices as well:

>>A = sparse(A);

A =

   (1,1) 3
   (2,1) 11
   (3,1) 9
   (1,2) 1
   (2,2) 4
   (3,2) 2
   (1,3) 7
   (2,3) 6

>> x = A \ b

x =

   -4.9750
   12.3875
    1.3625

I don't understand why you have problems.

Sterren

> I don't understand why you have problems.
>
> Sterren

The previus message was for Tim Devis that asked to upload
my matrix to his web site.

The problem is that my matrix is 1e+6 x 1e+5 ill conditioned....

"Sterren Latsky" <stez17@hotmail.com> wrote in message

> Honestly, it's not difficult at all. Let's say you have a
> set of simultaneous equations
...
>
> >> x = A \ b
>

That's what I meant by "use sparse QR", since if A is sparse
and rectangular, A\b uses sparse QR. The problem is
manifold: (1) A\b for this case uses the colmmd preordering,
which for this matrix gives an R with 261e6 nonzeros (R
would require about 3GB of memory), (2) sparse QR in MATLAB
is slow (chol is much faster), and (3) the matrix is
ill-conditioned so chol(A'*A) would give a bad answer, (4)
the matrix is huge. x=A\b just doesn't cut it here, for
those reasons.

With METIS(A,'col'), R has 115 million nonzeros.
METIS(A'*A) gives R with 128 million
COLAMD(A) gives R with 181 million
AMD(A'*A) gives R with 191 million
COLMMD(A) gives R with 261 million

so it would seem that metis(A,'col') gives the best column
ordering. You could try

p = metis(A,'col') ;
C = A(:,p) ;
x = (C'*C)\(C'*b) ;
x (p) = x ;

to use a direct method. This will be faster and use half
the memory as x=A\b. However, it might be inaccurate since
A is ill-conditioned (I'm trying it now; details to follow).

I will post the matrix in my sparse matrix collection,
shortly. I would need b as well, to try an iterative method.

I was able to solve it with a direct method, on my Pentium 4
which has 4GB of memory and a 3.2GHz clock cycle.

[x stats] = cholmod2 (A'*A, A'*b) ;
r = norm (A'*(A*x) - A'*b) ;

where b=rand(size(A,1),1) ;

CHOLMOD in this case is using the best of AMD or METIS. You
need to install CHOLMOD from the SuiteSparse on the File
Exchange to do this; METIS is not in MATLAB itself.

CHOLMOD got an rcond of 4e-5, so this is not too
ill-conditioned to use the normal equations (CHOLMOD's
condition estimate is 1/4e-5 or about 2e4). nnz(L) was 115
million. The flop count was 283 billion. Time was 94.7
seconds (CHOLMOD got 3 Gflops, or just less than 50% of the
theoretical peak. Total memory usage was 1283 MB, or about
1.3 GB. So if you have a machine, with say 2GB of main
memory, that should be sufficient.

Your results may vary ... please send me your b vector.
This problem is incomplete without it, and iterative solvers
cannot be tested without b.

Using AMD within CHOLMOD isn't sufficient because the fillin
is worse and the memory usage is too high.

This matrix is now posted at:

http://www.cise.ufl.edu/research/sparse/matrices/Rucci/Rucci1.html