$\newcommand{\W}[1]{ \; #1 \; } \newcommand{\R}[1]{ {\rm #1} } \newcommand{\B}[1]{ {\bf #1} } \newcommand{\D}[2]{ \frac{\partial #1}{\partial #2} } \newcommand{\DD}[3]{ \frac{\partial^2 #1}{\partial #2 \partial #3} } \newcommand{\Dpow}[2]{ \frac{\partial^{#1}}{\partial {#2}^{#1}} } \newcommand{\dpow}[2]{ \frac{ {\rm d}^{#1}}{{\rm d}\, {#2}^{#1}} }$
# include <cppad/example/cppad_eigen.hpp> # include <cppad/speed/det_by_minor.hpp> # include <Eigen/Dense> bool eigen_det(void) { bool ok = true; using CppAD::AD; using CppAD::NearEqual; using Eigen::Matrix; using Eigen::Dynamic; // typedef Matrix< double , Dynamic, Dynamic > matrix; typedef Matrix< AD<double> , Dynamic, Dynamic > a_matrix; // typedef Matrix< double , Dynamic , 1> vector; typedef Matrix< AD<double> , Dynamic , 1> a_vector; // some temporary indices size_t i, j; // domain and range space vectors size_t size = 3, n = size * size, m = 1; a_vector a_x(n), a_y(m); vector x(n); // set and declare independent variables and start tape recording for(i = 0; i < size; i++) { for(j = 0; j < size; j++) { // lower triangular matrix a_x[i * size + j] = x[i * size + j] = 0.0; if( j <= i ) a_x[i * size + j] = x[i * size + j] = double(1 + i + j); } } CppAD::Independent(a_x); // copy independent variable vector to a matrix a_matrix a_X(size, size); matrix X(size, size); for(i = 0; i < size; i++) { for(j = 0; j < size; j++) { X(i, j) = x[i * size + j]; // If we used a_X(i, j) = X(i, j), a_X would not depend on a_x. a_X(i, j) = a_x[i * size + j]; } } // Compute the log of determinant of X a_y[0] = log( a_X.determinant() ); // create f: x -> y and stop tape recording CppAD::ADFun<double> f(a_x, a_y); // check function value double eps = 100. * CppAD::numeric_limits<double>::epsilon(); CppAD::det_by_minor<double> det(size); ok &= NearEqual(Value(a_y[0]) , log(det(x)), eps, eps); // compute the derivative of y w.r.t x using CppAD vector jac = f.Jacobian(x); // check the derivative using the formula // d/dX log(det(X)) = transpose( inv(X) ) matrix inv_X = X.inverse(); for(i = 0; i < size; i++) { for(j = 0; j < size; j++) ok &= NearEqual(jac[i * size + j], inv_X(j, i), eps, eps); } return ok; }