[Octave-bug-tracker] [bug #33503] null.m: tolerance for zero in basis elements

[Olaf Till]

Follow-up Comment #5, bug #33503 (project octave):



(Tried the nomarkup tag, hope the text will be readable ...)

Here is the demanded explanation.

For the assessment of numeric examples, consider an 3-by-m matrix
A=[A1;A2] that has rank 2, but the 2-by-m matrix A2 has only rank
1. If A.'*x==A1.'*x(1)+A2.'*x(2:3)==y, since rank(A2)<2, different linear
combinations of A2(:,1) and A2(:,2) could yield the same result, but since by
including the single row of A1 into A the rank is increased by 1, no linear
combination of A2(:,1) and A2(:,2) can be equal to non_null_constant*A1. So
x(1) is uniquely determined (and x(2:3) is not). So all basis vectors of
null(A.')  should have their first element equal to zero.

Since the numeric problems seem to get worse with the following, I would also
like to show that, if A*A.'*x==b, then x(1) is uniquely defined (and so the
first elements of all basis-vectors returned by null(A*A.') should be zero).
I'll base this on optimization theory:

In the linear model function y=A.'*x, for a given y, x(1) is uniquely defined
as shown above. Now let y=A.'*x be the (unique) one that minimizes sumsq(z-y)
for a given z. Then, x is also given by A*A.'x=A*z. This completes the prove
(I hope without mistakes ...).

As a numeric example for A, take

octave:1> A_correct = [log(1:10); 1:10; 2*(1:10)];

and

octave:2> A_wrong = A_correct([2, 3, 1], :);

A_correct has rank 2, but A_correct([2, 3],:) has only rank
1. A_wrong only has changed order of rows, so that A_wrong([1, 2],:) has rank
1; this changes the above argumentation so that x(3) is uniquely defined and
the third element of the basis vectors should be zero.

Now

octave:3> null (A_correct.')
ans =

   0.000000000000000
   0.894427190999916
  -0.447213595499958

and

octave:4> null (A_correct*A_correct.')
ans =

   0.000000000000000
  -0.894427190999916
   0.447213595499958

work correctly,

but

octave:5> null (A_wrong.')
ans =

   8.94427190999916e-01
  -4.47213595499958e-01
  -7.28583859910259e-16

and

octave:6> null (A_wrong*A_wrong.')
ans =

   8.94427190999916e-01
  -4.47213595499959e-01
   6.21724893790088e-15

don't.

If I use the scaling

octave:7> A = A_wrong * diag (sqrt (.7) * ones (1, 10));

I get even larger values of the third element:

octave:8> null (A.')
ans =

  -8.94427190999916e-01
   4.47213595499958e-01
  -3.57353036051222e-15

and

octave:9> null (A * A.')
ans =

   8.94427190999915e-01
  -4.47213595499960e-01
   1.48908663177849e-14

To summarize, the case null(A.') is clearly not treated as it
should. But although in null(A*A.') a part of the inaccuracy probably stems
from the multiplication, null() is likely to be used in such a way by some
codes. So I suggest the case null(A*A.') should also be treated so that it
yields the expected result.

Still I don't know how to find the correct tolerance, except that I suppose
that it should be a fixed value, since the basis vectors are normalized.

-n


    _______________________________________________________

Reply to this item at:

  <http://savannah.gnu.org/bugs/?33503>

_______________________________________________
  Message sent via/by Savannah
  http://savannah.gnu.org/