[Octave-bug-tracker] [bug #53560] pinv BIST test failing occasionally

[Rik]

URL:
  <http://savannah.gnu.org/bugs/?53560>

                 Summary: pinv BIST test failing occasionally
                 Project: GNU Octave
            Submitted by: rik5
            Submitted on: Wed 04 Apr 2018 12:01:46 PM PDT
                Category: Test Suite
                Severity: 3 - Normal
                Priority: 5 - Normal
              Item Group: Unexpected Error
                  Status: None
             Assigned to: None
         Originator Name: 
        Originator Email: 
             Open/Closed: Open
         Discussion Lock: Any
                 Release: dev
        Operating System: GNU/Linux

    _______________________________________________________

Details:

The BIST test in pinv.cc is failing occasionally.  I noticed this first with
the buildbots
(http://buildbot.octave.org:8010/builders/gcc-fedora/builds/1448/steps/test/logs/stdio),
but I can reproduce it with the following code:


for i = 1:1e3
> bm(i) = test ('pinv.cc');
> endfor
sum (bm)
ans =  989


The BIST test is


%!shared a, b, tol, hitol, d, u, x, y
%! a = reshape (rand*[1:16], 4, 4);  # Rank 2 matrix
%! b = pinv (a);
%! tol = 4e-14;
%! hitol = 40*sqrt (eps);
%! d = diag ([rand, rand, hitol, hitol]);
%! u = rand (4);                     # Could be singular by freak accident
%! x = inv (u)*d*u;
%! y = pinv (x, sqrt (eps));
%!
%!assert (a*b*a, a, tol)
%!assert (b*a*b, b, tol)
%!assert ((b*a)', b*a, tol)
%!assert ((a*b)', a*b, tol)
%!assert (x*y*x, x, -hitol)
%!assert (y*x*y, y, -hitol)
%!assert ((x*y)', x*y, hitol)
%!assert ((y*x)', y*x, hitol)


Just for fun, I fixed the random seed at the top of the test and then I can
get 1000 successful tests.

Should we just fix the random seed?

I used the following code in a script tst_pinv.m


N = 1e3;

tol = 4e-14;
for i = 1:N
  a = reshape (rand*[1:16], 4, 4);  # Rank 2 matrix
  b = pinv (a);
  assert (a*b*a, a, tol);
  assert (b*a*b, b, tol);
  assert ((b*a)', b*a, tol);
  assert ((a*b)', a*b, tol);
endfor


When it errors out I can look at the results.  For the failing case the
condition number of a is greater than 2e17.  Should we really be using such a
difficult matrix for testing?  Even without the randomization, the condition
number is quite high.


cond (reshape(1:16, 4, 4))
ans =    9.3757e+16


This is because the columns are not linearly independent.  You can add the
value 4 to a column to turn it in to the adjacent column.

As a reference, the condition number of a random 4x4 matrix is much smaller.


cond (rand (4))
ans =  42.840
cond (rand (4))
ans =  11.335


Maybe we should be doing something like this instead?


a = rand (4, 4);
while (cond (a) > 100)
  a = rand (4,4);
endwhile


Hmmphh.  Just tried that and it didn't work either.





    _______________________________________________________

Reply to this item at:

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

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