Bugfix: Early exit in eigen_solve method

[jcs15c]

Summary

• This PR is a bugfix. It removes an unnecessary early termination criterion in the numerics::eigen_solve method.
• A test for the primal::OrientedBoundingBox( const PointType* pts ) constructor has been added, which utilizes this method.

The numerics::eigen_solve method uses the power method to find eigenvalues and eigenvectors for an NxN matrix. At each step, a random vector is created and then made orthonormal to previously found eigenvectors. If this randomly created vector is (numerically) already in the span of previously found eigenvectors, it has zero magnitude, and the algorithm will terminate early with an error status.

While this event has mathematically zero chance of occurring, it occurs numerically in roughly 0.05% of runs of the existing numerics_eigen_solve test, resulting in a fail state.

This behavior has been modified so that if the resulting vector has (near-) zero magnitude, it is instead discarded and a new one is generated.

Discussion

This behavior was originally observed through repeated calls to the primal::OrientedBoundingBox( const PointType* pts ) constructor as a part of a winding number algorithm, which would occasionally return erroneous results on Release or terminate entirely on Debug. Further investigation revealed that no test in the OBB testing suite used this constructor. Although the actual likelihood of failing such a test in this way is very low, as was the case with the existing numerics_eigen_solve tests, a test has been added.

[kennyweiss]

Thanks @jcs15c for this nice bugfix!

I have a few questions/concern though, possibly/likely outside the scope of this bugfix:

• if the point distribution truly has fewer than k dimensions, e.g. we're given 3D points that lie on a plane (2D affine space) or a line (1D affine space) and the user requests K=3 eigenvalues/vectors, will the code still find the right answer? Seems like it should, but perhaps we should add a test since one use case is to build an OBB over a triangle.
• similarly, there is a second normalize and check in step 3 -- does that always do the right thing for data that lives in a lower dimensional subspace? I think so too, but it might also be worth checking
• (almost definitely outside the scope of this PR): I noticed that we're only testing this function for k==N. Should we have tests for $0 < k<N$ or $k>N$? I suppose the higher level question is whether the parameter k is useful, or if it can always be assumed to be N
  • EDIT: I think we do want to keep the k parameter for eigen_solve, but I also think we should have tests for this. It does not have to be in this PR though. I created Improve testing of eigen_solve #1414 to track this.

Minor: There is a typo in the doxygen comment for eigen_solve:

- * \param [out] lambdas pointer to k eigenvales in order by size
+ * \param [out] lambdas pointer to k eigenvalues in order by size

Please also update the RELEASE_NOTES to reflect this bugfix

[kennyweiss]

Nice fix!

Bonus points: this was also a memory leak for temp!

[jcs15c]

@kennyweiss I've added some discussion of the eigen_solve method in the relevant issue, but your idea to check lower-dimensional subspaces for the OBB is a good one. Fortunately, degenerate axes of the OBB are already properly handled by the current implementation. I've updated the test accordingly.