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

[Olaf Till]

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

The point was to show that the numerical error of "eps" assumed by zero() to
set basis elements to zero is too low, i.e. the real numerical error is
higher. To show this, among others I had to prove that the exact result would
be zero (so the deviation from zero is numerical error). I also gave an
example for such a deviation, as an argument why the tolerance of zero()
should be set to something else.

Letting the user handle numerical error is only an option if the error-bound
is indeed fixed; otherwise, the user would have to repeat the singular value
decomposition. Even if fixed, the user would have to know the error bound ---
if the programmer knows it, it should be already considered in the function
IMO.

It turned out that LAPACK itself documents the error bound:

http://www.netlib.org/lapack/lug/node96.html

, and that the error bound is indeed not fixed, but must be computed with the
result of the singular value decomposition. I attach a respective patch
(against current unstable). It defines an interface function to LAPACKs SDISNA
and DDISNA and uses this function in null() to compute the error bound for the
angle between the computed basis vectors and the exact basis vectors. From
this angle, elements settable to zero are derived.

There is still some inaccuracy in considering these bounds. The influence of
removing some small elements on angle is so small that it is just computed to
be zero. But one can conclude that the effectively used error bound is in the
magnitude of order of the square of eps (i.e. much larger than before (just
eps)).



(file #23521)
    _______________________________________________________

Additional Item Attachment:

File name: null-space-zero-elements.changeset Size:7 KB


    _______________________________________________________

Reply to this item at:

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

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