Am I Leaking? [Eigen autodiff]

#1

I have been trying for the last week to calculate the logarithm of the absolute value of the determinant of the Jacobian matrix of a transformation from one square matrix to another and am attempting to verify my math. Has anyone used the unsupported AutoDiff module in Eigen with stan::math::var as the underlying type? Here is a simple example file testJ.cpp (1.1 KB)
yields the correct answer for a matrix inverse transformation, but I don’t know if I can hook something like this into a Stan program without having to do some memory management stuff to prevent memory leaks. Download the attachment and then do in R

library(Rcpp)
sourceCpp("testJ.cpp")
K <- 3L
A <- matrix(rnorm(K ^ 2), K, K)
test_Jacobian(A) # yields same value as the next line
-2 * K * log(abs(det(A)))
1 Like
#2

Did you mean using it as an externally defined function? What you’d want to do is execute nested autodiff which will allocate all thos vars on the top of the stack which can then be cleared off later.

Here’s the code in testJ.cpp for comparison:

double test_Jacobian(Eigen::MatrixXd A) {
  stan::math::matrix_v A_var = stan::math::to_var(A);
  typedef Eigen::AutoDiffScalar<stan::math::vector_v> ADScalar;
  typedef Eigen::Matrix<ADScalar,Eigen::Dynamic,Eigen::Dynamic> ADMatrix;
  ADMatrix A_(A.rows(), A.cols());
  for (int i = 0; i < A.size(); i++) A_(i).value() = A_var(i);
  const int derivative_num = A.size();
  int derivative_idx = 0;
  for (int i = 0; i < A.size(); i++) {
    A_(i).derivatives() = 
      Eigen::VectorXd::Unit(derivative_num, derivative_idx);
    derivative_idx++;
  }
  ADMatrix A_inv = A_.inverse();
  stan::math::matrix_v J(derivative_num, derivative_num);
  int pos = 0;
  for (int i = 0; i < A.rows(); i++) for (int j = 0; j < A.cols(); j++) {
    J.col(pos++) = A_inv(i, j).derivatives();
  }
  stan::math::var log_abs_det_J = stan::math::log_determinant(J);
  return log_abs_det_J.val();
}

Is ADMatrix defined to be Eigen::Matrix<var, -1, -1>?

#3

The ADMatrix typedef is closer to a matrix of fvar<var>. But the whole idea only works for a couple of examples. Although the documentation
http://eigen.tuxfamily.org/dox/unsupported/classEigen_1_1AutoDiffScalar.html
says functions like std::abs, std::sqrt, std::pow, std::exp, std::log, std::sin, std::cos, etc. are defined, if any of those get called, then it does not compile the way we have defined these in rev.

#4

If they’re defined and we’re including the source, they should get picked up with argument-dependent lookup. That’ll certainly happen for any stan::math::var if rev/mat.hpp is included.

What’s the compiler error you get?

#5

If, for example, you change the line that says ADMatrix A_inv = A_.inverse(); to ADMatrix U = A_.jacobiSvd(Eigen::ComputeFullU).matrixU(); (and change A_inv to U four lines down) you get calls to sqrt and hence a compiler error:

In file included from /home/ben/r-devel/library/StanHeaders/include/stan/math.hpp:4:
In file included from /home/ben/r-devel/library/StanHeaders/include/stan/math/rev/mat.hpp:11:
In file included from /home/ben/r-devel/library/StanHeaders/include/stan/math/rev/mat/vectorize/apply_scalar_unary.hpp:4:
/home/ben/r-devel/library/StanHeaders/include/stan/math/prim/mat/vectorize/apply_scalar_unary.hpp:39:37: error: implicit instantiation of undefined template
      'Eigen::internal::traits<Eigen::AutoDiffScalar<Eigen::CwiseBinaryOp<Eigen::internal::scalar_sum_op<stan::math::var, stan::math::var>, const Eigen::Matrix<stan::math::var, -1, 1, 0,
      -1, 1>, const Eigen::Matrix<stan::math::var, -1, 1, 0, -1, 1> > > >'
  typedef typename Eigen::internal::traits<T>::Scalar scalar_t;
                                    ^
/home/ben/r-devel/library/StanHeaders/include/stan/math/prim/mat/fun/sqrt.hpp:31:17: note: in instantiation of template class 'stan::math::apply_scalar_unary<stan::math::sqrt_fun,
      Eigen::AutoDiffScalar<Eigen::CwiseBinaryOp<Eigen::internal::scalar_sum_op<stan::math::var, stan::math::var>, const Eigen::Matrix<stan::math::var, -1, 1, 0, -1, 1>, const
      Eigen::Matrix<stan::math::var, -1, 1, 0, -1, 1> > > >' requested here
inline typename apply_scalar_unary<sqrt_fun, T>::return_t sqrt(const T& x) {
                ^
/home/ben/r-devel/library/RcppEigen/include/Eigen/src/Householder/Householder.h:88:12: note: while substituting deduced template arguments into function template 'sqrt' [with T =
      Eigen::AutoDiffScalar<Eigen::CwiseBinaryOp<Eigen::internal::scalar_sum_op<stan::math::var, stan::math::var>, const Eigen::Matrix<stan::math::var, -1, 1, 0, -1, 1>, const
      Eigen::Matrix<stan::math::var, -1, 1, 0, -1, 1> > >]
    beta = sqrt(numext::abs2(c0) + tailSqNorm);
           ^
/home/ben/r-devel/library/RcppEigen/include/Eigen/src/Householder/Householder.h:45:3: note: in instantiation of function template specialization
      'Eigen::MatrixBase<Eigen::Block<Eigen::Block<Eigen::Matrix<Eigen::AutoDiffScalar<Eigen::Matrix<stan::math::var, -1, 1, 0, -1, 1> >, -1, -1, 0, -1, -1>, -1, 1, true>, -1, 1, false>
      >::makeHouseholder<Eigen::VectorBlock<Eigen::Block<Eigen::Block<Eigen::Matrix<Eigen::AutoDiffScalar<Eigen::Matrix<stan::math::var, -1, 1, 0, -1, 1> >, -1, -1, 0, -1, -1>, -1, 1,
      true>, -1, 1, false>, -1> >' requested here
  makeHouseholder(essentialPart, tau, beta);
  ^
/home/ben/r-devel/library/RcppEigen/include/Eigen/src/QR/ColPivHouseholderQR.h:538:30: note: in instantiation of member function
      'Eigen::MatrixBase<Eigen::Block<Eigen::Block<Eigen::Matrix<Eigen::AutoDiffScalar<Eigen::Matrix<stan::math::var, -1, 1, 0, -1, 1> >, -1, -1, 0, -1, -1>, -1, 1, true>, -1, 1, false>
      >::makeHouseholderInPlace' requested here
    m_qr.col(k).tail(rows-k).makeHouseholderInPlace(m_hCoeffs.coeffRef(k), beta);
                             ^
/home/ben/r-devel/library/RcppEigen/include/Eigen/src/QR/ColPivHouseholderQR.h:475:3: note: in instantiation of member function
      'Eigen::ColPivHouseholderQR<Eigen::Matrix<Eigen::AutoDiffScalar<Eigen::Matrix<stan::math::var, -1, 1, 0, -1, 1> >, -1, -1, 0, -1, -1> >::computeInPlace' requested here
  computeInPlace();
  ^
/home/ben/r-devel/library/RcppEigen/include/Eigen/src/SVD/JacobiSVD.h:230:12: note: in instantiation of function template specialization
      'Eigen::ColPivHouseholderQR<Eigen::Matrix<Eigen::AutoDiffScalar<Eigen::Matrix<stan::math::var, -1, 1, 0, -1, 1> >, -1, -1, 0, -1, -1>
      >::compute<Eigen::Matrix<Eigen::AutoDiffScalar<Eigen::Matrix<stan::math::var, -1, 1, 0, -1, 1> >, -1, -1, 0, -1, -1> >' requested here
      m_qr.compute(m_adjoint);
           ^
/home/ben/r-devel/library/RcppEigen/include/Eigen/src/SVD/JacobiSVD.h:684:27: note: in instantiation of member function
      'Eigen::internal::qr_preconditioner_impl<Eigen::Matrix<Eigen::AutoDiffScalar<Eigen::Matrix<stan::math::var, -1, 1, 0, -1, 1> >, -1, -1, 0, -1, -1>, 2, 0, true>::run' requested here
    m_qr_precond_morecols.run(*this, m_scaledMatrix);
                          ^
/home/ben/r-devel/library/RcppEigen/include/Eigen/src/SVD/JacobiSVD.h:548:7: note: in instantiation of member function
      'Eigen::JacobiSVD<Eigen::Matrix<Eigen::AutoDiffScalar<Eigen::Matrix<stan::math::var, -1, 1, 0, -1, 1> >, -1, -1, 0, -1, -1>, 2>::compute' requested here
      compute(matrix, computationOptions);
      ^
/home/ben/r-devel/library/RcppEigen/include/Eigen/src/SVD/JacobiSVD.h:799:10: note: in instantiation of member function
      'Eigen::JacobiSVD<Eigen::Matrix<Eigen::AutoDiffScalar<Eigen::Matrix<stan::math::var, -1, 1, 0, -1, 1> >, -1, -1, 0, -1, -1>, 2>::JacobiSVD' requested here
  return JacobiSVD<PlainObject>(*this, computationOptions);
         ^
testJ2.cpp:26:19: note: in instantiation of member function 'Eigen::MatrixBase<Eigen::Matrix<Eigen::AutoDiffScalar<Eigen::Matrix<stan::math::var, -1, 1, 0, -1, 1> >, -1, -1, 0, -1, -1>
      >::jacobiSvd' requested here
  ADMatrix U = A_.jacobiSvd(Eigen::ComputeFullU).matrixU();
                  ^
/home/ben/r-devel/library/RcppEigen/include/Eigen/src/Core/util/ForwardDeclarations.h:17:29: note: template is declared here
template<typename T> struct traits;
                            ^
In file included from testJ2.cpp:4:
In file included from /home/ben/r-devel/library/StanHeaders/include/stan/math.hpp:4:
In file included from /home/ben/r-devel/library/StanHeaders/include/stan/math/rev/mat.hpp:11:
In file included from /home/ben/r-devel/library/StanHeaders/include/stan/math/rev/mat/vectorize/apply_scalar_unary.hpp:4:
/home/ben/r-devel/library/StanHeaders/include/stan/math/prim/mat/vectorize/apply_scalar_unary.hpp:45:38: error: no member named 'RowsAtCompileTime' in
      'Eigen::AutoDiffScalar<Eigen::CwiseBinaryOp<Eigen::internal::scalar_sum_op<stan::math::var, stan::math::var>, const Eigen::Matrix<stan::math::var, -1, 1, 0, -1, 1>, const
      Eigen::Matrix<stan::math::var, -1, 1, 0, -1, 1> > >'
  typedef Eigen::Matrix<scalar_t, T::RowsAtCompileTime, T::ColsAtCompileTime>
                                  ~~~^
/home/ben/r-devel/library/StanHeaders/include/stan/math/prim/mat/vectorize/apply_scalar_unary.hpp:39:37: error: implicit instantiation of undefined template
      'Eigen::internal::traits<Eigen::AutoDiffScalar<Eigen::Matrix<stan::math::var, -1, 1, 0, -1, 1> > >'
  typedef typename Eigen::internal::traits<T>::Scalar scalar_t;
                                    ^
/home/ben/r-devel/library/StanHeaders/include/stan/math/prim/mat/fun/sqrt.hpp:31:17: note: in instantiation of template class 'stan::math::apply_scalar_unary<stan::math::sqrt_fun,
      Eigen::AutoDiffScalar<Eigen::Matrix<stan::math::var, -1, 1, 0, -1, 1> > >' requested here
inline typename apply_scalar_unary<sqrt_fun, T>::return_t sqrt(const T& x) {
                ^
/home/ben/r-devel/library/RcppEigen/include/Eigen/src/Core/MathFunctions.h:1020:10: note: while substituting deduced template arguments into function template 'sqrt' [with T =
      Eigen::AutoDiffScalar<Eigen::Matrix<stan::math::var, -1, 1, 0, -1, 1> >]
  return sqrt(x);
         ^
/home/ben/r-devel/library/RcppEigen/include/Eigen/src/QR/ColPivHouseholderQR.h:568:52: note: in instantiation of function template specialization
      'Eigen::numext::sqrt<Eigen::AutoDiffScalar<Eigen::Matrix<stan::math::var, -1, 1, 0, -1, 1> > >' requested here
          m_colNormsUpdated.coeffRef(j) *= numext::sqrt(temp);
                                                   ^
/home/ben/r-devel/library/RcppEigen/include/Eigen/src/QR/ColPivHouseholderQR.h:475:3: note: in instantiation of member function
      'Eigen::ColPivHouseholderQR<Eigen::Matrix<Eigen::AutoDiffScalar<Eigen::Matrix<stan::math::var, -1, 1, 0, -1, 1> >, -1, -1, 0, -1, -1> >::computeInPlace' requested here
  computeInPlace();
  ^
/home/ben/r-devel/library/RcppEigen/include/Eigen/src/SVD/JacobiSVD.h:230:12: note: in instantiation of function template specialization
      'Eigen::ColPivHouseholderQR<Eigen::Matrix<Eigen::AutoDiffScalar<Eigen::Matrix<stan::math::var, -1, 1, 0, -1, 1> >, -1, -1, 0, -1, -1>
      >::compute<Eigen::Matrix<Eigen::AutoDiffScalar<Eigen::Matrix<stan::math::var, -1, 1, 0, -1, 1> >, -1, -1, 0, -1, -1> >' requested here
      m_qr.compute(m_adjoint);
           ^
/home/ben/r-devel/library/RcppEigen/include/Eigen/src/SVD/JacobiSVD.h:684:27: note: in instantiation of member function
      'Eigen::internal::qr_preconditioner_impl<Eigen::Matrix<Eigen::AutoDiffScalar<Eigen::Matrix<stan::math::var, -1, 1, 0, -1, 1> >, -1, -1, 0, -1, -1>, 2, 0, true>::run' requested here
    m_qr_precond_morecols.run(*this, m_scaledMatrix);
                          ^
/home/ben/r-devel/library/RcppEigen/include/Eigen/src/SVD/JacobiSVD.h:548:7: note: in instantiation of member function
      'Eigen::JacobiSVD<Eigen::Matrix<Eigen::AutoDiffScalar<Eigen::Matrix<stan::math::var, -1, 1, 0, -1, 1> >, -1, -1, 0, -1, -1>, 2>::compute' requested here
      compute(matrix, computationOptions);
      ^
/home/ben/r-devel/library/RcppEigen/include/Eigen/src/SVD/JacobiSVD.h:799:10: note: in instantiation of member function
      'Eigen::JacobiSVD<Eigen::Matrix<Eigen::AutoDiffScalar<Eigen::Matrix<stan::math::var, -1, 1, 0, -1, 1> >, -1, -1, 0, -1, -1>, 2>::JacobiSVD' requested here
  return JacobiSVD<PlainObject>(*this, computationOptions);
         ^
testJ2.cpp:26:19: note: in instantiation of member function 'Eigen::MatrixBase<Eigen::Matrix<Eigen::AutoDiffScalar<Eigen::Matrix<stan::math::var, -1, 1, 0, -1, 1> >, -1, -1, 0, -1, -1>
      >::jacobiSvd' requested here
  ADMatrix U = A_.jacobiSvd(Eigen::ComputeFullU).matrixU();
                  ^
/home/ben/r-devel/library/RcppEigen/include/Eigen/src/Core/util/ForwardDeclarations.h:17:29: note: template is declared here
template<typename T> struct traits;
                            ^
In file included from testJ2.cpp:4:
In file included from /home/ben/r-devel/library/StanHeaders/include/stan/math.hpp:4:
In file included from /home/ben/r-devel/library/StanHeaders/include/stan/math/rev/mat.hpp:11:
In file included from /home/ben/r-devel/library/StanHeaders/include/stan/math/rev/mat/vectorize/apply_scalar_unary.hpp:4:
/home/ben/r-devel/library/StanHeaders/include/stan/math/prim/mat/vectorize/apply_scalar_unary.hpp:45:38: error: no member named 'RowsAtCompileTime' in
      'Eigen::AutoDiffScalar<Eigen::Matrix<stan::math::var, -1, 1, 0, -1, 1> >'
  typedef Eigen::Matrix<scalar_t, T::RowsAtCompileTime, T::ColsAtCompileTime>
                                  ~~~^
1 warning and 4 errors generated.
#6

The root cause here seems to be this:

while substituting deduced template arguments into function template 'sqrt' 
[with T =
 Eigen::AutoDiffScalar
   <Eigen::CwiseBinaryOp
     <Eigen::internal::scalar_sum_op<stan::math::var,
                                     stan::math::var>, 
      const Eigen::Matrix<stan::math::var, -1, 1, 0, -1, 1>, 
      const Eigen::Matrix<stan::math::var, -1, 1, 0, -1, 1> > >]

beta = sqrt(numext::abs2(c0) + tailSqNorm);

Any idea why the AutoDiffScalar type here is taking a CwiseBinaryOp template parameter? And it looks like that binary op’s defined over vectors.

#7

I think that is because the forward mode derivatives are being held in a vector (of var in this case).

#8

It’d be interesting to see if it’s our var messing it up or if this would work with double in place of stan::math::var. It doesn’t look like a Stan problem. It looks like sqrt isn’t defined for this AutodiffScalar type.

#9

This double-version compiles (and runs and does not crash)

// [[Rcpp::depends(BH)]]
// [[Rcpp::depends(RcppEigen)]]
// [[Rcpp::depends(StanHeaders)]]
#include <stan/math.hpp>
#include <Rcpp.h>
#include <RcppEigen.h>
#include <unsupported/Eigen/AutoDiff>
#include <vector>

// [[Rcpp::export]]
Eigen::MatrixXd test_Jacobian(Eigen::MatrixXd A) {
  typedef Eigen::AutoDiffScalar<stan::math::vector_d> ADScalar;
  typedef Eigen::Matrix<ADScalar,Eigen::Dynamic,Eigen::Dynamic> ADMatrix;
  ADMatrix A_(A.rows(), A.cols());
  const int derivative_num = A.size();
  int derivative_idx = 0;
  for (int i = 0; i < A.size(); i++) {
    A_(i).value() = A(i);
    A_(i).derivatives() =
      Eigen::VectorXd::Unit(derivative_num, derivative_idx);
    derivative_idx++;
  }
  ADMatrix U = A_.fullPivHouseholderQr().matrixQ();
  stan::math::matrix_d J(derivative_num, derivative_num);
  int pos = 0;
  for (int i = 0; i < A.rows(); i++) for (int j = 0; j < A.cols(); j++) {
    J.col(pos++) = U(i, j).derivatives();
  }
  return J;
}

while the var version does not compile

// [[Rcpp::depends(BH)]]
// [[Rcpp::depends(RcppEigen)]]
// [[Rcpp::depends(StanHeaders)]]
#include <stan/math.hpp>
#include <Rcpp.h>
#include <RcppEigen.h>
#include <unsupported/Eigen/AutoDiff>
#include <vector>

// [[Rcpp::export]]
std::vector<double> test_Jacobian(Eigen::MatrixXd A) {
  typedef Eigen::AutoDiffScalar<stan::math::vector_v> ADScalar;
  typedef Eigen::Matrix<ADScalar,Eigen::Dynamic,Eigen::Dynamic> ADMatrix;
  ADMatrix A_(A.rows(), A.cols());
  std::vector<stan::math::var> theta;
  std::vector<double> gradient;
  const int derivative_num = A.size();
  int derivative_idx = 0;
  for (int i = 0; i < A.size(); i++) {
    theta.push_back(stan::math::to_var(A(i)));
    A_(i).value() = theta.back();
    A_(i).value() = A(i);
    A_(i).derivatives() =
      Eigen::VectorXd::Unit(derivative_num, derivative_idx);
    derivative_idx++;
  }
  ADMatrix Q = A_.fullPivHouseholderQr().matrixQ();
  stan::math::matrix_v J(derivative_num, derivative_num);
  int pos = 0;
  for (int i = 0; i < A.rows(); i++) for (int j = 0; j < A.cols(); j++) {
    J.col(pos++) = Q(i, j).derivatives();
  }
  stan::math::var log_abs_det_J = stan::math::log_determinant(J);
  log_abs_det_J.grad(theta, gradient);
  return gradient;
}

due to the sqrt thing

In file included from /home/ben/R/x86_64-pc-linux-gnu-library/3.3/StanHeaders/include/stan/math.hpp:4:
In file included from /home/ben/R/x86_64-pc-linux-gnu-library/3.3/StanHeaders/include/stan/math/rev/mat.hpp:11:
In file included from /home/ben/R/x86_64-pc-linux-gnu-library/3.3/StanHeaders/include/stan/math/rev/mat/vectorize/apply_scalar_unary.hpp:4:
/home/ben/R/x86_64-pc-linux-gnu-library/3.3/StanHeaders/include/stan/math/prim/mat/vectorize/apply_scalar_unary.hpp:39:37: error: implicit instantiation of undefined template
      'Eigen::internal::traits<Eigen::AutoDiffScalar<Eigen::CwiseBinaryOp<Eigen::internal::scalar_sum_op<stan::math::var, stan::math::var>, const Eigen::Matrix<stan::math::var, -1, 1, 0,
      -1, 1>, const Eigen::Matrix<stan::math::var, -1, 1, 0, -1, 1> > > >'
  typedef typename Eigen::internal::traits<T>::Scalar scalar_t;
                                    ^
/home/ben/R/x86_64-pc-linux-gnu-library/3.3/StanHeaders/include/stan/math/prim/mat/fun/sqrt.hpp:31:17: note: in instantiation of template class
      'stan::math::apply_scalar_unary<stan::math::sqrt_fun, Eigen::AutoDiffScalar<Eigen::CwiseBinaryOp<Eigen::internal::scalar_sum_op<stan::math::var, stan::math::var>, const
      Eigen::Matrix<stan::math::var, -1, 1, 0, -1, 1>, const Eigen::Matrix<stan::math::var, -1, 1, 0, -1, 1> > > >' requested here
inline typename apply_scalar_unary<sqrt_fun, T>::return_t sqrt(const T& x) {
                ^
/home/ben/R/x86_64-pc-linux-gnu-library/3.3/RcppEigen/include/Eigen/src/Householder/Householder.h:88:12: note: while substituting deduced template arguments into function template 'sqrt'
      [with T = Eigen::AutoDiffScalar<Eigen::CwiseBinaryOp<Eigen::internal::scalar_sum_op<stan::math::var, stan::math::var>, const Eigen::Matrix<stan::math::var, -1, 1, 0, -1, 1>, const
      Eigen::Matrix<stan::math::var, -1, 1, 0, -1, 1> > >]
    beta = sqrt(numext::abs2(c0) + tailSqNorm);
           ^
/home/ben/R/x86_64-pc-linux-gnu-library/3.3/RcppEigen/include/Eigen/src/Householder/Householder.h:45:3: note: in instantiation of function template specialization
      'Eigen::MatrixBase<Eigen::Block<Eigen::Block<Eigen::Matrix<Eigen::AutoDiffScalar<Eigen::Matrix<stan::math::var, -1, 1, 0, -1, 1> >, -1, -1, 0, -1, -1>, -1, 1, true>, -1, 1, false>
      >::makeHouseholder<Eigen::VectorBlock<Eigen::Block<Eigen::Block<Eigen::Matrix<Eigen::AutoDiffScalar<Eigen::Matrix<stan::math::var, -1, 1, 0, -1, 1> >, -1, -1, 0, -1, -1>, -1, 1,
      true>, -1, 1, false>, -1> >' requested here
  makeHouseholder(essentialPart, tau, beta);
  ^
/home/ben/R/x86_64-pc-linux-gnu-library/3.3/RcppEigen/include/Eigen/src/QR/FullPivHouseholderQR.h:521:30: note: in instantiation of member function
      'Eigen::MatrixBase<Eigen::Block<Eigen::Block<Eigen::Matrix<Eigen::AutoDiffScalar<Eigen::Matrix<stan::math::var, -1, 1, 0, -1, 1> >, -1, -1, 0, -1, -1>, -1, 1, true>, -1, 1, false>
      >::makeHouseholderInPlace' requested here
    m_qr.col(k).tail(rows-k).makeHouseholderInPlace(m_hCoeffs.coeffRef(k), beta);
                             ^
/home/ben/R/x86_64-pc-linux-gnu-library/3.3/RcppEigen/include/Eigen/src/QR/FullPivHouseholderQR.h:452:3: note: in instantiation of member function
      'Eigen::FullPivHouseholderQR<Eigen::Matrix<Eigen::AutoDiffScalar<Eigen::Matrix<stan::math::var, -1, 1, 0, -1, 1> >, -1, -1, 0, -1, -1> >::computeInPlace' requested here
  computeInPlace();
  ^
/home/ben/R/x86_64-pc-linux-gnu-library/3.3/RcppEigen/include/Eigen/src/QR/FullPivHouseholderQR.h:136:7: note: in instantiation of function template specialization
      'Eigen::FullPivHouseholderQR<Eigen::Matrix<Eigen::AutoDiffScalar<Eigen::Matrix<stan::math::var, -1, 1, 0, -1, 1> >, -1, -1, 0, -1, -1>
      >::compute<Eigen::Matrix<Eigen::AutoDiffScalar<Eigen::Matrix<stan::math::var, -1, 1, 0, -1, 1> >, -1, -1, 0, -1, -1> >' requested here
      compute(matrix.derived());
      ^
/home/ben/R/x86_64-pc-linux-gnu-library/3.3/RcppEigen/include/Eigen/src/QR/FullPivHouseholderQR.h:671:10: note: in instantiation of function template specialization
      'Eigen::FullPivHouseholderQR<Eigen::Matrix<Eigen::AutoDiffScalar<Eigen::Matrix<stan::math::var, -1, 1, 0, -1, 1> >, -1, -1, 0, -1, -1>
      >::FullPivHouseholderQR<Eigen::Matrix<Eigen::AutoDiffScalar<Eigen::Matrix<stan::math::var, -1, 1, 0, -1, 1> >, -1, -1, 0, -1, -1> >' requested here
  return FullPivHouseholderQR<PlainObject>(eval());
         ^
testJ2.cpp:27:19: note: in instantiation of member function 'Eigen::MatrixBase<Eigen::Matrix<Eigen::AutoDiffScalar<Eigen::Matrix<stan::math::var, -1, 1, 0, -1, 1> >, -1, -1, 0, -1, -1>
      >::fullPivHouseholderQr' requested here
  ADMatrix Q = A_.fullPivHouseholderQr().matrixQ();
                  ^
/home/ben/R/x86_64-pc-linux-gnu-library/3.3/RcppEigen/include/Eigen/src/Core/util/ForwardDeclarations.h:17:29: note: template is declared here
template<typename T> struct traits;
                            ^
In file included from testJ2.cpp:4:
In file included from /home/ben/R/x86_64-pc-linux-gnu-library/3.3/StanHeaders/include/stan/math.hpp:4:
In file included from /home/ben/R/x86_64-pc-linux-gnu-library/3.3/StanHeaders/include/stan/math/rev/mat.hpp:11:
In file included from /home/ben/R/x86_64-pc-linux-gnu-library/3.3/StanHeaders/include/stan/math/rev/mat/vectorize/apply_scalar_unary.hpp:4:
/home/ben/R/x86_64-pc-linux-gnu-library/3.3/StanHeaders/include/stan/math/prim/mat/vectorize/apply_scalar_unary.hpp:45:38: error: no member named 'RowsAtCompileTime' in
      'Eigen::AutoDiffScalar<Eigen::CwiseBinaryOp<Eigen::internal::scalar_sum_op<stan::math::var, stan::math::var>, const Eigen::Matrix<stan::math::var, -1, 1, 0, -1, 1>, const
      Eigen::Matrix<stan::math::var, -1, 1, 0, -1, 1> > >'
  typedef Eigen::Matrix<scalar_t, T::RowsAtCompileTime, T::ColsAtCompileTime>
                                  ~~~^
2 errors generated.
make: *** [testJ2.o] Error 1
#10

It looks like you should be including our sqrt through stan/math.hpp.

What is Eigen::AutoDiffScalar and how doid that get involved here:

/home/ben/R/x86_64-pc-linux-gnu-library/3.3/StanHeaders/include/stan/math/prim/mat/fun/sqrt.hpp:31:17: note: in instantiation of template class
      'stan::math::apply_scalar_unary<stan::math::sqrt_fun, Eigen::AutoDiffScalar<Eigen::CwiseBinaryOp<Eigen::internal::scalar_sum_op<stan::math::var, stan::math::var>, const
      Eigen::Matrix<stan::math::var, -1, 1, 0, -1, 1>, const Eigen::Matrix<stan::math::var, -1, 1, 0, -1, 1> > > >' requested here
inline typename apply_scalar_unary<sqrt_fun, T>::return_t sqrt(const T& x) {
                ^

It looks like it’s trying to instantiate sqrt(T) with T = Eigen::AutoDiffScalar<Eigen::CwiseBinaryOp< ...>>.

To make this work, that’s going to have to be defined.

#11

That may be true, but if we have to add a define sqrt and every function like it, then it just becomes yet another partially complete forward mode autodiff implementation. And it isn’t something we can expect users to do when they define external C++ functions.

But I don’t think it should be the case that

  • Differentiating works if we typedef Eigen::AutoDiffScalar<stan::math::vector_d> ADScalar;
  • Differenatiating also works if we typedef Eigen::AutoDiffScalar<stan::math::vector_v> ADScalar; as long as we only use +, -, *, and /
  • But it does compile anymore if we use C++ math functions like sqrt
#12

Moreover, I can get it to compile if I don’t #include <stan/math.hpp> and instead #include <stan/math/rev/arr.hpp>

// [[Rcpp::depends(BH)]]
// [[Rcpp::depends(RcppEigen)]]
// [[Rcpp::depends(StanHeaders)]]
#include <Rcpp.h>
#include <RcppEigen.h>
#include <unsupported/Eigen/AutoDiff>
#include <stan/math/rev/arr.hpp>
#include <vector>

// [[Rcpp::export]]
std::vector<double> test_Jacobian(Eigen::MatrixXd A) {
  typedef Eigen::Matrix<stan::math::var, Eigen::Dynamic, 1> vector_v;
  typedef Eigen::Matrix<stan::math::var, Eigen::Dynamic, Eigen::Dynamic>
    matrix_v;
  typedef Eigen::AutoDiffScalar<vector_v> ADScalar;
  typedef Eigen::Matrix<ADScalar,Eigen::Dynamic,Eigen::Dynamic> ADMatrix;
  ADMatrix A_(A.rows(), A.cols());
  std::vector<stan::math::var> theta;
  std::vector<double> gradient;
  const int derivative_num = A.size();
  int derivative_idx = 0;
  for (int i = 0; i < A.size(); i++) {
    theta.push_back(stan::math::to_var(A(i)));
    A_(i).value() = theta.back();
    A_(i).value() = A(i);
    A_(i).derivatives() =
      Eigen::VectorXd::Unit(derivative_num, derivative_idx);
    derivative_idx++;
  }
  ADMatrix Q = A_.householderQr().householderQ();
  matrix_v J(derivative_num, derivative_num);
  int pos = 0;
  for (int i = 0; i < A.rows(); i++) for (int j = 0; j < A.cols(); j++) {
    J.col(pos++) = Q(i, j).derivatives();
  }
  stan::math::var y = J(0, 0);
  y.grad(theta, gradient);
  return gradient;
}
#13

rev/arr.hpp is the header that doesn’t include Eigen. So it sounds like maybe there’s some kind of interference of definitions somewhere.