LDLT produce wrong results with AutoDiffScalar

Summary

The derivatives of the result of LDLT.solve() via AutoDiffScalar is incorrect.

Environment

• Operating System : Linux
• Architecture : x64
• Eigen Version : 3.4.0
• Compiler Version : Gcc9.4.0/Clang10.0.0
• Compile Flags : None
• Vector Extension : None

Minimal Example

#include <iostream>
#include <Eigen/Dense>
#include <unsupported/Eigen/AutoDiff>

template <typename Scalar, int Rows>
using Vector = Eigen::Matrix<Scalar, Rows, 1>;

using Eigen::Vector2d;
using Eigen::VectorXd;
using Eigen::Dynamic;
using std::cout;
using std::endl;
using AutoDiff = Eigen::AutoDiffScalar<Eigen::VectorXd>;
using Eigen::Matrix;

int main()
{
  constexpr int kNumDerivatives = 3;
  constexpr int kNumVariables = 2;
  // Set up rhs b and derivatives.
  Vector<AutoDiff, kNumVariables> b(0,0);
  VectorXd db0(kNumDerivatives);
  db0 << 1,2,3;
  VectorXd db1(kNumDerivatives);
  db1 << 4,5,6;
  b[0].derivatives() = db0;
  b[1].derivatives() = db1;

  // Set up matrix A.
  Matrix<AutoDiff, kNumVariables, kNumVariables> A;
  A << 2, 1, 1, 2;

  // Solve A*x = b.
  Vector<AutoDiff, Dynamic> x = A.ldlt().solve(b);
  Vector<AutoDiff, Dynamic> Ax = A * x;

  cout << "d(A*x): " << endl;
  cout << Ax[0].derivatives().transpose() << endl;
  // Expect 4,5,6 printed here. But got 4.5, 6, 7.5.
  cout << Ax[1].derivatives().transpose() << endl;
  cout << endl;

  return 0;
}

Steps to reproduce

Save the above reproducing code to test.cc and g++ -I path/to/eigen/3.4.0 test.cc && ./a.out

What is the current bug behavior?

When solving A*x = b. The derivative of A*x is different from that of b.

What is the expected correct behavior?

The should match.

Anything else that might help

1. This didn't happen in Eigen 3.3.
2. Interestingly, if the line Vector<AutoDiff, Dynamic> x = A.ldlt().solve(b); is changed to Vector<AutoDiff, 2> x = A.ldlt().solve(b);, the problem goes away.
3. LLT suffers from the same issue.