From e9e510335050d231f3ce75075147c89e9f3a6f39 Mon Sep 17 00:00:00 2001 From: Zachary Ferguson Date: Sat, 18 Jul 2026 14:25:40 -0400 Subject: [PATCH 1/5] Fix missing PSD projection of mollified m<=0 hessian blocks MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit NormalPotential::hessian() early-returns the mollified block (weight·f)·∇²m when the mollifier m == 0 (exactly parallel edges), but did so WITHOUT PSD projection while every other path projects. For positive weights the block is PSD so this was harmless, but IMPROVED_MAX_APPROX produces negative-weight collisions, making the block negative-(semi)definite and the assembled "PSD-projected" hessian non-PSD. Project the block like the other paths. Since m == 0 is a global minimum of the mollifier, ∇²m is PSD and f = f(d) > 0, so the block is a scalar multiple of a PSD matrix and its projection reduces to projecting the scalar weight·f (clamp/abs) -- no eigendecomposition needed. Add a regression test asserting the assembled CLAMP hessian is PSD for a cube under IMPROVED_MAX_APPROX (which exercises the m == 0 branch). Co-Authored-By: Claude Opus 4.8 --- src/ipc/potentials/normal_potential.cpp | 16 +++++++- .../potential/test_barrier_potential.cpp | 41 +++++++++++++++++++ 2 files changed, 56 insertions(+), 1 deletion(-) diff --git a/src/ipc/potentials/normal_potential.cpp b/src/ipc/potentials/normal_potential.cpp index 079cbd3c3..de116e68c 100644 --- a/src/ipc/potentials/normal_potential.cpp +++ b/src/ipc/potentials/normal_potential.cpp @@ -8,6 +8,9 @@ #include #include +#include +#include + namespace ipc { // -- Cumulative methods ------------------------------------------------------- @@ -183,9 +186,20 @@ MatrixMax12d NormalPotential::hessian( // m(x) before evaluating barrier derivatives. const double m = collision.mollifier(positions); // m(x) if (collision.is_mollified() && m <= 0) { + // The mollified hessian is w·f·∇²m here (∇m = 0 when m = 0). Because + // m == 0 (⟺ ‖cross‖² == 0, exactly parallel edges) is a global minimum + // of the mollifier, hess_m = ∇²m is PSD and f = f(d) > 0. The block is + // therefore a scalar multiple of a PSD matrix, so its PSD projection + // reduces to projecting the scalar w·f -- no eigendecomposition needed. const double f = (*this)(d, collision.dmin); const MatrixMax12d hess_m = collision.mollifier_hessian(positions); - return (collision.weight * f) * hess_m; + double scale = collision.weight * f; + if (project_hessian_to_psd == PSDProjectionMethod::CLAMP) { + scale = std::max(scale, 0.0); // NSD (w < 0) projects to zero + } else if (project_hessian_to_psd == PSDProjectionMethod::ABS) { + scale = std::abs(scale); + } + return scale * hess_m; } // ∇d(x) diff --git a/tests/src/tests/potential/test_barrier_potential.cpp b/tests/src/tests/potential/test_barrier_potential.cpp index 2bb56f509..e90a816dd 100644 --- a/tests/src/tests/potential/test_barrier_potential.cpp +++ b/tests/src/tests/potential/test_barrier_potential.cpp @@ -466,6 +466,47 @@ TEST_CASE( } } +TEST_CASE( + "Mollified m<=0 hessian block is PSD", + "[potential][barrier_potential][hessian]") +{ + // The cube's edges are axis-aligned, so IMPROVED_MAX_APPROX's near-parallel + // negative-weight collisions have mollifier m == 0 exactly. That exercises + // NormalPotential::hessian()'s early-return branch, which returns the block + // (weight * f) * hess_m WITHOUT PSD projection. In debug builds an assert + // there checks the returned block is actually PSD -- if it fires, the CPU + // is returning a non-PSD block unprojected (a latent bug). + Eigen::MatrixXd vertices; + Eigen::MatrixXi edges, faces; + REQUIRE(tests::load_mesh("cube.ply", vertices, edges, faces)); + const CollisionMesh mesh(vertices, edges, faces); + + const double dhat = std::sqrt(2.0); + NormalCollisions collisions; + collisions.set_use_area_weighting(true); + collisions.set_collision_set_type( + NormalCollisions::CollisionSetType::IMPROVED_MAX_APPROX); + collisions.build(mesh, vertices, dhat); + REQUIRE(!collisions.empty()); + + const BarrierPotential barrier_potential(dhat, /*stiffness=*/1.0); + + // Triggers the m <= 0 early-return (and its debug PSD assert) for every + // mollified collision that reaches it. + const Eigen::SparseMatrix H = barrier_potential.hessian( + collisions, mesh, vertices, PSDProjectionMethod::CLAMP); + REQUIRE(H.nonZeros() > 0); + + // With per-block CLAMP projection (including the m <= 0 blocks), the + // assembled hessian must be PSD. Before the fix, the m <= 0 branch returned + // negative-definite blocks unprojected, making this fail. + const Eigen::MatrixXd H_dense(H); + const Eigen::SelfAdjointEigenSolver eigensolver(H_dense); + const double min_eig = eigensolver.eigenvalues().minCoeff(); + const double scale = eigensolver.eigenvalues().cwiseAbs().maxCoeff(); + CHECK(min_eig >= -1e-10 * std::max(scale, 1.0)); +} + // -- Benchmarking ------------------------------------------------------------ TEST_CASE( From 9adf40e0181da37f6840b7934dfb0d7c7eac2221 Mon Sep 17 00:00:00 2001 From: Zachary Ferguson Date: Sat, 18 Jul 2026 14:32:47 -0400 Subject: [PATCH 2/5] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- tests/src/tests/potential/test_barrier_potential.cpp | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/tests/src/tests/potential/test_barrier_potential.cpp b/tests/src/tests/potential/test_barrier_potential.cpp index e90a816dd..e2fa5f26b 100644 --- a/tests/src/tests/potential/test_barrier_potential.cpp +++ b/tests/src/tests/potential/test_barrier_potential.cpp @@ -472,10 +472,11 @@ TEST_CASE( { // The cube's edges are axis-aligned, so IMPROVED_MAX_APPROX's near-parallel // negative-weight collisions have mollifier m == 0 exactly. That exercises - // NormalPotential::hessian()'s early-return branch, which returns the block - // (weight * f) * hess_m WITHOUT PSD projection. In debug builds an assert - // there checks the returned block is actually PSD -- if it fires, the CPU - // is returning a non-PSD block unprojected (a latent bug). + // NormalPotential::hessian()'s m <= 0 early-return branch. + // + // Before the fix, that branch returned (weight * f) * hess_m without applying + // the requested PSD projection, so negative-weight collisions could introduce + // non-PSD blocks into the assembled "PSD-projected" hessian. Eigen::MatrixXd vertices; Eigen::MatrixXi edges, faces; REQUIRE(tests::load_mesh("cube.ply", vertices, edges, faces)); From 532c9332e0198d8c8fb159091e7c71d8a3c617a3 Mon Sep 17 00:00:00 2001 From: Zachary Ferguson Date: Sat, 18 Jul 2026 14:33:29 -0400 Subject: [PATCH 3/5] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- src/ipc/potentials/normal_potential.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/ipc/potentials/normal_potential.cpp b/src/ipc/potentials/normal_potential.cpp index de116e68c..c5761e0d1 100644 --- a/src/ipc/potentials/normal_potential.cpp +++ b/src/ipc/potentials/normal_potential.cpp @@ -188,9 +188,9 @@ MatrixMax12d NormalPotential::hessian( if (collision.is_mollified() && m <= 0) { // The mollified hessian is w·f·∇²m here (∇m = 0 when m = 0). Because // m == 0 (⟺ ‖cross‖² == 0, exactly parallel edges) is a global minimum - // of the mollifier, hess_m = ∇²m is PSD and f = f(d) > 0. The block is - // therefore a scalar multiple of a PSD matrix, so its PSD projection - // reduces to projecting the scalar w·f -- no eigendecomposition needed. + // of the mollifier, hess_m = ∇²m is PSD. The block is therefore a scalar + // multiple of a PSD matrix, so its PSD projection reduces to projecting + // the scalar w·f -- no eigendecomposition needed. const double f = (*this)(d, collision.dmin); const MatrixMax12d hess_m = collision.mollifier_hessian(positions); double scale = collision.weight * f; From d910dd42e089281015f042fa987736fcebfa8603 Mon Sep 17 00:00:00 2001 From: Zachary Ferguson Date: Sat, 18 Jul 2026 14:34:19 -0400 Subject: [PATCH 4/5] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- .../potential/test_barrier_potential.cpp | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) diff --git a/tests/src/tests/potential/test_barrier_potential.cpp b/tests/src/tests/potential/test_barrier_potential.cpp index e2fa5f26b..977ba8c28 100644 --- a/tests/src/tests/potential/test_barrier_potential.cpp +++ b/tests/src/tests/potential/test_barrier_potential.cpp @@ -490,6 +490,25 @@ TEST_CASE( collisions.build(mesh, vertices, dhat); REQUIRE(!collisions.empty()); + // Ensure the test actually exercises NormalPotential::hessian()'s m <= 0 + // early-return branch (and the negative-weight case that motivated the fix). + int m_le_0_count = 0; + int m_le_0_negative_weight_count = 0; + for (const auto& c : collisions) { + if (!c.is_mollified()) { + continue; + } + const double m = c.mollifier(c.dof(vertices, mesh.edges(), mesh.faces())); + if (m <= 0) { + m_le_0_count++; + if (c.weight < 0) { + m_le_0_negative_weight_count++; + } + } + } + REQUIRE(m_le_0_count > 0); + REQUIRE(m_le_0_negative_weight_count > 0); + const BarrierPotential barrier_potential(dhat, /*stiffness=*/1.0); // Triggers the m <= 0 early-return (and its debug PSD assert) for every From f9e4ef21e9ce80bf5bba0811e17956147fdefe44 Mon Sep 17 00:00:00 2001 From: Zachary Ferguson Date: Sat, 18 Jul 2026 15:23:49 -0400 Subject: [PATCH 5/5] Fix CI: formatting and test build errors - Reword the m<=0 comment in ASCII with shorter lines so clang-format 20 (used in CI) leaves it unchanged. - test_barrier_potential.cpp: NormalCollisions has no begin()/end(), so the range-based for over collisions did not compile; use an index loop over operator[]/size(). Also include for SelfAdjointEigenSolver and use SparseMatrix::toDense() instead of the sparse->dense constructor. Co-Authored-By: Claude Opus 4.8 --- src/ipc/potentials/normal_potential.cpp | 11 +++++----- .../potential/test_barrier_potential.cpp | 20 ++++++++++++------- 2 files changed, 19 insertions(+), 12 deletions(-) diff --git a/src/ipc/potentials/normal_potential.cpp b/src/ipc/potentials/normal_potential.cpp index c5761e0d1..96b1b0e6b 100644 --- a/src/ipc/potentials/normal_potential.cpp +++ b/src/ipc/potentials/normal_potential.cpp @@ -186,11 +186,12 @@ MatrixMax12d NormalPotential::hessian( // m(x) before evaluating barrier derivatives. const double m = collision.mollifier(positions); // m(x) if (collision.is_mollified() && m <= 0) { - // The mollified hessian is w·f·∇²m here (∇m = 0 when m = 0). Because - // m == 0 (⟺ ‖cross‖² == 0, exactly parallel edges) is a global minimum - // of the mollifier, hess_m = ∇²m is PSD. The block is therefore a scalar - // multiple of a PSD matrix, so its PSD projection reduces to projecting - // the scalar w·f -- no eigendecomposition needed. + // The mollified hessian is weight * f * hess_m here (the mollifier + // gradient is zero when m = 0). At m == 0 the edges are exactly + // parallel -- a global minimum of the mollifier -- so hess_m is PSD, + // and f = f(d) > 0. The block is therefore a scalar multiple of a PSD + // matrix, so its PSD projection reduces to projecting the scalar + // weight * f (no eigendecomposition needed). const double f = (*this)(d, collision.dmin); const MatrixMax12d hess_m = collision.mollifier_hessian(positions); double scale = collision.weight * f; diff --git a/tests/src/tests/potential/test_barrier_potential.cpp b/tests/src/tests/potential/test_barrier_potential.cpp index 977ba8c28..5a82f0403 100644 --- a/tests/src/tests/potential/test_barrier_potential.cpp +++ b/tests/src/tests/potential/test_barrier_potential.cpp @@ -15,6 +15,8 @@ #include #include +#include + using namespace ipc; TEST_CASE( @@ -474,9 +476,10 @@ TEST_CASE( // negative-weight collisions have mollifier m == 0 exactly. That exercises // NormalPotential::hessian()'s m <= 0 early-return branch. // - // Before the fix, that branch returned (weight * f) * hess_m without applying - // the requested PSD projection, so negative-weight collisions could introduce - // non-PSD blocks into the assembled "PSD-projected" hessian. + // Before the fix, that branch returned (weight * f) * hess_m without + // applying the requested PSD projection, so negative-weight collisions + // could introduce non-PSD blocks into the assembled "PSD-projected" + // hessian. Eigen::MatrixXd vertices; Eigen::MatrixXi edges, faces; REQUIRE(tests::load_mesh("cube.ply", vertices, edges, faces)); @@ -491,14 +494,17 @@ TEST_CASE( REQUIRE(!collisions.empty()); // Ensure the test actually exercises NormalPotential::hessian()'s m <= 0 - // early-return branch (and the negative-weight case that motivated the fix). + // early-return branch (and the negative-weight case that motivated the + // fix). int m_le_0_count = 0; int m_le_0_negative_weight_count = 0; - for (const auto& c : collisions) { + for (size_t i = 0; i < collisions.size(); i++) { + const NormalCollision& c = collisions[i]; if (!c.is_mollified()) { continue; } - const double m = c.mollifier(c.dof(vertices, mesh.edges(), mesh.faces())); + const double m = + c.mollifier(c.dof(vertices, mesh.edges(), mesh.faces())); if (m <= 0) { m_le_0_count++; if (c.weight < 0) { @@ -520,7 +526,7 @@ TEST_CASE( // With per-block CLAMP projection (including the m <= 0 blocks), the // assembled hessian must be PSD. Before the fix, the m <= 0 branch returned // negative-definite blocks unprojected, making this fail. - const Eigen::MatrixXd H_dense(H); + const Eigen::MatrixXd H_dense = H.toDense(); const Eigen::SelfAdjointEigenSolver eigensolver(H_dense); const double min_eig = eigensolver.eigenvalues().minCoeff(); const double scale = eigensolver.eigenvalues().cwiseAbs().maxCoeff();