From 71072c025d21fa48b429287b613d1cb9b401ed34 Mon Sep 17 00:00:00 2001 From: Dhruvil Date: Sun, 17 May 2026 13:54:47 -0400 Subject: [PATCH] Honor lower-triangle contract for H in mju_boxQP. mju_boxQP documents that only the lower triangle of the Hessian H is read, but the gradient and search-direction updates inside mju_boxQPoption still called the dense mju_mulMatVec, which reads the upper triangle as well. This violated the documented contract and prevented callers from safely leaving the upper triangle uninitialized. Add a file-local mulMatVecSym helper that computes res = H*vec while reading only the lower triangle of H (mirroring the convention of the existing mulVecMatVecSym quadratic-form helper), and use it in place of mju_mulMatVec in both call sites. Extend the BoxQP test suite with UpperTrianglePoisoned, which fills the strict upper triangle of H with NaN and verifies that the solver produces the same result as on the clean symmetric input. Reported by @lshdlut. Fixes #3275 --- src/engine/engine_util_solve.c | 19 +++++++++- test/engine/engine_util_solve_test.cc | 53 +++++++++++++++++++++++++++ 2 files changed, 70 insertions(+), 2 deletions(-) diff --git a/src/engine/engine_util_solve.c b/src/engine/engine_util_solve.c index a6e5541c..6da65ac0 100644 --- a/src/engine/engine_util_solve.c +++ b/src/engine/engine_util_solve.c @@ -1407,6 +1407,21 @@ static mjtNum mulVecMatVecSym(const mjtNum* vec, const mjtNum* mat, int n) { } +// multiply symmetric matrix with vector: res = mat*vec +// assumes symmetry of mat, ignores upper triangle +// res must not alias vec +static void mulMatVecSym(mjtNum* res, const mjtNum* mat, const mjtNum* vec, int n) { + for (int i=0; i < n; i++) { + // diagonal + strict lower triangle: res[i] = sum_{j<=i} mat[i,j] * vec[j] + res[i] = mat[n*i+i] * vec[i] + mju_dot(mat+n*i, vec, i); + // strict upper mirror contribution: res[k] += mat[i,k] * vec[i] for k < i + for (int k=0; k < i; k++) { + res[k] += mat[n*i+k] * vec[i]; + } + } +} + + // minimize 0.5*x'*H*x + x'*g s.t. lower <= x <=upper, explicit options // additional arguments to mju_boxQP (see mju_boxQP documentation): // maxiter maximum number of iterations @@ -1512,7 +1527,7 @@ int mju_boxQPoption(mjtNum* res, mjtNum* R, int* index, // outputs oldvalue = value; // compute gradient - mju_mulMatVec(grad, H, res, n, n); + mulMatVecSym(grad, H, res, n); mju_addTo(grad, g, n); // find clamped dimensions @@ -1555,7 +1570,7 @@ int mju_boxQPoption(mjtNum* res, mjtNum* R, int* index, // outputs for (int i=0; i < n; i++) { temp[i] = clamped[i] ? res[i] : 0; } - mju_mulMatVec(search, H, temp, n, n); + mulMatVecSym(search, H, temp, n); mju_addTo(search, g, n); // search = compress_free(search) diff --git a/test/engine/engine_util_solve_test.cc b/test/engine/engine_util_solve_test.cc index d80abde1..ed224a68 100644 --- a/test/engine/engine_util_solve_test.cc +++ b/test/engine/engine_util_solve_test.cc @@ -19,6 +19,7 @@ #include #include #include +#include #include #include #include @@ -228,6 +229,58 @@ TEST_F(BoxQPTest, AsymmetricUpperIgnored) { EXPECT_MJTNUM_EQ(res[1], lower[1]); } +// verify mju_boxQP reads only the lower triangle of H by poisoning the upper +// triangle with NaN and comparing to a clean symmetric solve (see issue #3275) +TEST_F(BoxQPTest, UpperTrianglePoisoned) { + int n = 30; + const mjtNum nan = std::numeric_limits::quiet_NaN(); + + // allocate on heap + mjtNum *H, *g, *lower, *upper; // inputs + mjtNum *res, *R; // outputs + int* index; // outputs + mju_boxQPmalloc(&res, &R, &index, &H, &g, n, &lower, &upper); + + // generate a symmetric SPD Hessian and bounded QP problem + randomBoxQP(n, H, g, lower, upper, /*seed=*/1); + + // solve with symmetric H to get the reference result + mju_zero(res, n); + int nfree_ref = mju_boxQP(res, R, index, H, g, n, lower, upper); + ASSERT_GT(nfree_ref, -1); + + // save reference + std::vector res_ref(res, res + n); + std::vector index_ref(index, index + n); + + // poison the strict upper triangle of H with NaN + for (int i=0; i < n; i++) { + for (int j=i+1; j < n; j++) { + H[n*i+j] = nan; + } + } + + // solve again; result must match because only lower triangle should be read + mju_zero(res, n); + int nfree_poisoned = mju_boxQP(res, R, index, H, g, n, lower, upper); + + EXPECT_EQ(nfree_poisoned, nfree_ref); + for (int i=0; i < n; i++) { + EXPECT_EQ(res[i], res_ref[i]) << "mismatch at index " << i; + } + for (int i=0; i < nfree_ref; i++) { + EXPECT_EQ(index[i], index_ref[i]) << "index mismatch at " << i; + } + + mju_free(res); + mju_free(R); + mju_free(index); + mju_free(H); + mju_free(g); + mju_free(lower); + mju_free(upper); +} + // test mju_boxQP on a single random bounded QP TEST_F(BoxQPTest, BoundedQP) { int n = 50; // problem size