1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
|
// SPDX-License-Identifier: EPL-2.0 OR GPL-2.0-or-later
// SPDX-FileCopyrightText: Bradley M. Bell <bradbell@seanet.com>
// SPDX-FileContributor: 2003-22 Bradley M. Bell
// ----------------------------------------------------------------------------
/*
{xrst_begin eigen_det.cpp}
Using Eigen To Compute Determinant: Example and Test
####################################################
{xrst_literal
// BEGIN C++
// END C++
}
{xrst_end eigen_det.cpp}
*/
// BEGIN C++
# 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;
using Eigen::Index;
//
typedef Matrix< double , Dynamic, Dynamic > matrix;
typedef Matrix< AD<double> , Dynamic, Dynamic > a_matrix;
//
typedef CppAD::eigen_vector<double> vector;
typedef CppAD::eigen_vector< AD<double> > a_vector;
//
// 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(size_t i = 0; i < size; i++)
{ for(size_t 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
Index Size = Index(size);
a_matrix a_X(Size, Size);
matrix X(Size, Size);
for(size_t i = 0; i < size; i++)
{ for(size_t j = 0; j < size; j++)
{ Index I = Index(i);
Index J = Index(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(size_t i = 0; i < size; i++)
{ for(size_t j = 0; j < size; j++)
{ Index I = Index(i);
Index J = Index(j);
ok &= NearEqual(jac[i * size + j], inv_X(J, I), eps, eps);
}
}
return ok;
}
// END C++
|
