bugGNU Octave - Bugs: bug #53242, tests: 'test subspace.m' fails...

 
 

bug #53242: tests: 'test subspace.m' fails randomly

Submitter:  Mike Miller <mtmiller>
Submitted:  Wed 28 Feb 2018 01:08:32 AM UTC
   
 
Category:  Test Suite Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Inaccurate Result
Status:  Fixed Assigned to:  None
Originator Name:  Open/Closed:  * Closed
Release:  * dev Operating System:  * GNU/Linux
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Tue 06 Mar 2018 04:20:18 PM UTC, comment #21: 

I applied a combined patch for the BIST tests here http://hg.savannah.gnu.org/hgweb/octave/rev/55021b03884d.

Closing report.

Rik <rik5>
Group administrator
Tue 06 Mar 2018 12:29:08 AM UTC, comment #20: 

Correct, merge or apply the one in comment #18 and follow up with the tests from Michael, whichever is more convenient.  Basically, all the previous tests should go.

Dan Sebald <sebald>
Tue 06 Mar 2018 12:19:28 AM UTC, comment #19: 

What needs to happen to close this bug report?  If I understand correctly, Michael's tests from comment #7 need to be merged with the patch from comment #18, and the whole thing given a commit message and checked in.

Rik <rik5>
Group administrator
Thu 01 Mar 2018 11:08:20 PM UTC, comment #18: 

Try that a third time, oy.

(file #43431)

Dan Sebald <sebald>
Thu 01 Mar 2018 11:04:01 PM UTC, comment #17: 

I forgot to change the comment


%! ## For random vectors


to


%! ## For small angle between subspaces


New patch attached.

(file #43430)

Dan Sebald <sebald>
Thu 01 Mar 2018 10:54:45 PM UTC, comment #16: 

Here's a new changeset with the test I suggested.  It removes both of the random tests, the second of which will be covered by M.L.'s more comprehensive examples.  Mike can add a new changeset, combine them into one, propose some other type of small-angle test, etc.

Instead of the 3*e length for the rotated vector, I set the length to rand(1) and saw no fails for 100000 tries.  I don't think it is worth randomizing anything here.

To illustrate that the test is checking what it's supposed to, here is the computation with and without the use of asin():


octave:2> theta = pi/200;
octave:3> Ry = [cos(theta) 0 sin(theta);0 1 0;-sin(theta) 0 cos(theta)];
octave:4> a = Ry*[3*e 0 0]';
octave:5> b = [1 1 0;1 -1 0]';
octave:6> theta - subspace (a, b)
ans =   -3.4694e-18
octave:7> eps
ans =    2.2204e-16
octave:8> A = a;
octave:9> B = b;
octave:10>   A = orth (A);
octave:11>   B = orth (B);
octave:12>   c = A'*B;
octave:13>   scos = min (svd (c));
octave:14>   ang = acos (scos);
octave:15> theta - ang
ans =    2.0511e-14


(file #43429)

Dan Sebald <sebald>
Thu 01 Mar 2018 08:08:47 PM UTC, comment #15: 

Dan - just noting that the provided patch does fix the problem with the current test approach, but is not yet ready to commit. Let me know if you have another changeset you want tested.

Mike Miller <mtmiller>
Group Member
Thu 01 Mar 2018 08:07:22 PM UTC, comment #14: 

I've filed bug #53254 for the separate issue that other tests are not preserving the random number generator state correctly.

Mike Miller <mtmiller>
Group Member
Thu 01 Mar 2018 07:58:39 PM UTC, comment #13: 

I could rewrite the test, but with the discussion with M.L. I think we are moving away from testing so many angles and instead boil it down to maybe a half-dozen to ten tests that exercise all the important aspects of a subspace comparison.  The randomness and so many tests in my original test script really don't add much because they are limited to 2D vector space, mostly have big angles, etc.  I just wanted to make sure all sectors of the space were tested.

Something like this


theta = pi/200;
Ry = [cos(theta) 0 sin(theta);0 1 0;-sin(theta) 0 cos(theta)];
assert(subspace(Ry*[1 0 0]',[1 1 0;1 -1 0]'), theta, eps);


is a small angle test that accomplishes the task.  Picking random vectors most times isn't small angle so doesn't test the core of the subspace algorithm.

We could make the hyperplane randomly oriented and then chose the line to be a small angle from that hyperplane, but I'm not sure that adds anything.  Maybe it would affect the orth() routine that is used by subspace(), but that's what "test orth" is for.

The one thing that intrigues me is why we are finding that the


assert(subspace([1 0 1]',[1 1 0;1 -1 0]'),pi/4,3*eps)


requires an increased tolerance of 3 times eps.  That's still not bad, but the funny thing is I implement the [1 0 1] as a rotation:


octave:2> theta = -pi/4;
octave:3> Ry = [cos(theta) 0 sin(theta);0 1 0;-sin(theta) 0 cos(theta)];
octave:4> Ry * [sqrt(2) 0 0]' - [1 0 1]'
ans =

   2.2204e-16
   0.0000e+00
   0.0000e+00

octave:5> assert(subspace(Ry*[sqrt(2) 0 0]',[1 1 0;1 -1 0]'),pi/4,eps)


and that little inaccuracy seems to compensate so the test passes using a tolerance of eps.  Just luck I guess.

Dan Sebald <sebald>
Thu 01 Mar 2018 07:28:44 PM UTC, comment #12: 

RE: comment #10 - If you look at several of the buildbot log files on that system, you will see that the subspace.m test is failing on the exact same value each time. Doesn't that suggest that maybe a prior test is resetting the random number generator to a particular fixed state (several tests do) and not resetting it to its original state? So I don't think the test is failing on that system because of the build configuration necessarily. But maybe exposing a logic bug in a separate unit test that is invalidating the randomness of all following tests.

For example, I see that scripts/sparse/eigs.m sets the random state but does not restore it. That should be fixed. But that can't be the culprit because it is tested after scripts/linear-algebra/subspace.m.

RE: comment #4 - Yes, Dan's test rewrite fixes the error for me, running locally on a fully built Octave with 1000 iterations of the test. It also drastically increases the test elapsed time from an average of about 4.5 ms to about 36.5 ms.

Dan - can you look at vectorizing the loops in your proposed test rewrite?

Mike Miller <mtmiller>
Group Member
Thu 01 Mar 2018 05:38:00 PM UTC, comment #11: 

Duh, not my day.  In comment #10, I meant to refer to comment #4.

John W. Eaton <jwe>
Group administrator
Thu 01 Mar 2018 05:36:34 PM UTC, comment #10: 

RE: comment #3, I believe that this configuration is always failing the test:

http://buildbot.octave.org:8010/builders/no-extras-debian

This is the one will all optional packages disabled.  Which one could affect this test specifically?  The configuration options are shown in this log file:

http://buildbot.octave.org:8010/builders/no-extras-debian/builds/1444/steps/configure/logs/stdio

John W. Eaton <jwe>
Group administrator
Thu 01 Mar 2018 11:39:43 AM UTC, comment #9: 

Actually not only hyperplanes but any linear subspaces of \mathbb{R}^n https://en.wikipedia.org/wiki/Flat_(geometry) (hyperplanes have dimension n-1).

Yes, as probably always there are two aspects to test: whether the algorithm is in principle correct, that is, disregarding inaccuracies and so on, the results make sense for all allowed input data, and second, whether the accuracy is okay. My tests address only the first part, and it would be good to also test the second part, where it is not so important to cover all qualitatively different inputs.

Michael Leitner <mleitner>
Thu 01 Mar 2018 10:39:27 AM UTC, comment #8: 

Ah, OK.  It's the angle between the hyperplanes formed by the vectors in each subspace.

Your tests are good--there is some higher dimensionality and some vectors in quadrants other than the first.

We do want to test some small angles though, which is the point of the higher accuracy algorithm.  (Yesterday I emailed the author listed in the script file for good tests regarding accuracy, but haven't heard anything.)  The problem with examples that are common fractions of pi is that often the trig function algorithm is perfectly accurate at that point because it could be a specially defined case rather than a numeric estimation.  An example of small angle might be


theta = pi/200;
Ry = [cos(theta) 0 sin(theta);0 1 0;-sin(theta) 0 cos(theta)];
assert(subspace(Ry*[1 0 0]',[1 1 0;1 -1 0]'), theta, eps);


Dan Sebald <sebald>
Thu 01 Mar 2018 08:39:11 AM UTC, comment #7: 

No, you misunderstood: The dimensionality of the subspaces can well be different, but this means that the number of columns would be different. However, they still have to be subspaces of a common vector space, specifically \mathbb{R}^n, thus they always need to have the same number of rows, in this case n. This is exactly what is checked in the input validation, and it also works:


subspace([1 1 0;1 -1 0]',[1 0 1]')


By the way, I think that this is something that is more worth being tested, i.e. to construct some such examples as above, with different column numbers of A and B, and row numbers larger than 2, and to compare with the exact results (of course pi/4 in the case above), rather than to test just the angles between two two-dimensional vectors, which is about the simplest conceivable task this very potent function can do. And also include degenerate test cases like


subspace([1,1)


I just took a few minutes and wrote the following tests


assert(subspace(1,1),0);
assert(subspace([1 0]',[1 1;0 1]'),0,3*eps);
assert(subspace([1 0 1]',[1 1 0;1 -1 0]'),pi/4,3*eps)
assert(subspace([1 5 0 0;-3 2 0 0]',[0 0 4 2;0 0 4 3]'),pi/2)
assert(subspace([1 1 1 1;1 2 3 4]',[1 -1 -1 1;]'),pi/2)


They should test also the more general behaviour and should not fail.

Michael Leitner <mleitner>
Wed 28 Feb 2018 06:38:54 PM UTC, comment #6: 

Here's the abstract for that article referenced in the subspace code:

http://epubs.siam.org/doi/pdf/10.1137/S1064827500377332


## reference:
## [1]  Andrew V. Knyazev, Merico E. Argentati:
## Principal Angles between Subspaces in an A-Based Scalar Product:
## Algorithms and Perturbation Estimates.
## SIAM Journal on Scientific Computing, Vol. 23 no. 6, pp. 2008-2040


One can see that the main point of the article is the small angle round-off errors.  So the test that was failing simply illustrates that problem with its use of acos(dot(a,b)/(a1*b1)).  Funny.  Anyway, we should be using the angle directly as in the changeset I provided.

As for the issue of different dimensionality A and B, would it be sufficient to simply extend the vectors by appending zeros?  I.e.,


--- a/scripts/linear-algebra/subspace.m
+++ b/scripts/linear-algebra/subspace.m
@@ -38,8 +38,10 @@ function ang = subspace (A, B)
     print_usage ();
   elseif (ndims (A) != 2 || ndims (B) != 2)
     error ("subspace: A and B must be 2-dimensional arrays");
-  elseif (rows (A) != rows (B))
-    error ("subspace: column dimensions of A and B must match");
+  elseif (rows (A) < rows (B))
+    A = [A; zeros(rows(B)-rows(A), size(A,2))];
+  elseif (rows (A) > rows (B))
+    B = [B; zeros(rows(A)-rows(B), size(B,2))];
   endif

   A = orth (A);


I'll think that one over.

Dan Sebald <sebald>
Wed 28 Feb 2018 05:08:12 PM UTC, comment #5: 

Something else I should point out.  Generally, the definition of subspace can be of different vector size, i.e., different space dimensionality (hence the word "sub"):

https://en.wikipedia.org/wiki/Angle#Angles_between_subspaces
https://en.wikipedia.org/wiki/Angles_between_flats

If I were to add a dimension to one of the vectors, say,

a(end+1)
subspace(a,b)

I think that the result should be the same.  But the current routine doesn't allow it:

error: subspace: column dimensions of A and B must match
error: called from
    subspace at line 42 column 5
    test_script3 at line 18 column 5

Dan Sebald <sebald>
Wed 28 Feb 2018 04:43:58 PM UTC, comment #4: 

Yes, I've noticed that the buildbots occasionally fail on this test, but had been too lazy to do something about it.

@Mike: Can you check Dan's patch in comment #3?

Rik <rik5>
Group administrator
Wed 28 Feb 2018 07:47:27 AM UTC, comment #3: 

Give the following changeset test a try.

As I see it, it doesn't make sense to compute theta in two different ways--one in the routine subspace() and one in the test script.  In the following hunk of code from subspace.m:


  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


It's basically computing exactly the same thing as the test script when scos^2 <= 1/2, i.e., ang = acos (scos).  One must assume the other case is to account for numerical issues with acos(), which of course the test script isn't doing.  So that right there is the flaw.  We could change the test script, but then it is as though we are replicating the function itself, which I don't see the point of.

A much more useful test seems to be starting with the angle known, then we are ostensibly testing the accuracy of the subspace() function, not two functions.  So, I used a rotation matrix (its angle theta is what's of interest) to generate a vector from the first random vector with some additional random scaling of the amplitude.  The angle goes around the unit circle so it is more than just the first quadrant that the vectors live in.  I also added some secondary jitter to be sure to test those often numerically problematic cases of near-0, near-pi, etc.

(file #43414)

Dan Sebald <sebald>
Wed 28 Feb 2018 04:14:46 AM UTC, comment #2: 

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)

Dan Sebald <sebald>
Wed 28 Feb 2018 01:11:40 AM UTC, comment #1: 

Edit: I just saw


***** test
 ## For random vectors
 a = rand (2,1);
 b = rand (2,1);
 a1 = norm (a,2);
 b1 = norm (b,2);
 theta = acos (dot (a,b)/(a1*b1));
 assert (theta, subspace (a, b), 100*eps);
!!!!! test failed
ASSERT errors for:  assert (theta,subspace (a, b),100 * eps)

  Location  |  Observed  |  Expected  |  Reason
     ()       4.7702e-05   4.7702e-05    Abs err 1.0036e-12 exceeds tol 2.2204e-14 by 1e-12


which is about 4520*eps.

Mike Miller <mtmiller>
Group Member
Wed 28 Feb 2018 01:08:32 AM UTC, original submission:  

The test of subspace with norm and acos fails relies on random data and fails randomly with a tolerance error some small percentage of the time. For example


***** test
 ## For random vectors
 a = rand (2,1);
 b = rand (2,1);
 a1 = norm (a,2);
 b1 = norm (b,2);
 theta = acos (dot (a,b)/(a1*b1));
 assert (theta, subspace (a, b), 100*eps);
!!!!! test failed
ASSERT errors for:  assert (theta,subspace (a, b),100 * eps)

  Location  |  Observed  |  Expected  |  Reason
     ()       0.0027191    0.0027191     Abs err 7.2089e-14 exceeds tol 2.2204e-14 by 5e-14


7.2089e-14 is the largest absolute error I've seen so far in testing this several hundred times, which is about 325*eps.

Any better guesses on a good upper bound on the tolerance?

Or do we want to seed the random generator to a fixed state like we have done for other unit tests?

Mike Miller <mtmiller>
Group Member

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:
   
   
Comment:
   

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by jwe (Posted a comment)
  • -email is unavailable- added by mleitner (Posted a comment)
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by sebald (Updated the item)
  • -email is unavailable- added by mtmiller (Submitted the item)
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only group members can vote.

     

    Follow 9 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2018-03-06 rik5 StatusIn Progress Fixed
        Open/ClosedOpen Closed
    2018-03-01 sebald Attached File- Added octave-subspace_angle_test_bug_53242-djs2018mar03.patch, #43431
    2018-03-01 sebald Attached File- Added octave-subspace_angle_test_bug_53242-djs2018mar02.patch, #43430
    2018-03-01 sebald Attached File- Added octave-subspace_angle_test_bug_53242-djs2018mar01.patch, #43429
    2018-03-01 mtmiller StatusConfirmed In Progress
    2018-02-28 rik5 StatusNone Confirmed
    2018-02-28 sebald Attached File- Added octave-subspace_angle_test_bug_53242-djs2018feb28.patch, #43414
    2018-02-28 sebald Attached File- Added test_script.m, #43411

    Back to the top

    Powered by Savane 3.13-02a9.
    Corresponding source code