Fix pca_numpy conditioning: SVD the centered data matrix, not the covariance matrix - #1524
Conversation
… matrix Forming C = Y.T @ Y / (n - 1) squares the condition number, so near- degenerate inputs (almost collinear point clouds) lose their smallest principal direction to floating-point rounding - bestfit_plane_numpy normals came back 11-37 degrees off on exactly planar sliver clouds. The right-singular vectors of Y are the same eigenvectors; eigenvalues are rescaled (s**2 / (n - 1)) to keep their variance meaning, so well- conditioned results are unchanged. Fixes compas-dev#1522
… no matmul operator)
|
Fixed the |
|
@tomvanmele ping ;) |
tomvanmele
left a comment
There was a problem hiding this comment.
nice, thanks for correcting this!
LGTM
|
@mengxihex perhaps add yourself to the author list :) |
|
thanks a lot @tomvanmele! Long-time admirer of compas, the abstraction layer AEC has been missing to scale. Glad to contribute. |
|
btw @mengxihex i am currently modernizing the codebase in this PR #1512 |
Fixes #1522.
What
pca_numpyformed the covariance matrixC = Y.T @ Y / (n - 1)and SVD'd it. FormingCsquares the condition number, so for near-degenerate inputs (an almost collinear point cloud whose smallest extent is ~1e-8 of its largest — still an exactly planar, perfectly posed fit) the smallest principal direction drowns at machine epsilon and comes back wrong:bestfit_plane_numpynormals were off by 11.6° (1e-8 aspect) to 37.3° (1e-9 aspect) on the reproduction in #1522.This PR runs the SVD on the centered data matrix
Ydirectly. The right-singular vectors ofYare exactly the eigenvectors ofC, and the singular values are rescaled (s**2 / (n - 1)) so the returned eigenvalues keep their meaning (variance along each principal direction). Well-conditioned inputs return the same results as before; only the ill-conditioned ones improve. Cost is equivalent (both routes are O(n·dim²)).Checks
tests/compas/geometry/test_pca_numpy.py: (1) well-conditioned parity with the covariance route (eigenvalues matcheigvalsh(C), eigenvectors diagonalizeC), (2) the near-collinear cloud recovers its smallest direction to <1e-4°, (3)bestfit_plane_numpyregression on the same cloud (fails on the previous implementation with ~11.6° error).tests/compas/geometrysuite: 519 passed (includes thepca_numpyconsumersbestfit_numpy,bbox_numpy,icp_numpy).One behavioural note for review: SVD sign indeterminacy means individual principal directions may flip sign relative to the old implementation on some inputs (the sign was already unspecified in both routes).