Got Questions? Get Answers.
Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
Don't detect singular sparse matrix

Subject: Don't detect singular sparse matrix

From: Sebastian

Date: 1 May, 2010 16:58:05

Message: 1 of 9

Why in the some cases the linear system solutions using sparse matrix (with mrdivide or inv), don't is indicated that the matrix is closed to singular?

By example in version 2010a:
a = speye(10,10);
a(1,5) = -1;
a(5,1) = -1;
b = randn(10,1);
sprank(a);
>> ans = 10
rank(full(a));
>> ans = 9
a\b;
(no warning)
inv(a)*b;
(no warning)
full(a)\b
Warning: Matrix is singular to working precision.

But in version 2007a:
a = speye(10,10);
a(1,5) = -1;
a(5,1) = -1;
b = randn(10,1);
sprank(a);
>> ans = 10
rank(full(a));
>> ans = 9
a\b;
Warning: Matrix is singular to working precision.
inv(a)*b;
Warning: Matrix is singular to working precision.
full(a)\b
Warning: Matrix is singular to working precision.

Do exist a way always is indicated that matrix is closed to singular?

Thanks.

Subject: Don't detect singular sparse matrix

From: John D'Errico

Date: 1 May, 2010 17:38:08

Message: 2 of 9

"Sebastian " <quit.spaseba@email.google> wrote in message <hrhmit$jl9$1@fred.mathworks.com>...
> Why in the some cases the linear system solutions using sparse matrix (with mrdivide or inv), don't is indicated that the matrix is closed to singular?
>
> By example in version 2010a:
> a = speye(10,10);
> a(1,5) = -1;
> a(5,1) = -1;
> b = randn(10,1);
> sprank(a);
> >> ans = 10
> rank(full(a));
> >> ans = 9
> a\b;
> (no warning)
> inv(a)*b;
> (no warning)
> full(a)\b
> Warning: Matrix is singular to working precision.
>
> But in version 2007a:
> a = speye(10,10);
> a(1,5) = -1;
> a(5,1) = -1;
> b = randn(10,1);
> sprank(a);
> >> ans = 10
> rank(full(a));
> >> ans = 9
> a\b;
> Warning: Matrix is singular to working precision.
> inv(a)*b;
> Warning: Matrix is singular to working precision.
> full(a)\b
> Warning: Matrix is singular to working precision.
>
> Do exist a way always is indicated that matrix is closed to singular?
>
> Thanks.

help cond
help condest
help rank

John

Subject: Don't detect singular sparse matrix

From: Bruno Luong

Date: 1 May, 2010 18:06:05

Message: 3 of 9

"John D'Errico" <woodchips@rochester.rr.com> wrote in message <hrhou0$it1$1@fred.mathworks.com>...
>
> help cond
> help condest
> help rank
>

I might add only CONDEST can really work on sparse matrix.

Personally I have a preference of using NORMEST/EIGS to estimate the condition number of sparse matrix.

The difference between Matlab version could be explained by the sparse QR factorization method introduced from recent Matlab versions.

Bruno

Subject: Don't detect singular sparse matrix

From: Sebastian

Date: 1 May, 2010 18:28:04

Message: 4 of 9

"John D'Errico" <woodchips@rochester.rr.com> wrote in message <hrhou0$it1$1@fred.mathworks.com>...
> "Sebastian " <quit.spaseba@email.google> wrote in message <hrhmit$jl9$1@fred.mathworks.com>...
> > Why in the some cases the linear system solutions using sparse matrix (with mrdivide or inv), don't is indicated that the matrix is closed to singular?
> >
> > By example in version 2010a:
> > a = speye(10,10);
> > a(1,5) = -1;
> > a(5,1) = -1;
> > b = randn(10,1);
> > sprank(a);
> > >> ans = 10
> > rank(full(a));
> > >> ans = 9
> > a\b;
> > (no warning)
> > inv(a)*b;
> > (no warning)
> > full(a)\b
> > Warning: Matrix is singular to working precision.
> >
> > But in version 2007a:
> > a = speye(10,10);
> > a(1,5) = -1;
> > a(5,1) = -1;
> > b = randn(10,1);
> > sprank(a);
> > >> ans = 10
> > rank(full(a));
> > >> ans = 9
> > a\b;
> > Warning: Matrix is singular to working precision.
> > inv(a)*b;
> > Warning: Matrix is singular to working precision.
> > full(a)\b
> > Warning: Matrix is singular to working precision.
> >
> > Do exist a way always is indicated that matrix is closed to singular?
> >
> > Thanks.
>
> help cond
> help condest
> help rank
>
> John

I known that they are estimators of the singularities (is it that I have infer of your contestation?), but why it change between versions? In my program it behavior is dangerous because don't is indicated that the linear system is wrong definided. Besides that return a valid result (No NaNs) and I don't know its interpretation.

Thanks

Subject: Don't detect singular sparse matrix

From: Bruno Luong

Date: 1 May, 2010 18:54:05

Message: 5 of 9

"Sebastian " <quit.spaseba@email.google> wrote in message <hrhrrk$pe1$1@fred.mathworks.com>...
>
>
> I known that they are estimators of the singularities (is it that I have infer of your contestation?), but why it change between versions?

As I wrote earlier, there is a new linear algebra algorithm implemented in recent version, which make the sparse backslash more robust.

For example try
  [q r e] = qr(a)

in both versions.

>In my program it behavior is dangerous because don't is indicated that the linear system is wrong definided. Besides that return a valid result (No NaNs) and I don't know its interpretation.

I always consider it is user's duty to make sure his matrix is not singular, and do some workaround otherwise. Matlab gives a warning, fine, but I wouldn't trust on the warning issued in order to catch the singularity of my problem.

Bruno

Subject: Don't detect singular sparse matrix

From: John D'Errico

Date: 1 May, 2010 20:26:06

Message: 6 of 9

"Sebastian " <quit.spaseba@email.google> wrote in message <hrhrrk$pe1$1@fred.mathworks.com>...
> "John D'Errico" <woodchips@rochester.rr.com> wrote in message <hrhou0$it1$1@fred.mathworks.com>...
> > "Sebastian " <quit.spaseba@email.google> wrote in message <hrhmit$jl9$1@fred.mathworks.com>...
> > > Why in the some cases the linear system solutions using sparse matrix (with mrdivide or inv), don't is indicated that the matrix is closed to singular?
> > >
> > > By example in version 2010a:
> > > a = speye(10,10);
> > > a(1,5) = -1;
> > > a(5,1) = -1;
> > > b = randn(10,1);
> > > sprank(a);
> > > >> ans = 10
> > > rank(full(a));
> > > >> ans = 9
> > > a\b;
> > > (no warning)
> > > inv(a)*b;
> > > (no warning)
> > > full(a)\b
> > > Warning: Matrix is singular to working precision.
> > >
> > > But in version 2007a:
> > > a = speye(10,10);
> > > a(1,5) = -1;
> > > a(5,1) = -1;
> > > b = randn(10,1);
> > > sprank(a);
> > > >> ans = 10
> > > rank(full(a));
> > > >> ans = 9
> > > a\b;
> > > Warning: Matrix is singular to working precision.
> > > inv(a)*b;
> > > Warning: Matrix is singular to working precision.
> > > full(a)\b
> > > Warning: Matrix is singular to working precision.
> > >
> > > Do exist a way always is indicated that matrix is closed to singular?
> > >
> > > Thanks.
> >
> > help cond
> > help condest
> > help rank
> >
> > John
>
> I known that they are estimators of the singularities (is it that I have infer of your contestation?), but why it change between versions? In my program it behavior is dangerous because don't is indicated that the linear system is wrong definided. Besides that return a valid result (No NaNs) and I don't know its interpretation.
>
> Thanks

As Bruno says, don't simply put blind trust in the
warning messages. Know your problem, or at least
learn how to determine the characteristics of your
problem.

John

Subject: Don't detect singular sparse matrix

From: Brian Borchers

Date: 1 May, 2010 21:34:51

Message: 7 of 9

I'm comfortable with \ returning some kind of reasonable least squares
solution, and I'm also comfortable with the way that MATLAB deals with
inv() acting on the dense matrix in this example:

>> A=eye(5); A(1,5)=-1; A(5,1)=-1;
>> inv(A)
Warning: Matrix is singular to working precision.

ans =

   Inf Inf Inf Inf Inf
   Inf Inf Inf Inf Inf
   Inf Inf Inf Inf Inf
   Inf Inf Inf Inf Inf
   Inf Inf Inf Inf Inf


However, I can't see how the behavior of inv() acting on the sparse
matrix can be considered anything other than a bug:

>> AS=sparse(A);
>> inv(AS)

ans =

   (1,1) 2
   (5,1) 2
   (2,2) 1
   (3,3) 1
   (4,4) 1
   (1,5) 2
   (5,5) 3

>> AS*inv(AS)

ans =

   (2,2) 1
   (3,3) 1
   (4,4) 1
   (1,5) -1
   (5,5) 1

This seems to be completely inconsistent with the documentation in
help inv

>> help inv
 INV Matrix inverse.
    INV(X) is the inverse of the square matrix X.
    A warning message is printed if X is badly scaled or
    nearly singular.

    See also SLASH, PINV, COND, CONDEST, LSQNONNEG, LSCOV.

Subject: Don't detect singular sparse matrix

From: Bruno Luong

Date: 2 May, 2010 06:40:05

Message: 8 of 9

Brian Borchers <borchers.brian@gmail.com> wrote in message <29f7d81e-3a71-41ab-93df-6cc4d1622ea8@v12g2000prb.googlegroups.com>...
>
> However, I can't see how the behavior of inv() acting on the sparse
> matrix can be considered anything other than a bug:
>
> >> AS=sparse(A);
> >> inv(AS)
>
> ans =
>
> (1,1) 2
> (5,1) 2
> (2,2) 1
> (3,3) 1
> (4,4) 1
> (1,5) 2
> (5,5) 3
>
> >> AS*inv(AS)
>
> ans =
>
> (2,2) 1
> (3,3) 1
> (4,4) 1
> (1,5) -1
> (5,5) 1
>

Yes, arguably it's a bug.

Bruno

Subject: Don't detect singular sparse matrix

From: Tim Davis

Date: 27 May, 2010 23:59:07

Message: 9 of 9

The matrix is square, so it's not sparse QR. You can ask MATLAB what it is using via:

spparms('spumoni',2)


Try this:


>> A=eye(5); A(1,5)=-1; A(5,1)=-1;
>> A

A =

     1 0 0 0 -1
     0 1 0 0 0
     0 0 1 0 0
     0 0 0 1 0
    -1 0 0 0 1

>> inv(A)
Warning: Matrix is singular to working precision.

ans =

   Inf Inf Inf Inf Inf
   Inf Inf Inf Inf Inf
   Inf Inf Inf Inf Inf
   Inf Inf Inf Inf Inf
   Inf Inf Inf Inf Inf

>> inv(sparse(A))

ans =

   (1,1) 2
   (5,1) 2
   (2,2) 1
   (3,3) 1
   (4,4) 1
   (1,5) 2
   (5,5) 3

>> spparms('spumoni',2)
>> inv(sparse(A))
sp\: bandwidth = 4+1+4.
sp\: is A diagonal? no.
sp\: is band density (0.28) > 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.
sp\: A is not positive definite.

CHOLMOD version 1.7.0, Sept 20, 2008: : status: warning, matrix not positive definite
  Architecture: unknown
    sizeof(int): 4
    sizeof(UF_long): 8
    sizeof(void *): 8
    sizeof(double): 8
    sizeof(Int): 8 (CHOLMOD's basic integer)
    sizeof(BLAS_INT): 8 (integer used in the BLAS)
  Results from most recent analysis:
    Cholesky flop count: 8
    Nonzeros in L: 6
  memory blocks in use: 13
  memory in use (MB): 0.0
  peak memory usage (MB): 0.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
        nnz(L): 6
  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: 6
    Forecast number of Integer entries in factors: 16
    Forecast maximum frontal size: 2
    Number of nodes in Assembly tree: 4
    Forecast length of FACT array (real)
    without numerical pivot: 39
    Forecast length of ifact array (integer)
    without numerical pivot 44
    Length of fact required for a successful
    completion of the numerical phase allowing
    data compression (without numerical pivoting) 39
    Length of ifact required for a successful
    completion of the numerical phase allowing
    data compression (without numerical pivoting) 44
    Number of data compresses 0
    Forecast number of floating point additions
    for the assembly processes 6.000000e+00
    Forecast number of floating point operations
    to perform the elimination operations
    counting multiply-add pair as two operations 8.000000e+00
    *** Return code from MA57AD: 1
sp\: is MA57's numeric Factorization successful? yes.
    Number of entries in factors: 6
    Storage for real data in factors: 6
    Storage for integer data in factors: 19
    Minimum length of fact required for a successful
    completion of the numerical phase: 39
    Minimum length of ifact required for a successful
    completion of the numerical phase: 44
    Minimum length of fact required for a successful
    completion of the numerical phase without
    data compression: 39
    Minimum length of ifact required for a successful
    completion of the numerical phase without
    data compression: 44
    Order of the largest frontal matrix: 2
    Number of 2x2 numerical pivots: 0
    Total number of fully-summed variables passed to
    the father node because of numerical pivot: 1
    Number of negative eigenvalues: 0
    Rank of factorization of the matrix: 4
    Pivot step where static pivot commences and
    is set to zero if no modification performed: 0
    Number of data compresses on real
    data structures: 0
    Number of data compresses on integer
    data structures: 0
    Number of block pivots in factors 5
    Number of zeros in triangle of factors: 0
    Number of zeros in rectangle of factors: 0
    Number of zero columns in rectangle of factors: 0
    Number of static pivots chosen: 0
    Number of floating point additions
    for the assembly processes: 6.000000e+00
    Number of floating point operations
    to perform the elimination operations
    counting multiply-add pair as two operations: 7.000000e+00
    Minimum value of the scaling factor (MC64): 1.000000e+00
    Maximum value of the scaling factor (MC64): 1.000000e+00
    Maximum modulus of matrix entry: 1.000000e+00
    *** Return code from MA57BD: 4
sp\: is MA57's solve successful? yes.

ans =

   (1,1) 2
   (5,1) 2
   (2,2) 1
   (3,3) 1
   (4,4) 1
   (1,5) 2
   (5,5) 3

>> full(ans)

ans =

     2 0 0 0 2
     0 1 0 0 0
     0 0 1 0 0
     0 0 0 1 0
     2 0 0 0 3

>>

From the MA57 spec sheet, at http://www.hsl.rl.ac.uk/specs/ma57.pdf this looks like a warning that the matrix is rank deficient:

+4 Matrix is rank deficient on exit from MA57B/BD. In this case, a decomposition will still have been produced that will enable the subsequent solution of consistent equations. INFO(25) will be set to the rank of the factorized matrix.

And sure enough, it is rank deficient:

>> rank(A)

ans =

     4

>> svd(A)

ans =

     2
     1
     1
     1
     0

(the rank(A) of 4 and the return code of 4 are coincidental). So the bug is not MA57, but the interface to it (MATLAB) that returns inv(A) in spite of the fact that the matrix is rank deficient.

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us