# Am I Leaking? [Eigen autodiff]

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

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>`?

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.

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?

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.
``````

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.

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

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.

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
``````

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.

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`

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;
}
``````

`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.