Re: diagnosing short rank matrices

[John W. Eaton]

On  1-Aug-2003, Mike Miller <address@hidden> wrote:

| There's probably a more straightforward
| way (maybe using eigenvalues and eigenvectors?)

Yes, you can use the null function.  Octave's implementation is based
on svd.

| You probably want to do this sort of thing:
| 
| max(abs(X(:,[[1:(k-1)],[(k+1):m]])*B-X(:,k))) < tiny_number
| 
| because it won't be *exactly* zero.

Right, you'll also need to do the same with the result returned from
null.  It will return an empty matrix for full rank matrices.  For
others, you will need to throw out elements that are nearly zero.  For
example, you could do something like this:

  octave:8> a = [1,2,3;2,4,6;4,2,17]'
  a =

     1   2   4
     2   4   2
     3   6  17

  octave:9> x = null (a)
  x =

      8.9443e-01
     -4.4721e-01
      4.3775e-18

  octave:10> x(abs (x) < eps) = 0
  x =

     0.89443
    -0.44721
     0.00000

  octave:11> a*x
  ans =

     -1.1102e-16
     -2.2204e-16
     -3.3307e-16


This tells you that the first two columns are linearly dependent, and
gives you the (normalized) coefficients for those columns that will
give you a (numerically) zero vector.

Note that there are some bugs in Octave's handling of empty matrices
that prevent the above code from working correctly if null returns an
Nx0 empty matrix (which it will, for an NxN matrix full rank).  The
first bug is that

  zeros (n, 0) < eps

returns [] instead of [](0x1).  This may simply be a Matlab
compatibility problem.  I'm not really sure why it should return an
empty column vector instead of a 0x0 empty matrix.  Perhaps because it
is not a row vector?  But then it also has no columns, so...

The second is that even if it returned the Matlab-compatible result,

  x = zeros (n, 0);
  x (zeros (n, 0) < eps)

returns [](0x1) instead of [](nx0), so the final a*x will still fail.

If not for these problems, you could use the above code without special
checks for empty matrices.  I will see about fixing them.

Thanks,

jwe



-------------------------------------------------------------
Octave is freely available under the terms of the GNU GPL.

Octave's home on the web:  http://www.octave.org
How to fund new projects:  http://www.octave.org/funding.html
Subscription information:  http://www.octave.org/archive.html
-------------------------------------------------------------