From 0b25f0a2c4ac75e830c057fba72690a83cf8381e Mon Sep 17 00:00:00 2001 From: Mengxi He Date: Mon, 6 Jul 2026 19:21:35 +0200 Subject: [PATCH 1/2] Fix pca_numpy conditioning: SVD the centered data, not the covariance 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 #1522 --- CHANGELOG.md | 1 + src/compas/geometry/pca_numpy.py | 34 ++++++------ tests/compas/geometry/test_pca_numpy.py | 70 +++++++++++++++++++++++++ 3 files changed, 87 insertions(+), 18 deletions(-) create mode 100644 tests/compas/geometry/test_pca_numpy.py diff --git a/CHANGELOG.md b/CHANGELOG.md index 48f846510721..51970fe268b3 100644 --- a/CHANGELOG.md +++ b/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`. diff --git a/src/compas/geometry/pca_numpy.py b/src/compas/geometry/pca_numpy.py index 9e2fcf13d3f4..e758ae82d8b3 100644 --- a/src/compas/geometry/pca_numpy.py +++ b/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 diff --git a/tests/compas/geometry/test_pca_numpy.py b/tests/compas/geometry/test_pca_numpy.py new file mode 100644 index 000000000000..8d059974c235 --- /dev/null +++ b/tests/compas/geometry/test_pca_numpy.py @@ -0,0 +1,70 @@ +import numpy as np +import pytest + +from compas.geometry import bestfit_plane_numpy +from compas.geometry import pca_numpy + + +@pytest.fixture +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.""" + 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 @ Rx + points = pts_local @ R.T + np.array([100.0, -40.0, 7.0]) + normal = R @ np.array([0.0, 0.0, 1.0]) + return points, normal + + +def test_pca_numpy_well_conditioned_matches_covariance_route(): + 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 @ 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 @ C @ 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(sliver_points): + 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 @ 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(sliver_points): + points, normal = sliver_points + + _, fitted_normal = bestfit_plane_numpy(points) + + fitted = np.asarray(fitted_normal, dtype=float) + fitted /= np.linalg.norm(fitted) + angle = np.degrees(np.arccos(np.clip(abs(fitted @ normal), -1.0, 1.0))) + assert angle < 1e-4 From 3f7ba38795919318d77745ef69c996a8d9b15015 Mon Sep 17 00:00:00 2001 From: Mengxi He Date: Tue, 7 Jul 2026 09:23:14 +0200 Subject: [PATCH 2/2] Guard test_pca_numpy for IronPython (py2-parseable, compas.IPY skips, no matmul operator) --- tests/compas/geometry/test_pca_numpy.py | 56 ++++++++++++++++--------- 1 file changed, 37 insertions(+), 19 deletions(-) diff --git a/tests/compas/geometry/test_pca_numpy.py b/tests/compas/geometry/test_pca_numpy.py index 8d059974c235..56333bb1cc3c 100644 --- a/tests/compas/geometry/test_pca_numpy.py +++ b/tests/compas/geometry/test_pca_numpy.py @@ -1,15 +1,12 @@ -import numpy as np -import pytest +import compas -from compas.geometry import bestfit_plane_numpy -from compas.geometry import pca_numpy - -@pytest.fixture -def sliver_points(): +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) @@ -19,13 +16,20 @@ def sliver_points(): 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 @ Rx - points = pts_local @ R.T + np.array([100.0, -40.0, 7.0]) - normal = R @ np.array([0.0, 0.0, 1.0]) + 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)) @@ -33,38 +37,52 @@ def test_pca_numpy_well_conditioned_matches_covariance_route(): # eigenvalues must still be the variances along the principal directions Y = points - points.mean(axis=0) - C = Y.T @ Y / (len(points) - 1) + 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 @ C @ V.T, np.diag(eigenvalues), atol=1e-12) + 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(sliver_points): - points, normal = sliver_points +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 @ normal), -1.0, 1.0))) + 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(sliver_points): - points, normal = sliver_points +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 /= np.linalg.norm(fitted) - angle = np.degrees(np.arccos(np.clip(abs(fitted @ normal), -1.0, 1.0))) + 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