BUG: invert a square sparse matrix.

[David Bateman]

Julien MICHOT wrote:
> hi,
> i have a bug when i want to invert a square sparse matrix.
> It says that it's not square, but it is !
>
> size(d)
>
> ans =
>
> 1800 1800
>
>
>
> unzip the file from ( http://dl.free.fr/oJ0cWW7nw/erroroctave.zip )
> *load -binary erroroctave.bin;
>   
You don't need the "-binary" here. Octave can figure that out for itself.

> k=spinv(d);*
>
> Error :
> error: inverse requires square matrix
>
> error: operator *: nonconformant arguments (op1 is 1800x1800, op2 is 0x0)
>
> error: operator *: nonconformant arguments (op1 is 0x0, op2 is 1800x1800)
>
> error: evaluating assignment expression near line 8, column 2
>
> Octave 3.0.0 under Windows XP 32.
>
> Maybe the matrix is too big, or it cannot invert it.
> It appends generally with big matrix.
>
>   
Yes there is a small bug, due to the fact that this is a near singular
hermitian matrix.. How Octave should work in this case is that it should
try a cholesky factorization and then if that fails fall back to an LU
factorization. What is happening in your case is that the fallback to
the LU factorization is not taking place (I'll figure out why and send a
patch later). You can see the matrix is singular if you force an LU
factorization with

load erroroctave.bin
d = matrix_trpe (d, 'Full');
e = inv (d);

Please use "inv" and not spinv, as spinv has disappeared in the 3.1.x
tree and the overloading for sparse matrices occurs within the existing
inv function. For the above I see

warning: inverse: matrix singular to machine precision, rcond = 1.41022e-27

and a probably incorrect inverse is returned.. Note that the inverse of
a sparse matrix is rarely a good idea... I can't think of many cases you
need the inverse, and the mldivide operator is probably better, and a
matrix inverse for strongly connected sparse matrices like your is a
great way of turning it into a dense one. You can identify a strong
connected sparse matrix with

[p,q,r,s] = dmperm (d);

where if the matrix is strongly connected the size of r and s are small
relative to the order of d. Check

octave:30> nnz(d)
ans =  906300
octave:31> nnz(e)
ans =  3240000
octave:32> prod(size(d))
ans =  3240000

Regards
David




-- 
David Bateman                                David.Bateman at motorola.com
Motorola Labs - Paris                        +33 1 69 35 48 04 (Ph) 
Parc Les Algorithmes, Commune de St Aubin    +33 6 72 01 06 33 (Mob) 
91193 Gif-Sur-Yvette FRANCE                  +33 1 69 35 77 01 (Fax) 

The information contained in this communication has been classified as: 

[x] General Business Information 
[ ] Motorola Internal Use Only 
[ ] Motorola Confidential Proprietary