Fix pca_numpy conditioning: SVD the centered data matrix, not the covariance matrix

CHANGELOG.md

@@ -15,6 +15,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Changed
 
+* Changed `pca_numpy` to run the SVD on the centered data matrix instead of the covariance matrix: forming the covariance squares the condition number, so near-degenerate inputs (e.g. almost collinear point clouds) returned wrong principal directions — `bestfit_plane_numpy` normals were off by 11–37 degrees on exactly planar sliver clouds (#1522). Well-conditioned results are unchanged (same eigenvectors; eigenvalues rescaled to keep their variance meaning).
 * Changed `Tolerance` class to no longer use singleton pattern. `Tolerance()` now creates independent instances instead of returning the global `TOL`. 
 * Renamed `Tolerance.units` to `Tolerance.unit` to better reflect the documented properties. Left `units` with deprecation warning.
 * Fixed `NotImplementedErorr` when calling `BrepLoop.vertices`.

src/compas/geometry/pca_numpy.py

@@ -54,34 +54,32 @@ def pca_numpy(data):
     # across all observations
     Y = X - mean
 
-    # covariance matrix of spread
-    # note: there is a covariance function in NumPy...
-    # the shape of the covariance matrix is dim x dim
-    # for example, if the data are 2D point coordinates, the shape of C is 2 x 2
-    # the diagonal of the covariance matrix contains the variance of each variable
-    # the off-diagonal elements of the covariance matrix contain the covariance
-    # of two independent variables
-    C = Y.T.dot(Y) / (n - 1)
-
-    # assert C.shape[0] == dim, "The shape of the covariance matrix is not correct."
-
-    # SVD of covariance matrix
-    u, s, vT = svd(C, full_matrices=False)
+    # SVD of the spread matrix directly — deliberately NOT of the covariance
+    # matrix C = Y.T @ Y / (n - 1).
+    # The right-singular vectors of Y are exactly the eigenvectors of C, but
+    # forming C squares the condition number: for near-degenerate data (e.g.
+    # an almost collinear point cloud, whose smallest extent is ~1e-8 of its
+    # largest) the smallest principal direction drowns in floating-point
+    # rounding and the returned directions are wrong, while the SVD of Y
+    # recovers them to full precision. See issue #1522.
+    u, s, vT = svd(Y, full_matrices=False)
 
     # eigenvectors
     # ------------
     # note: the eigenvectors are normalized
     # note: vT is exactly what it says it will be => the transposed eigenvectors
     # => take the rows of vT, or the columns of v
-    # the right-singular vectors of C (the columns of V or the rows of Vt)
-    # are the eigenvectors of CtC
+    # the right-singular vectors of Y (the columns of V or the rows of Vt)
+    # are the eigenvectors of C = Y.T @ Y / (n - 1)
     eigenvectors = vT
 
     # eigenvalues
     # -----------
-    # the nonzero singular values of C are the square roots
-    # of the nonzero eigenvalues of CtC and CCt
-    eigenvalues = s
+    # the singular values of Y are the square roots of the (n - 1)-scaled
+    # eigenvalues of the covariance matrix C; rescale so the returned values
+    # keep their meaning (the variance of the data along each principal
+    # direction), identical to what the SVD of C used to return
+    eigenvalues = (s**2) / (n - 1)
 
     # return
     return mean[0], eigenvectors, eigenvalues

tests/compas/geometry/test_pca_numpy.py

@@ -0,0 +1,88 @@
+import compas
+
+
+def _sliver_points():
+    """An exactly planar but near-collinear cloud: long axis ~1, short in-plane
+    axis ~1e-8, rotated off-axis. The plane (and its normal) are perfectly
+    well defined; only an ill-conditioned fit loses them."""
+    import numpy as np
+
+    rng = np.random.default_rng(7)
+    n = 200
+    t = rng.uniform(-1.0, 1.0, n)
+    s = rng.uniform(-1e-8, 1e-8, n)
+    pts_local = np.column_stack([t, s, np.zeros(n)])
+
+    a, b = 0.6, -1.1
+    Rx = np.array([[1, 0, 0], [0, np.cos(a), -np.sin(a)], [0, np.sin(a), np.cos(a)]])
+    Rz = np.array([[np.cos(b), -np.sin(b), 0], [np.sin(b), np.cos(b), 0], [0, 0, 1]])
+    R = Rz.dot(Rx)
+    points = pts_local.dot(R.T) + np.array([100.0, -40.0, 7.0])
+    normal = R.dot(np.array([0.0, 0.0, 1.0]))
+    return points, normal
+
+
+def test_pca_numpy_well_conditioned_matches_covariance_route():
+    if compas.IPY:
+        return
+
+    import numpy as np
+
+    from compas.geometry import pca_numpy
+
+    rng = np.random.default_rng(42)
+    points = rng.uniform(-10.0, 10.0, (100, 3))
+
+    mean, eigenvectors, eigenvalues = pca_numpy(points)
+
+    # eigenvalues must still be the variances along the principal directions
+    Y = points - points.mean(axis=0)
+    C = Y.T.dot(Y) / (len(points) - 1)
+    expected = np.sort(np.linalg.eigvalsh(C))[::-1]
+    assert np.allclose(eigenvalues, expected)
+
+    # eigenvectors diagonalize the covariance matrix
+    V = np.asarray(eigenvectors)
+    assert np.allclose(V.dot(C).dot(V.T), np.diag(eigenvalues), atol=1e-12)
+
+    assert np.allclose(mean, points.mean(axis=0))
+
+
+def test_pca_numpy_near_collinear_recovers_smallest_direction():
+    if compas.IPY:
+        return
+
+    import numpy as np
+
+    from compas.geometry import pca_numpy
+
+    points, normal = _sliver_points()
+
+    _, eigenvectors, eigenvalues = pca_numpy(points)
+
+    # the smallest principal direction is the plane normal (the cloud is
+    # exactly planar); with the covariance route this was off by ~12 degrees
+    smallest = np.asarray(eigenvectors)[2]
+    angle = np.degrees(np.arccos(np.clip(abs(smallest.dot(normal)), -1.0, 1.0)))
+    assert angle < 1e-4
+
+    # eigenvalues are returned in descending order and non-negative
+    assert eigenvalues[0] >= eigenvalues[1] >= eigenvalues[2] >= 0.0
+
+
+def test_bestfit_plane_numpy_near_collinear():
+    if compas.IPY:
+        return
+
+    import numpy as np
+
+    from compas.geometry import bestfit_plane_numpy
+
+    points, normal = _sliver_points()
+
+    _, fitted_normal = bestfit_plane_numpy(points)
+
+    fitted = np.asarray(fitted_normal, dtype=float)
+    fitted = fitted / np.linalg.norm(fitted)
+    angle = np.degrees(np.arccos(np.clip(abs(fitted.dot(normal)), -1.0, 1.0)))
+    assert angle < 1e-4