Solving linear systems with AutoDiffScalar produces wrong results with LLT and LDLT solver

Summary

Environment

• Operating System : Linux
• Architecture : x64
• Eigen Version : 3.4.0
• Compiler Version : Gcc8.4
• Compile Flags :
• Vector Extension :

Minimal Example

//show your code here
#include <Eigen/Core>
#include <Eigen/Dense>
#include <iostream>
#include <unsupported/Eigen/AutoDiff>

using namespace Eigen;
template <typename LinearSolverType>
void CheckLinearSolver(
    const LinearSolverType& linear_solver,
    const Eigen::Matrix<Eigen::AutoDiffScalar<Eigen::VectorXd>, 2, 2>& A_ad) {
  Eigen::Matrix<Eigen::AutoDiffScalar<Eigen::VectorXd>, 2, 1> b_ad;
  b_ad(0).value() = 3;
  b_ad(1).value() = 5;
  b_ad(0).derivatives() = Eigen::Vector3d(0, 0, 0);
  b_ad(1).derivatives() = Eigen::Vector3d(0, 0, 0);

  const auto x = linear_solver.solve(b_ad);
  const Matrix<Eigen::AutoDiffScalar<Eigen::VectorXd>, 2, 1> Ax = A_ad * x;
  std::cout << "A*x value: " << Ax(0).value() << ", " << Ax(1).value() << "\n";
  std::cout << "b value: " << b_ad(0).value() << ", " << b_ad(1).value()
            << "\n";
  // The derivative of A * x should be 0, since the right-hand side b_ad
  // contains no gradient.
  std::cout << "A*x derivatives:\n"
            << Ax(0).derivatives().transpose() << "\n"
            << Ax(1).derivatives().transpose() << "\n";
  std::cout << "b derivatives:\n"
            << b_ad(0).derivatives().transpose() << "\n"
            << b_ad(1).derivatives().transpose() << "\n";

  std::cout << "--------------------\n";
}

int main() {
  using TypeA = Eigen::Matrix<Eigen::AutoDiffScalar<Eigen::VectorXd>, 2, 2>;
  TypeA A_ad;
  A_ad(0, 0).value() = 1;
  A_ad(0, 1).value() = 3;
  A_ad(1, 0).value() = 3;
  A_ad(1, 1).value() = 10;
  A_ad(0, 0).derivatives() = Eigen::Vector3d(1, 2, 3);
  A_ad(0, 1).derivatives() = Eigen::Vector3d(4, 5, 6);
  A_ad(1, 0).derivatives() = Eigen::Vector3d(7, 8, 9);
  A_ad(1, 1).derivatives() = Eigen::Vector3d(10, 11, 12);

  const Eigen::LLT<TypeA> llt_solver(A_ad);
  std::cout << "LLT solver\n";
  CheckLinearSolver(llt_solver, A_ad);

  const Eigen::LDLT<TypeA> ldlt_solver(A_ad);
  std::cout << "LDLT solver\n";
  CheckLinearSolver(ldlt_solver, A_ad);

  const Eigen::PartialPivLU<TypeA> partial_piv_lu_solver(A_ad);
  std::cout << "PartialPivLU solver\n";
  CheckLinearSolver(partial_piv_lu_solver, A_ad);

  const Eigen::FullPivLU<TypeA> full_piv_lu_solver(A_ad);
  std::cout << "FullPivLU solver\n";
  CheckLinearSolver(full_piv_lu_solver, A_ad);

  const Eigen::HouseholderQR<TypeA> householder_qr_solver(A_ad);
  std::cout << "HouseholderQR solver\n";
  CheckLinearSolver(householder_qr_solver, A_ad);

  const Eigen::ColPivHouseholderQR<TypeA> col_piv_householder_qr_solver(A_ad);
  std::cout << "ColPivHouseholderQR solver\n";
  CheckLinearSolver(col_piv_householder_qr_solver, A_ad);

  const Eigen::FullPivHouseholderQR<TypeA> full_piv_householder_qr_solver(A_ad);
  std::cout << "FullPivHouseholderQR solver\n";
  CheckLinearSolver(full_piv_householder_qr_solver, A_ad);

  const Eigen::CompleteOrthogonalDecomposition<TypeA>
      complete_orthogonal_decomposition_solver(A_ad);
  std::cout << "CompleteOrthogonalDecomposition solver\n";
  CheckLinearSolver(complete_orthogonal_decomposition_solver, A_ad);

  // BDCSVD and JacobiSVD don't support autodiffing.
  return 0;
}

Steps to reproduce

1. compile the code as gcc -Ipath_to_eigen test.cc
2. run the generated binary

What is the current bug behavior?

When using LLT or LDLT solver, the computed solution doesn't have the right gradient, that A*x have non-zero derivatives, while the right-hand side b has zero derivatives.

What is the expected correct behavior?

Each solver should compute a result, such that A*x has zero derivatives.

Relevant logs

Anything else that might help

• Have a plan to fix this issue.