[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: diagnosing short rank matrices
From: 
John W. Eaton 
Subject: 
Re: diagnosing short rank matrices 
Date: 
Fri, 1 Aug 2003 20:20:34 0500 
On 1Aug2003, 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:(k1)],[(k+1):m]])*BX(:,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.9443e01
4.4721e01
4.3775e18
octave:10> x(abs (x) < eps) = 0
x =
0.89443
0.44721
0.00000
octave:11> a*x
ans =
1.1102e16
2.2204e16
3.3307e16
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 Matlabcompatible 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
