Fix dropped derivatives in Eigen Cholesky solves of AutoDiff matrices

[bigboateng]

Towards #17037.

The problem

#17037 tracked wrong derivatives from dynamic-size llt() / ldlt() solves of AutoDiff matrices: Eigen's triangular solvers skip the divide-and-propagate step for any right-hand-side component that compares equal to zero, and for autodiff types == 0 sees only the value. A component whose value is exactly zero but whose derivatives are nonzero (e.g., velocities when linearizing about a fixed point) therefore has its derivatives silently dropped.

Two things have changed since that discussion, and their combination reintroduces the bug on current master:

1. Eigen's fix (libeigen/eigen!976) shipped in Eigen 5.0.0+, but was not backported to the 3.4 branch (3.4.1 still reproduces the original problem). The fix works through a scalar customization point, Eigen::internal::is_identically_zero, which Eigen specializes only for its own AutoDiffScalar.
2. Drake has since replaced Eigen::AutoDiffScalar with drake::ad::AutoDiff, which does not specialize that customization point — so the value-only fallback applies and derivatives are dropped again, even with the Eigen 5.0.1 that Drake now pins.

Minimal repro on current master (Eigen 5.0.1): solving diag(4, 9) x = b with b0 = 0 (∂b0/∂b0 = 1), b1 = 5 (∂b1/∂b1 = 1) via dynamic-size llt().solve(b) yields ∂x0/∂b0 == 1 instead of 0.25 — the derivative is never divided by the Cholesky diagonal.

All current in-tree call sites happen to be unaffected (they use math::LinearSolver, double matrices, or fixed-size solves, which take Eigen's unrolled triangular path that has no zero-skip). Dynamic-size solves in downstream user code are exposed.

The fix

Specialize Eigen::internal::is_identically_zero for AutoDiff so the check is derivative-aware (true only when the value and all partials are zero), mirroring the specialization Eigen ships for its own AutoDiffScalar.

The specialization is guarded by EIGEN_VERSION_AT_LEAST(5, 0, 0): Eigen 3.4.x has no such customization point, so builds against a system Eigen 3.4.x remain subject to #17037. The new regression test skips on 3.4.x with a pointer to the issue, and the long-standing guidance to prefer drake::math::LinearSolver over raw Eigen factorizations for AutoDiff still applies there.

Performance: for components with nonzero values the check costs a single value comparison, as before. The derivative scan only runs when a component's value is exactly zero — the case that previously computed wrong results.

Testing

New regression test common/ad/test/eigen_cholesky_test.cc:

• dynamic llt() solve — exercises the non-unit-diagonal division skip;
• dynamic ldlt() solve — exercises the unit-diagonal propagation skip;
• fixed-size llt() solve — control case documenting the unaffected boundary.

Both dynamic cases fail without the specialization and pass with it. Expected gradients are the analytic M⁻¹ (constant matrix, so ∂x/∂b = M⁻¹ exactly).

---

This change is 

[jwnimmer-tri]

@drake-jenkins-bot test this please

[jwnimmer-tri]

Thanks for sending this fix!

Just a few small notes. I opened a PR with my suggested fixes.

bigboateng#1

@jwnimmer-tri reviewed 3 files and all commit messages, and made 5 comments.
Reviewable status: 4 unresolved discussions, LGTM missing from assignee jwnimmer-tri(platform), needs at least two assigned reviewers (waiting on bigboateng).

---

common/ad/internal/eigen_specializations.h line 73 at r1 (raw file):

// Eigen 3.4.x remain subject to drake#17037 and should use
// drake::math::LinearSolver instead of calling Eigen's factorizations directly
// on AutoDiff matrices.

nit This comment explanation is good for a commit message or regression test, but is too much irrelevant detail for this header. Please remove it.

If we had a citation for the thing we're specializing (is_identically_zero) we could hyperlink to that in a comment (like we did with NumTraits above), but I don't immediately find anything by searching.

---

common/ad/test/eigen_cholesky_test.cc line 26 at r1 (raw file):

// Returns the gradient of v with respect to the two independent variables.
MatrixXd ExtractGradient(const VectorX<AutoDiff>& v) {

nit We can use the Drake function for this. It's OK for a unit test in //common to depend on the //math package, especially when it's a small, standalone regression test like this one.

My branch contains the suggested fixup.

---

common/ad/test/eigen_cholesky_test.cc line 49 at r1 (raw file):

// diagonal. x[0] has value zero with nonzero derivatives.
GTEST_TEST(EigenCholeskyTest, LltDynamicSolve) {
#if !EIGEN_VERSION_AT_LEAST(5, 0, 0)

nit Even if the answer is wrong, there is still value in running the test to completion. It can at least check that the test doesn't crash, and the value() part should at least still be correct.

My branch contains the suggested fixup.

---

common/ad/test/eigen_cholesky_test.cc line 97 at r1 (raw file):

// on all supported Eigen versions).
GTEST_TEST(EigenCholeskyTest, LltFixedSizeSolve) {
  Eigen::Matrix<AutoDiff, 2, 2> M;

nit Drake has compact alises for things like this (Matrix2<T>).

My contains the suggested fixup.

[bigboateng]

Thanks for sending this fix!

Just a few small notes. I opened a PR with my suggested fixes.

bigboateng#1

@jwnimmer-tri reviewed 3 files and all commit messages, and made 5 comments.
Reviewable status: 4 unresolved discussions, LGTM missing from assignee jwnimmer-tri(platform), needs at least two assigned reviewers (waiting on bigboateng).

common/ad/internal/eigen_specializations.h line 73 at r1 (raw file):

// Eigen 3.4.x remain subject to drake#17037 and should use
// drake::math::LinearSolver instead of calling Eigen's factorizations directly
// on AutoDiff matrices.

nit This comment explanation is good for a commit message or regression test, but is too much irrelevant detail for this header. Please remove it.

If we had a citation for the thing we're specializing (is_identically_zero) we could hyperlink to that in a comment (like we did with NumTraits above), but I don't immediately find anything by searching.

common/ad/test/eigen_cholesky_test.cc line 26 at r1 (raw file):

// Returns the gradient of v with respect to the two independent variables.
MatrixXd ExtractGradient(const VectorX<AutoDiff>& v) {

nit We can use the Drake function for this. It's OK for a unit test in //common to depend on the //math package, especially when it's a small, standalone regression test like this one.

My branch contains the suggested fixup.

common/ad/test/eigen_cholesky_test.cc line 49 at r1 (raw file):

// diagonal. x[0] has value zero with nonzero derivatives.
GTEST_TEST(EigenCholeskyTest, LltDynamicSolve) {
#if !EIGEN_VERSION_AT_LEAST(5, 0, 0)

nit Even if the answer is wrong, there is still value in running the test to completion. It can at least check that the test doesn't crash, and the value() part should at least still be correct.

My branch contains the suggested fixup.

common/ad/test/eigen_cholesky_test.cc line 97 at r1 (raw file):

// on all supported Eigen versions).
GTEST_TEST(EigenCholeskyTest, LltFixedSizeSolve) {
  Eigen::Matrix<AutoDiff, 2, 2> M;

nit Drake has compact alises for things like this (Matrix2<T>).

My contains the suggested fixup.

Thanks! I’ve applied all four suggested changes by merging bigboateng#1. PR #24796 now points to 068f3a7. The updated head is awaiting CI.

[jwnimmer-tri]

feature review.

+a:@rpoyner-tri for platform review per schedule on Monday, please.

@jwnimmer-tri reviewed 3 files and all commit messages, made 1 comment, and resolved 4 discussions.
Reviewable status: LGTM missing from assignee rpoyner-tri(platform), commits need curation (https://drake.mit.edu/reviewable.html#curated-commits) (waiting on rpoyner-tri).