Re: how to get a rational basis for the null space

[Thomas Shores]

There may well be an existing octave function somewhere that does the 
trick.  If so, ignore the following, which is probably far too long 
winded, but nonetheless entertaining (to me), so here it is:   Suppose 
you have a matrix whose columns you believe to be multiples of rational 
vectors. If they are unit vectors, as in the returned vectors from the 
svd, then each vector is a normalized vector, hence likely involves an 
irrational square root  factor.  How do you recover rational multiples 
of each column?  Let's do an example:

octave:1> a = [1 2 7; 1 2 7; 1 2 7]
a =

   1   2   7
   1   2   7
   1   2   7

octave:2> v = null(a)
v =

   0.653213   0.744845
   0.698262  -0.662085
  -0.292820   0.082761

Now we know that *some* multiple of each column is rational.  So if we 
multiply by a scalar and a single rational entry results, then all 
entries must be rational. Hence, we know that the following matrix w has 
rational entries, to machine precision:

octave:3> w = v./repmat(norm(v,inf,'cols'),size(v,1),1)
w =

   0.93548   1.00000
   1.00000  -0.88889
  -0.41935   0.11111

Sometimes it's obvious at this point what the fractions are.  Do you see 
them? If not, we need some help.  What we want are the convergents of a 
continued fraction approximation.  I've wanted these from time to time, 
so I wrote the function 'printrat.m' that I'm attaching to this message 
with a little bit of documentation added in.  Put it to work to find the 
numerators (num) and denominators (den) of the fractional forms of 
entries of v:

octave:4> [num,den]=printrat(w)

29/31  1/1
1/1     -8/9
-13/31  1/9

num =

   29    1
    1   -8
  -13    1

den =

   31    1
     1    9
   31    9


Now suppose we only wanted integer entries.  In this case inspection 
tells us what to multiply each vector by, but if we want to automate it 
somewhat, do this:

octave:5> x=(num./den)*diag([lcm(den(:,1)),lcm(den(:,2))])
x =

   29    9
   31   -8
  -13    1

(If you wanted a general one-liner without loops and are satisfied with 
non-optimal numbers, use "x=(num./den)*diag(prod(den,1))".)  Now check 
that these vectors are orthogonal null vectors for a:

octave:6> a*x
ans =

   0   0
   0   0
   0   0

octave:7> x'*x
ans =

   1971      0
      0    146

Yup, it checks.



On 06/14/2010 12:16 PM, John W. Eaton wrote:
> Sorry, this message wasn't supposed to appear to come from me.  It was
> a botched forward to the list.  I will try again and see if I can get
> the From: address to be correct...
>
> On 12-Jun-2010, John W. Eaton wrote:
>
> | Now I use octave to get the orthonormal basis for the null space of A,but I
> | found I can't get a "rational" basis for the null space.The function null
> | (A,'r') in matlab can do this,I wonder if there is a function like null(A,'r')
> | in octave. Thank you!
> |
> |
> | ----------------------------------------------------------------------
> | _______________________________________________
> | Help-octave mailing list
> |address@hidden
> |https://www-old.cae.wisc.edu/mailman/listinfo/help-octave
> _______________________________________________
> Help-octave mailing list
> address@hidden
> https://www-old.cae.wisc.edu/mailman/listinfo/help-octave
>

printrat.m
Description: application/vnd.wolfram.mathematica.package