Prev Next Index-> contents reference index search external Up-> CppAD AD ADValued atomic atomic_base atomic_eigen_cholesky.cpp ADValued-> Arithmetic unary_standard_math binary_math CondExp Discrete numeric_limits atomic atomic-> checkpoint atomic_base atomic_base-> atomic_ctor atomic_option atomic_afun atomic_forward atomic_reverse atomic_for_sparse_jac atomic_rev_sparse_jac atomic_for_sparse_hes atomic_rev_sparse_hes atomic_base_clear atomic_get_started.cpp atomic_norm_sq.cpp atomic_reciprocal.cpp atomic_set_sparsity.cpp atomic_tangent.cpp atomic_eigen_mat_mul.cpp atomic_eigen_mat_inv.cpp atomic_eigen_cholesky.cpp atomic_mat_mul.cpp atomic_eigen_cholesky.cpp-> cholesky_theory atomic_eigen_cholesky.hpp Headings-> Description Contents Use Atomic Function ---..Constructor

$\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}} }$
Atomic Eigen Cholesky Factorization: Example and Test

Description
The ADFun function object f for this example is $$f(x) = \R{chol} \left( \begin{array}{cc} x_0 & x_1 \\ x_1 & x_2 \end{array} \right) = \frac{1}{ \sqrt{x_0} } \left( \begin{array}{cc} x_0 & 0 \\ x_1 & \sqrt{ x_0 x_2 - x_1 x_1 } \end{array} \right)$$ where the matrix is positive definite; i.e., $x_0 > 0$, $x_2 > 0$ and $x_0 x_2 - x_1 x_1 > 0$.

Contents
 cholesky_theory AD Theory for Cholesky Factorization atomic_eigen_cholesky.hpp Atomic Eigen Cholesky Factorization Class

Use Atomic Function
# include <cppad/cppad.hpp>

bool eigen_cholesky(void)
{
typedef double scalar;
//
bool ok    = true;
scalar eps = 10. * std::numeric_limits<scalar>::epsilon();
//

Constructor
     // -------------------------------------------------------------------
// object that computes cholesky factor of a matrix
atomic_eigen_cholesky<scalar> cholesky;
// -------------------------------------------------------------------
// declare independent variable vector x
size_t n = 3;
// -------------------------------------------------------------------
// A = [ x[0]  x[1] ]
//     [ x[1]  x[2] ]
size_t nr  = 2;
// -------------------------------------------------------------------
// use atomic operation to L such that A = L * L^T
// -------------------------------------------------------------------
// declare the dependent variable vector y
size_t m = 3;
// -------------------------------------------------------------------
// check zero order forward mode
x[0] = 2.0;
x[1] = 0.5;
x[2] = 5.0;
y   = f.Forward(0, x);
scalar check;
check = std::sqrt( x[0] );
ok   &= NearEqual(y[0], check, eps, eps);
check = x[1] / std::sqrt( x[0] );
ok   &= NearEqual(y[1], check, eps, eps);
check = std::sqrt( x[2] - x[1] * x[1] / x[0] );
ok   &= NearEqual(y[2], check, eps, eps);
// -------------------------------------------------------------------
// check first order forward mode
//
// partial w.r.t. x[0]
x1[0] = 1.0;
x1[1] = 0.0;
x1[2] = 0.0;
//
y1    = f.Forward(1, x1);
check = 1.0 / (2.0 * std::sqrt( x[0] ) );
ok   &= NearEqual(y1[0], check, eps, eps);
//
check = - x[1] / (2.0 * x[0] * std::sqrt( x[0] ) );
ok   &= NearEqual(y1[1], check, eps, eps);
//
check = std::sqrt( x[2] - x[1] * x[1] / x[0] );
check = x[1] * x[1] / (x[0] * x[0] * 2.0 * check);
ok   &= NearEqual(y1[2], check, eps, eps);
//
// partial w.r.t. x[1]
x1[0] = 0.0;
x1[1] = 1.0;
x1[2] = 0.0;
//
y1    = f.Forward(1, x1);
ok   &= NearEqual(y1[0], 0.0, eps, eps);
//
check = 1.0 / std::sqrt( x[0] );
ok   &= NearEqual(y1[1], check, eps, eps);
//
check = std::sqrt( x[2] - x[1] * x[1] / x[0] );
check = - 2.0 * x[1] / (2.0 * check * x[0] );
ok   &= NearEqual(y1[2], check, eps, eps);
//
// partial w.r.t. x[2]
x1[0] = 0.0;
x1[1] = 0.0;
x1[2] = 1.0;
//
y1    = f.Forward(1, x1);
ok   &= NearEqual(y1[0], 0.0, eps, eps);
ok   &= NearEqual(y1[1], 0.0, eps, eps);
//
check = std::sqrt( x[2] - x[1] * x[1] / x[0] );
check = 1.0 / (2.0 * check);
ok   &= NearEqual(y1[2], check, eps, eps);
// -------------------------------------------------------------------
// check second order forward mode
//
// second partial w.r.t x[2]
x2[0] = 0.0;
x2[1] = 0.0;
x2[2] = 0.0;
y2    = f.Forward(2, x2);
ok   &= NearEqual(y2[0], 0.0, eps, eps);
ok   &= NearEqual(y2[1], 0.0, eps, eps);
//
check = std::sqrt( x[2] - x[1] * x[1] / x[0] );  // funciton value
check = - 1.0 / ( 4.0 * check * check * check ); // second derivative
check = 0.5 * check;                             // taylor coefficient
ok   &= NearEqual(y2[2], check, eps, eps);
// -------------------------------------------------------------------
// check first order reverse mode
w[0] = 0.0;
w[1] = 0.0;
w[2] = 1.0;
d1w  = f.Reverse(1, w);
//
// partial of f[2] w.r.t x[0]
scalar f2    = std::sqrt( x[2] - x[1] * x[1] / x[0] );
scalar f2_x0 = x[1] * x[1] / (2.0 * f2 * x[0] * x[0] );
ok          &= NearEqual(d1w[0], f2_x0, eps, eps);
//
// partial of f[2] w.r.t x[1]
scalar f2_x1 = - x[1] / (f2 * x[0] );
ok          &= NearEqual(d1w[1], f2_x1, eps, eps);
//
// partial of f[2] w.r.t x[2]
scalar f2_x2 = 1.0 / (2.0 * f2 );
ok          &= NearEqual(d1w[2], f2_x2, eps, eps);
// -------------------------------------------------------------------
// check second order reverse mode
d2w  = f.Reverse(2, w);
//
// check first order results
ok &= NearEqual(d2w[0 * 2 + 0], f2_x0, eps, eps);
ok &= NearEqual(d2w[1 * 2 + 0], f2_x1, eps, eps);
ok &= NearEqual(d2w[2 * 2 + 0], f2_x2, eps, eps);
//
// check second order results
scalar f2_x2_x0 = - 0.5 * f2_x0 / (f2 * f2 );
ok             &= NearEqual(d2w[0 * 2 + 1], f2_x2_x0, eps, eps);
scalar f2_x2_x1 = - 0.5 * f2_x1 / (f2 * f2 );
ok             &= NearEqual(d2w[1 * 2 + 1], f2_x2_x1, eps, eps);
scalar f2_x2_x2 = - 0.5 * f2_x2 / (f2 * f2 );
ok             &= NearEqual(d2w[2 * 2 + 1], f2_x2_x2, eps, eps);
// -------------------------------------------------------------------
// check third order reverse mode
d3w  = f.Reverse(3, w);
//
// check first order results
ok &= NearEqual(d3w[0 * 3 + 0], f2_x0, eps, eps);
ok &= NearEqual(d3w[1 * 3 + 0], f2_x1, eps, eps);
ok &= NearEqual(d3w[2 * 3 + 0], f2_x2, eps, eps);
//
// check second order results
ok             &= NearEqual(d3w[0 * 3 + 1], f2_x2_x0, eps, eps);
ok             &= NearEqual(d3w[1 * 3 + 1], f2_x2_x1, eps, eps);
ok             &= NearEqual(d3w[2 * 3 + 1], f2_x2_x2, eps, eps);
// -------------------------------------------------------------------
scalar f2_x2_x2_x0 = - 0.5 * f2_x2_x0 / (f2 * f2);
f2_x2_x2_x0 += f2_x2 * f2_x0 / (f2 * f2 * f2);
ok          &= NearEqual(d3w[0 * 3 + 2], 0.5 * f2_x2_x2_x0, eps, eps);
scalar f2_x2_x2_x1 = - 0.5 * f2_x2_x1 / (f2 * f2);
f2_x2_x2_x1 += f2_x2 * f2_x1 / (f2 * f2 * f2);
ok          &= NearEqual(d3w[1 * 3 + 2], 0.5 * f2_x2_x2_x1, eps, eps);
scalar f2_x2_x2_x2 = - 0.5 * f2_x2_x2 / (f2 * f2);
f2_x2_x2_x2 += f2_x2 * f2_x2 / (f2 * f2 * f2);
ok          &= NearEqual(d3w[2 * 3 + 2], 0.5 * f2_x2_x2_x2, eps, eps);
return ok;
}

Input File: example/atomic/eigen_cholesky.cpp