Last change on this file since 3829 was 3829, checked in by bradbell, 3 years ago

merge to branch: trunk
start hash code: 02e2cde4b6486747926917328764ff776ec7c252

Date: Tue Sep 13 14:15:14 2016 -0700

Date: Tue Sep 13 14:01:28 2016 -0700

eigen_cholesky.hpp: fix bug in reverse mode for orders >= 3.
eigen_cholesky.cpp: test third order reverse mode.

commit 0337e557eef27710c445fc1cf9d000ed0234c98e
Date: Tue Aug 30 20:34:59 2016 -0700

Test and fix atomic_eigen_cholesky reverse mode calculation.

commit b2e44a78fc40b39d4451514a735ee42f9f4b51a5
Date: Tue Aug 30 07:48:43 2016 -0700

2. Add the atomic_eigen_cholesky example (under construction).
3. In Eigen atomic examples, mention that we are not checking if certain reverse mode terms need be recalculated.

eigen_mat_inv.hpp: change f_sum_ -> f_sum (not member variable).

File size: 12.1 KB
Line
1// $Id$
4
5/* --------------------------------------------------------------------------
7
9the terms of the
10                    Eclipse Public License Version 1.0.
11
12A copy of this license is included in the COPYING file of this distribution.
14-------------------------------------------------------------------------- */
15
16/*
17$begin atomic_eigen_mat_inv.hpp$$18spell 19 Eigen 20 Taylor 21$$ 22 23$section Atomic Eigen Matrix Inversion Class$$24 25head Purpose$$
26Construct an atomic operation that computes the matrix inverse
27$latex R = A^{-1}$$28for any positive integer latex p$$ 29and invertible matrix$latex A \in \B{R}^{p \times p}$$. 30 31head Matrix Dimensions$$
32This example puts the matrix dimension $latex p$$33in the atomic function arguments, 34instead of the cref/constructor/atomic_ctor/$$, 35so it can be different for different calls to the atomic function. 36 37$head Theory$$38 39subhead Forward$$
40The zero order forward mode Taylor coefficient is give by
41$latex $42 R_0 = A_0^{-1} 43$$$44For latex k = 1 , \ldots$$, 45the$th k$$order Taylor coefficient of latex A R$$ is given by
46$latex $47 0 = \sum_{\ell=0}^k A_\ell R_{k-\ell} 48$ $$49Solving for latex R_k$$ in terms of the coefficients 50for$latex A$$and the lower order coefficients for latex R$$ we have
51$latex $52 R_k = - R_0 \left( \sum_{\ell=1}^k A_\ell R_{k-\ell} \right) 53$ $$54Furthermore, once we have latex R_k$$ we can compute the sum using 55$latex $56 A_0 R_k = - \left( \sum_{\ell=1}^k A_\ell R_{k-\ell} \right) 57$ $$58 59 60subhead Product of Three Matrices$$
61Suppose $latex \bar{E}$$is the derivative of the 62scalar value function latex s(E)$$ with respect to$latex E$$; i.e., 63latex $64 \bar{E}_{i,j} = \frac{ \partial s } { \partial E_{i,j} } 65$$$
66Also suppose that $latex t$$is a scalar valued argument and 67latex $68 E(t) = B(t) C(t) D(t) 69$$$ 70It follows that 71$latex $72 E'(t) = B'(t) C(t) D(t) + B(t) C'(t) D(t) + B(t) C(t) D'(t) 73$ $$74 75latex $76 (s \circ E)'(t) 77 = 78 \R{tr} [ \bar{E}^\R{T} E'(t) ] 79$$$
80$latex $81 = 82 \R{tr} [ \bar{E}^\R{T} B'(t) C(t) D(t) ] + 83 \R{tr} [ \bar{E}^\R{T} B(t) C'(t) D(t) ] + 84 \R{tr} [ \bar{E}^\R{T} B(t) C(t) D'(t) ] 85$ $$86latex $87 = 88 \R{tr} [ B(t) D(t) \bar{E}^\R{T} B'(t) ] + 89 \R{tr} [ D(t) \bar{E}^\R{T} B(t) C'(t) ] + 90 \R{tr} [ \bar{E}^\R{T} B(t) C(t) D'(t) ] 91$$$ 92$latex $93 \bar{B} = \bar{E} (C D)^\R{T} \W{,} 94 \bar{C} = B^\R{T} \bar{E} D^\R{T} \W{,} 95 \bar{D} = (B C)^\R{T} \bar{E} 96$ $$97 98subhead Reverse$$
99For $latex k > 0$$, reverse mode 100eliminates latex R_k$$ and expresses the function values 101$latex s$$in terms of the coefficients of latex A$$
102and the lower order coefficients of $latex R$$. 103The effect on latex \bar{R}_0$$ 104(of eliminating$latex R_k$$) is 105latex $106\bar{R}_0 107= \bar{R}_0 - \bar{R}_k \left( \sum_{\ell=1}^k A_\ell R_{k-\ell} \right)^\R{T} 108= \bar{R}_0 + \bar{R}_k ( A_0 R_k )^\R{T} 109$$$
110For $latex \ell = 1 , \ldots , k$$, 111the effect on latex \bar{R}_{k-\ell}$$ and$latex A_\ell$$112(of eliminating latex R_k$$) is
113$latex $114\bar{A}_\ell = \bar{A}_\ell - R_0^\R{T} \bar{R}_k R_{k-\ell}^\R{T} 115$ $$116latex $117\bar{R}_{k-\ell} = \bar{R}_{k-\ell} - ( R_0 A_\ell )^\R{T} \bar{R}_k 118$$$ 119We note that 120$latex $121 R_0 '(t) A_0 (t) + R_0 (t) A_0 '(t) = 0 122$ $$123latex $124 R_0 '(t) = - R_0 (t) A_0 '(t) R_0 (t) 125$$$
126The reverse mode formula that eliminates $latex R_0$$is 127latex $128 \bar{A}_0 129 = \bar{A}_0 - R_0^\R{T} \bar{R}_0 R_0^\R{T} 130$$$ 131 132$nospell
133
134$head Start Class Definition$$135srccode%cpp% */ 136# include <cppad/cppad.hpp> 137# include <Eigen/Core> 138# include <Eigen/LU> 139 140 141 142/* %$$ 143$head Public$$144 145subhead Types$$
146$srccode%cpp% */ 147namespace { // BEGIN_EMPTY_NAMESPACE 148 149template <class Base> 150class atomic_eigen_mat_inv : public CppAD::atomic_base<Base> { 151public: 152 // ----------------------------------------------------------- 153 // type of elements during calculation of derivatives 154 typedef Base scalar; 155 // type of elements during taping 156 typedef CppAD::AD<scalar> ad_scalar; 157 // type of matrix during calculation of derivatives 158 typedef Eigen::Matrix< 159 scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> matrix; 160 // type of matrix during taping 161 typedef Eigen::Matrix< 162 ad_scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > ad_matrix; 163/* %$$164subhead Constructor$$ 165$srccode%cpp% */
166        // constructor
172/* %$$173subhead op$$
174$srccode%cpp% */ 175 // use atomic operation to invert an AD matrix 176 ad_matrix op(const ad_matrix& arg) 177 { size_t nr = size_t( arg.rows() ); 178 size_t ny = nr * nr; 179 size_t nx = 1 + ny; 180 assert( nr == size_t( arg.cols() ) ); 181 // ------------------------------------------------------------------- 182 // packed version of arg 183 CPPAD_TESTVECTOR(ad_scalar) packed_arg(nx); 184 packed_arg[0] = ad_scalar( nr ); 185 for(size_t i = 0; i < ny; i++) 186 packed_arg[1 + i] = arg.data()[i]; 187 // ------------------------------------------------------------------- 188 // packed version of result = arg^{-1}. 189 // This is an atomic_base function call that CppAD uses to 190 // store the atomic operation on the tape. 191 CPPAD_TESTVECTOR(ad_scalar) packed_result(ny); 192 (*this)(packed_arg, packed_result); 193 // ------------------------------------------------------------------- 194 // unpack result matrix 195 ad_matrix result(nr, nr); 196 for(size_t i = 0; i < ny; i++) 197 result.data()[i] = packed_result[i]; 198 return result; 199 } 200 /* %$$201head Private$$ 202 203$subhead Variables$$204srccode%cpp% */ 205private: 206 // ------------------------------------------------------------- 207 // one forward mode vector of matrices for argument and result 208 CppAD::vector<matrix> f_arg_, f_result_; 209 // one reverse mode vector of matrices for argument and result 210 CppAD::vector<matrix> r_arg_, r_result_; 211 // ------------------------------------------------------------- 212/* %$$
213$subhead forward$$214srccode%cpp% */ 215 // forward mode routine called by CppAD 216 virtual bool forward( 217 // lowest order Taylor coefficient we are evaluating 218 size_t p , 219 // highest order Taylor coefficient we are evaluating 220 size_t q , 221 // which components of x are variables 222 const CppAD::vector<bool>& vx , 223 // which components of y are variables 224 CppAD::vector<bool>& vy , 225 // tx [ j * (q+1) + k ] is x_j^k 226 const CppAD::vector<scalar>& tx , 227 // ty [ i * (q+1) + k ] is y_i^k 228 CppAD::vector<scalar>& ty 229 ) 230 { size_t n_order = q + 1; 231 size_t nr = size_t( CppAD::Integer( tx[ 0 * n_order + 0 ] ) ); 232 size_t ny = nr * nr; 233 size_t nx = 1 + ny; 234 assert( vx.size() == 0 || nx == vx.size() ); 235 assert( vx.size() == 0 || ny == vy.size() ); 236 assert( nx * n_order == tx.size() ); 237 assert( ny * n_order == ty.size() ); 238 // 239 // ------------------------------------------------------------------- 240 // make sure f_arg_ and f_result_ are large enough 241 assert( f_arg_.size() == f_result_.size() ); 242 if( f_arg_.size() < n_order ) 243 { f_arg_.resize(n_order); 244 f_result_.resize(n_order); 245 // 246 for(size_t k = 0; k < n_order; k++) 247 { f_arg_[k].resize(nr, nr); 248 f_result_[k].resize(nr, nr); 249 } 250 } 251 // ------------------------------------------------------------------- 252 // unpack tx into f_arg_ 253 for(size_t k = 0; k < n_order; k++) 254 { // unpack arg values for this order 255 for(size_t i = 0; i < ny; i++) 256 f_arg_[k].data()[i] = tx[ (1 + i) * n_order + k ]; 257 } 258 // ------------------------------------------------------------------- 259 // result for each order 260 // (we could avoid recalculting f_result_[k] for k=0,...,p-1) 261 // 262 f_result_[0] = f_arg_[0].inverse(); 263 for(size_t k = 1; k < n_order; k++) 264 { // initialize sum 265 matrix f_sum = matrix::Zero(nr, nr); 266 // compute sum 267 for(size_t ell = 1; ell <= k; ell++) 268 f_sum -= f_arg_[ell] * f_result_[k-ell]; 269 // result_[k] = arg_[0]^{-1} * sum_ 270 f_result_[k] = f_result_[0] * f_sum; 271 } 272 // ------------------------------------------------------------------- 273 // pack result_ into ty 274 for(size_t k = 0; k < n_order; k++) 275 { for(size_t i = 0; i < ny; i++) 276 ty[ i * n_order + k ] = f_result_[k].data()[i]; 277 } 278 // ------------------------------------------------------------------- 279 // check if we are computing vy 280 if( vx.size() == 0 ) 281 return true; 282 // ------------------------------------------------------------------ 283 // This is a very dumb algorithm that over estimates which 284 // elements of the inverse are variables (which is not efficient). 285 bool var = false; 286 for(size_t i = 0; i < ny; i++) 287 var |= vx[1 + i]; 288 for(size_t i = 0; i < ny; i++) 289 vy[i] = var; 290 return true; 291 } 292/* %$$ 293$subhead reverse$$294srccode%cpp% */ 295 // reverse mode routine called by CppAD 296 virtual bool reverse( 297 // highest order Taylor coefficient that we are computing derivative of 298 size_t q , 299 // forward mode Taylor coefficients for x variables 300 const CppAD::vector<double>& tx , 301 // forward mode Taylor coefficients for y variables 302 const CppAD::vector<double>& ty , 303 // upon return, derivative of G[ F[ {x_j^k} ] ] w.r.t {x_j^k} 304 CppAD::vector<double>& px , 305 // derivative of G[ {y_i^k} ] w.r.t. {y_i^k} 306 const CppAD::vector<double>& py 307 ) 308 { size_t n_order = q + 1; 309 size_t nr = size_t( CppAD::Integer( tx[ 0 * n_order + 0 ] ) ); 310 size_t ny = nr * nr; 311 size_t nx = 1 + ny; 312 // 313 assert( nx * n_order == tx.size() ); 314 assert( ny * n_order == ty.size() ); 315 assert( px.size() == tx.size() ); 316 assert( py.size() == ty.size() ); 317 // ------------------------------------------------------------------- 318 // make sure f_arg_ is large enough 319 assert( f_arg_.size() == f_result_.size() ); 320 // must have previous run forward with order >= n_order 321 assert( f_arg_.size() >= n_order ); 322 // ------------------------------------------------------------------- 323 // make sure r_arg_, r_result_ are large enough 324 assert( r_arg_.size() == r_result_.size() ); 325 if( r_arg_.size() < n_order ) 326 { r_arg_.resize(n_order); 327 r_result_.resize(n_order); 328 // 329 for(size_t k = 0; k < n_order; k++) 330 { r_arg_[k].resize(nr, nr); 331 r_result_[k].resize(nr, nr); 332 } 333 } 334 // ------------------------------------------------------------------- 335 // unpack tx into f_arg_ 336 for(size_t k = 0; k < n_order; k++) 337 { // unpack arg values for this order 338 for(size_t i = 0; i < ny; i++) 339 f_arg_[k].data()[i] = tx[ (1 + i) * n_order + k ]; 340 } 341 // ------------------------------------------------------------------- 342 // unpack py into r_result_ 343 for(size_t k = 0; k < n_order; k++) 344 { for(size_t i = 0; i < ny; i++) 345 r_result_[k].data()[i] = py[ i * n_order + k ]; 346 } 347 // ------------------------------------------------------------------- 348 // initialize r_arg_ as zero 349 for(size_t k = 0; k < n_order; k++) 350 r_arg_[k] = matrix::Zero(nr, nr); 351 // ------------------------------------------------------------------- 352 // matrix reverse mode calculation 353 // 354 for(size_t k1 = n_order; k1 > 1; k1--) 355 { size_t k = k1 - 1; 356 // bar{R}_0 = bar{R}_0 + bar{R}_k (A_0 R_k)^T 357 r_result_[0] += 358 r_result_[k] * f_result_[k].transpose() * f_arg_[0].transpose(); 359 // 360 for(size_t ell = 1; ell <= k; ell++) 361 { // bar{A}_l = bar{A}_l - R_0^T bar{R}_k R_{k-l}^T 362 r_arg_[ell] -= f_result_[0].transpose() 363 * r_result_[k] * f_result_[k-ell].transpose(); 364 // bar{R}_{k-l} = bar{R}_{k-1} - (R_0 A_l)^T bar{R}_k 365 r_result_[k-ell] -= f_arg_[ell].transpose() 366 * f_result_[0].transpose() * r_result_[k]; 367 } 368 } 369 r_arg_[0] -= 370 f_result_[0].transpose() * r_result_[0] * f_result_[0].transpose(); 371 // ------------------------------------------------------------------- 372 // pack r_arg into px 373 for(size_t k = 0; k < n_order; k++) 374 { for(size_t i = 0; i < ny; i++) 375 px[ (1 + i) * n_order + k ] = r_arg_[k].data()[i]; 376 } 377 // 378 return true; 379 } 380/* %$$
381$head End Class Definition$$382srccode%cpp% */ 383}; // End of atomic_eigen_mat_inv class 384 385} // END_EMPTY_NAMESPACE 386/* %$$ 387$$comment end nospell$$ 388$end