source: trunk/cppad/example/eigen_mat_inv.hpp @ 3829

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

merge to branch: trunk
from repository: https://github.com/coin-or/CppAD
start hash code: 02e2cde4b6486747926917328764ff776ec7c252
end hash code: f99038275e49c519b8136bbc8ecab93ed709cad4

commit f99038275e49c519b8136bbc8ecab93ed709cad4
Author: Brad Bell <bradbell@…>
Date: Tue Sep 13 14:15:14 2016 -0700

  1. Advance version to cppad-20160913.

commit 941b6adbe202e6e3c9a7add6fd83a4de4858bef2
Author: Brad Bell <bradbell@…>
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
Author: Brad Bell <bradbell@…>
Date: Tue Aug 30 20:34:59 2016 -0700

Test and fix atomic_eigen_cholesky reverse mode calculation.

commit b2e44a78fc40b39d4451514a735ee42f9f4b51a5
Author: Brad Bell <bradbell@…>
Date: Tue Aug 30 07:48:43 2016 -0700

  1. Advance to cppad-20160830.
  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$
2# ifndef CPPAD_EXAMPLE_EIGEN_MAT_INV_HPP
3# define CPPAD_EXAMPLE_EIGEN_MAT_INV_HPP
4
5/* --------------------------------------------------------------------------
6CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-16 Bradley M. Bell
7
8CppAD is distributed under multiple licenses. This distribution is under
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.
13Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
14-------------------------------------------------------------------------- */
15
16/*
17$begin atomic_eigen_mat_inv.hpp$$
18$spell
19        Eigen
20        Taylor
21$$
22
23$section Atomic Eigen Matrix Inversion Class$$
24
25$head 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
31$head 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
39$subhead 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
60$subhead 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.,
63$latex \[
64        \bar{E}_{i,j} = \frac{ \partial s } { \partial E_{i,j} }
65\] $$
66Also suppose that $latex t$$ is a scalar valued argument and
67$latex \[
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
75$latex \[
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\] $$
86$latex \[
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
98$subhead 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
105$latex \[
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\] $$
116$latex \[
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\] $$
123$latex \[
124        R_0 '(t) = - R_0 (t) A_0 '(t) R_0 (t)
125\] $$
126The reverse mode formula that eliminates $latex R_0$$ is
127$latex \[
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$$
135$srccode%cpp% */
136# include <cppad/cppad.hpp>
137# include <Eigen/Core>
138# include <Eigen/LU>
139
140
141
142/* %$$
143$head Public$$
144
145$subhead 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/* %$$
164$subhead Constructor$$
165$srccode%cpp% */
166        // constructor
167        atomic_eigen_mat_inv(void) : CppAD::atomic_base<Base>(
168                "atom_eigen_mat_inv"                             ,
169                CppAD::atomic_base<Base>::set_sparsity_enum
170        )
171        { }
172/* %$$
173$subhead 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        /* %$$
201$head Private$$
202
203$subhead Variables$$
204$srccode%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$$
214$srccode%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$$
294$srccode%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$$
382$srccode%cpp% */
383}; // End of atomic_eigen_mat_inv class
384
385}  // END_EMPTY_NAMESPACE
386/* %$$
387$$ $comment end nospell$$
388$end
389*/
390
391
392# endif
Note: See TracBrowser for help on using the repository browser.