Bug fixes:
1. computeDirect 3x3: eigenvalues from computeRoots() are theoretically
sorted via the trigonometric formula, but floating-point rounding
(especially in float) can break the ordering. Add a 3-element sorting
network after computeRoots() to guarantee sorted output.
2. compute(): Inf input silently returns Success with garbage results,
unlike NaN which correctly returns NoConvergence. Add an isfinite()
check on the scaling factor (which is maxCoeff of the matrix) to
detect Inf/NaN early and return NoConvergence.
3. computeFromTridiagonal(): lacks the equilibration scaling that
compute() applies, causing overflow and NoConvergence for tridiagonal
matrices with large entries. Add the same scale-to-[-1,1] pattern.
New tests:
- Eigenvalue sorting verification (both iterative and direct solvers)
- Repeated/degenerate eigenvalues (all equal, multiplicity n-1, two
clusters, nearly repeated separated by O(epsilon))
- Extreme eigenvalue ranges (high condition number spanning many orders
of magnitude, near-underflow, near-overflow, mixed positive/negative,
rank-deficient with zero eigenvalue)
- computeFromTridiagonal with large and tiny values
- Diagonal matrices (eigenvalues must match sorted diagonal)
- operatorInverseSqrt accuracy (sqrtA*invSqrtA=I, invSqrtA*A*invSqrtA=I,
symmetry)
- RowMajor storage for computeDirect (2x2, 3x3, float, double) and
iterative solver (dynamic RowMajor)
- Inf input detection
- Tridiagonal structure verification (off-tridiagonal entries are zero)
- Direct solver stress tests: 3x3 (near-planar covariance, triple
eigenvalue, double eigenvalue, large off-diagonal, nearly singular)
and 2x2 (equal eigenvalues, tiny off-diagonal, huge diagonal ratio,
anti-diagonal dominant, negative entries)
- Tightened unitary tolerance from fixed 32*test_precision to
4*n*epsilon (scales with matrix size)
- Fixed typo: "expponential" -> "exponential"
Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>
Fix a bug in the 3x3 direct eigensolver's Gram-Schmidt orthogonalization
for near-degenerate eigenvalues. The code was subtracting the projection
onto eivecs.col(l) (itself) instead of onto eivecs.col(k):
// Before (bug): subtracts scalar multiple of self — does nothing useful
eivecs.col(l) -= eivecs.col(k).dot(eivecs.col(l)) * eivecs.col(l);
// After (fix): removes component along eivecs.col(k)
eivecs.col(l) -= eivecs.col(k).dot(eivecs.col(l)) * eivecs.col(k);
This path is taken when two of three eigenvalues are nearly equal, which
is common for covariance matrices of near-planar point clouds.
Also add comprehensive small fixed-size matrix benchmarks covering the
operations that dominate robotics/CV inner loops: matmul, matvec,
inverse, determinant, LLT, LDLT, PartialPivLU, ColPivHouseholderQR,
JacobiSVD, SelfAdjointEigenSolver (iterative and direct) for sizes
2x2 through 8x9.
Note: the direct 3x3 eigensolver (computeDirect) is 3x faster than
the iterative solver but has 5-6 orders of magnitude worse residuals
for near-degenerate eigenvalues. This is inherent to the closed-form
algorithm, not a consequence of the Gram-Schmidt bug. Users should
prefer compute() when accuracy matters and computeDirect() only when
speed is critical and eigenvalues are well-separated.
Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>
