[Octave-bug-tracker] [bug #53242] tests: 'test subspace.m' fails randomly

[Dan Sebald]

Follow-up Comment #2, bug #53242 (project octave):

I'm attaching a short test script that will break when the assert fails so
that one can investigate.  As I suspected, the assert fails when the two
random vectors are nearly coincident, i.e., very small angle between them. 
This doesn't mean that the values of the vectors are small or anything else
special about them, just that the two vectors point in nearly the same
direction.  E.g.,


octave:42> a
a =

   0.356464193870818
   0.134784034792930

octave:43> b
b =

   0.1672006432659599
   0.0630492086110193

octave:44> a'*b
ans =  0.0680990692435769
octave:45> a(1)/a(2)
ans =  2.64470635872013
octave:46> b(1)/b(2)
ans =  2.65190708891369


and


octave:49> a
a =

   0.985524629933309
   0.915307704770172

octave:50> b
b =

   0.858594918710192
   0.805314467477197

octave:51> a'*b
ans =  1.58327697636925
octave:53> a(1)/a(2)
ans =  1.07671401081538
octave:54> b(1)/b(2)
ans =  1.06616105060164


The best characterization for the fail condition is that the follow expression
is always very near but less than 1.0:


  a1 = norm (a,2);
  b1 = norm (b,2);
  dot (a,b)/(a1*b1)


The subspace() routine isn't too complicated, so we can probably pick this
apart.


  A = orth (A);
  B = orth (B);
  c = A'*B;
  scos = min (svd (c));
  if (scos^2 > 1/2)
    if (columns (A) >= columns (B))
      c = B - A*c;
    else
      c = A - B*c';
    endif
    ssin = max (svd (c));
    ang = asin (min (ssin, 1));
  else
    ang = acos (scos);
  endif


I've tested up to the point of 


c=A'*B


and it is exactly accurate with the normalized dot product:


octave:80> dot (a,b)/(a1*b1) - c
ans =    2.22044604925031e-16
octave:81> eps
ans =    2.22044604925031e-16


There are some SVDs in the code, but these are like degenerate SVDs because of
the dimensionality of 1.  So I suspect they don't do anything to the value. 
But notice that the subspace() routine takes the difference of the vectors
followed by arcsine, and therein lies the issue, most likely.  To see that,
compare the graphs:


subplot(2,1,1);
plot([0:0.01:1.0], asin(0:0.01:1.0));
title('arcsine');
subplot(2,1,2);
plot([0:0.01:1.0], acos(0:0.01:1.0));
title('arccosine');


where the inverse-trig function to get the value approximately equal to 0
means operating in two different "realms", as it were, for asin() and acos().

I've tested some trig formula variations rather than direct acos() and the
same result occurs.  The idea was to try using asin(), just as the subspace
routine uses asin().  Ultimately, I think this test boils down to comparing
the trig relationship between acos() and asin().

Do you want to pursue further finding a formula that results in better
accuracy when comparing two vectors that are nearly coincident?  (Might end up
with basically the same formula outside subspace() as inside subspace().) 
Ensure that the random vectors selected are outside that "near-extremal-value"
zone?  Relax the tolerance?

If we are testing random vectors, we should probably test more than just once.
 For example, the following doesn't seem to fail:


for i=1:500
  a = rand (2,1);
  b = rand (2,1);
  a1 = norm (a,2);
  b1 = norm (b,2);
  arg = dot (a,b)/(a1*b1);
  theta = asin (sqrt(1 - (dot (a,b)/(a1*b1))^2));
  if (arg < 0.999 && abs (subspace (a, b) - theta) > 100*eps)
    abs (subspace (a, b) - theta)
    break;
  end
end


(file #43411)
    _______________________________________________________

Additional Item Attachment:

File name: test_script.m                  Size:0 KB


    _______________________________________________________

Reply to this item at:

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

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