Repository navigation
Conversation
SymmetricEigen3::new always computed the first eigenvector from the smallest eigenvalue with eigenvector1, which requires a multiplicity of 1. When the smallest eigenvalue is repeated (e.g. the inertia tensor of an axisymmetric body), its eigenvector is dominated by rounding noise, so the decomposition of a rotated such matrix doesn't reconstruct it. Start from the extreme eigenvalue that is the farthest from the middle one, as in Eberly's algorithm.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Problem
SymmetricEigen3::newreturns wrong eigenvectors for a rotated matrix whose smallest eigenvalue is repeated. Its eigenvalues are correct, butV · diag(λ) · Vᵀdoesn't reconstruct the input.For example,
R · diag(0.02, 0.02, 0.035) · RᵀwithRa rotation of 0.9 rad around(1, 2, 3)is reconstructed with a ~24% relative error. A distinct or fully repeated spectrum is fine.This is the inertia tensor of any axisymmetric body. Downstream,
parry3d::MassProperties::with_inertia_matrix(used when summing the mass properties of several colliders) gets wrong principal axes. In rapier3d 0.36, adding a zero-mass collider to a rigid-body whose only mass carrier has an axisymmetric, rotated inertia changes the body's inertia tensor by ~30%. rapier3d-compat 0.20, which still used nalgebra's eigensolver, was off by ~5e-5 on the same input.Cause
newalways computes the first eigenvector from the smallest eigenvalue witheigenvector1, which (as its documentation says) requires an eigenvalue with a multiplicity of 1. When the smallest eigenvalue is repeated,mat - λIhas rank 1, so the cross productseigenvector1relies on are pure rounding noise and give an arbitrary direction. Without a rotation the noise is exactly zero and thed_max < 1e-20fallback happens to return the right axis, which is why the diagonal tests pass.Fix
As in Eberly's algorithm, start from the extreme eigenvalue that is the farthest from the middle one, which is the one with a multiplicity of 1 when two eigenvalues are equal:
λ1 - λ0 >= λ2 - λ1:v0 = eigenvector1(λ0),v1 = eigenvector2(v0, λ1),v2 = v0 × v1(unchanged),v2 = eigenvector1(λ2),v1 = eigenvector2(v2, λ1),v0 = v1 × v2.Eigenvalues stay in ascending order, eigenvector columns match them, and the basis stays right-handed.
Tests
eigen_3x3_repeated_eigenvalue: rotated matrices with a repeated smallest or largest eigenvalue, forSymmetricEigen3andSymmetricEigen3A. It checks the eigenvalues, orthonormality, right-handedness and the reconstruction. It fails without the fix and passes with it.eigen_3x3hard-coded the sign of each eigenvector, which isn't defined. With this change, the first and last ones come out negated (still a valid right-handed basis), so the test now compares each eigenvector up to its sign and checks the determinant.cargo test --all-features,cargo test --no-default-features --features libm,cargo clippy --all-features -- -D warningsandcargo fmt --checkpass.MassProperties::with_inertia_matrixon the example above goes from a 2.4e-1 to a 1.1e-7 reconstruction error. In a rapier3d JS build, the rigid-body inertia error mentioned above goes from 3e-1 to 1e-7.AI assistance (Claude) was used to investigate, write the fix and the tests, and I reviewed the results.