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

Last change on this file since 3810 was 3810, checked in by bradbell, 4 years ago

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

commit f77cd2e443b8d14f354f97efc4c26369227dfe81
Author: Brad Bell <bradbell@…>
Date: Sat Mar 26 06:47:33 2016 -0700

  1. Add reverse mode to atomic_eigen_mat_inv.
  2. Advance version to cppad-20160326.


eigen_mat_mul.cpp: impove reverse mode comments.

commit 94c81399191d06ae9cef81c271f520b0e78f079d
Author: Brad Bell <bradbell@…>
Date: Fri Mar 25 12:16:03 2016 -0700

whats_new_16.omh: add comment about changes to eigen_mat_mul.

commit d647288c9e022f16ef94229db612a774808a5de3
Author: Brad Bell <bradbell@…>
Date: Fri Mar 25 11:59:45 2016 -0700

Replace pack and unpack by member function that inverts matrices.

commit 6fba95d90e50590abd438946f47025e0a6e34c2c
Author: Brad Bell <bradbell@…>
Date: Fri Mar 25 11:45:43 2016 -0700

Replace pack and unpack by member function that multiplies matrices.

commit cac97c115abad9fba0d0c3706f4bc0e5a7f617f0
Author: Brad Bell <bradbell@…>
Date: Fri Mar 25 11:11:16 2016 -0700

Change atomic_eigen_mat_div -> atomic_eigen_mat_inv.

commit 1f7e52a4d96986f7be2ad2d6abeafc086c4f7622
Author: Brad Bell <bradbell@…>
Date: Fri Mar 25 08:09:39 2016 -0700

  1. Implement and test forward for atomic_eigen_mat_div.
  2. Advance version to 20160325.


eigen_mat_mul.hpp: edit comments and documentation.

commit e953f669d8e9e122dc2f182c28c051c9ba75c705
Author: Brad Bell <bradbell@…>
Date: Fri Mar 25 04:32:45 2016 -0700

eigen_mat_mul.cpp: No longer under construction.

commit 3a62044eaebb6393676276e5b843fbe5db4f45a2
Author: Brad Bell <bradbell@…>
Date: Fri Mar 25 04:07:18 2016 -0700

eigen_mat_mul.hpp: add assert for f_left_ and f_right_ in reverse.

commit 721d64783547487fe69100c7190628efab7149c3
Author: Brad Bell <bradbell@…>
Date: Fri Mar 25 03:59:05 2016 -0700

eigen_mat_mul.hpp: typo Publice -> Public.

File size: 11.2 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$$
26For fixed positive integer $latex p$$,
27construct and atomic operation that solve the computes the matrix inverse
28$latex R = A^{-1}$$
29for any invertible $latex A \in \B{R}^{p \times p}$$.
30
31$head Theory$$
32
33$head Forward$$
34For $latex k = 0 , \ldots$$, the $th k$$ order Taylor coefficient is given by
35$latex \[
36        I_k = \sum_{\ell=0}^k A_\ell R_{k-\ell}
37\] $$
38$latex \[
39        R_k = A_0^{-1} \left( I_k - \sum_{\ell=1}^k A_\ell R_{k-\ell} \right)
40\] $$
41where $latex I_k$$ is the $th k$$ order Taylor coefficient for the
42identity matrix. i.e., $latex I_k$$ is the identity matrix if $latex k = 0$$
43and is the zero matrix if $latex k \neq 0$$.
44
45$head Reverse$$
46We use $latex \bar{R}_k$$ for the partial of the scalar final result
47with respect to $latex R_k$$.
48Note that $latex R_0 = A_0^{-1}$$.
49The back-propagation algorithm that eliminates $latex R_k$$,
50for $latex k > 0$$, is
51for $latex \ell = 1 , \ldots , k$$
52$latex \[
53\bar{A}_\ell = \bar{A}_\ell - R_0 \bar{R}_k R_{k-\ell}^\R{T}
54\] $$
55$latex \[
56\bar{R}_{k-\ell} = \bar{R}_{k-\ell} - R_0 A_\ell^\R{T} \bar{R}_k
57\] $$
58$latex \[
59\bar{R}_0  = \bar{R}_0 + \bar{R}_k
60        \left( I_k - \sum_{\ell=1}^k A_\ell R_{k-\ell} \right)^\R{T}
61\] $$
62$latex \[
63\bar{R}_0  = \bar{R}_0 + \bar{R}_k ( A_0 R_k )^\R{T}
64\] $$
65The back-propagation algorithm that eliminates $latex R_0$$ is
66$latex \[
67        \bar{A}_0 = \bar{A}_0 - R_0^\R{T} \bar{R}_0 R_0^\R{T}
68\]$$
69
70$nospell
71
72$head Start Class Definition$$
73$srccode%cpp% */
74# include <cppad/cppad.hpp>
75# include <Eigen/Core>
76# include <Eigen/LU>
77
78
79
80/* %$$
81$head Public$$
82
83$subhead Types$$
84$srccode%cpp% */
85namespace { // BEGIN_EMPTY_NAMESPACE
86
87template <class Base>
88class atomic_eigen_mat_inv : public CppAD::atomic_base<Base> {
89public:
90        // -----------------------------------------------------------
91        // type of elements during calculation of derivatives
92        typedef Base              scalar;
93        // type of elements during taping
94        typedef CppAD::AD<scalar> ad_scalar;
95        // type of matrix during calculation of derivatives
96        typedef Eigen::Matrix<
97                scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>     matrix;
98        // type of matrix during taping
99        typedef Eigen::Matrix<
100                ad_scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > ad_matrix;
101/* %$$
102$subhead Constructor$$
103$srccode%cpp% */
104        // constructor
105        atomic_eigen_mat_inv(
106                // number of rows and columns in argument
107                const size_t nr
108        ) :
109        CppAD::atomic_base<Base>(
110                "atom_eigen_mat_inv"                             ,
111                CppAD::atomic_base<Base>::set_sparsity_enum
112        ) ,
113        nr_(  nr  )       ,
114        nx_( nr * nr )    ,
115        f_sum_( nr, nr )
116        { }
117/* %$$
118$subhead op$$
119$srccode%cpp% */
120        // use atomic operation to invert an AD matrix
121        ad_matrix op(const ad_matrix& arg)
122        {       assert( nr_ == size_t( arg.rows() ) );
123                assert( nr_ == size_t( arg.cols() ) );
124
125                // -------------------------------------------------------------------
126                // packed version of arg
127                CPPAD_TESTVECTOR(ad_scalar) packed_arg(nx_);
128                for(size_t i = 0; i < nx_; i++)
129                        packed_arg[i] = arg.data()[i];
130
131                // -------------------------------------------------------------------
132                // packed version of result = arg^{-1}
133                CPPAD_TESTVECTOR(ad_scalar) packed_result(nx_);
134                (*this)(packed_arg, packed_result);
135
136                // -------------------------------------------------------------------
137                // unpack result matrix
138                ad_matrix result(nr_, nr_);
139                for(size_t i = 0; i < nx_; i++)
140                        result.data()[i] = packed_result[ i ];
141
142                return result;
143        }
144/* %$$
145$head Private$$
146
147$subhead Variables$$
148$srccode%cpp% */
149private:
150        // -------------------------------------------------------------
151        // number of of rows in argument and result
152        const size_t nr_;
153        // size of the domain (and range) space
154        const size_t nx_;
155        // -------------------------------------------------------------
156        // one forward mode vector of matrices for argument and result
157        CppAD::vector<matrix> f_arg_, f_result_;
158        // matrix used for forward mode summation
159        matrix f_sum_;
160        // one reverse mode vector of matrices for argument and result
161        CppAD::vector<matrix> r_arg_, r_result_;
162        // -------------------------------------------------------------
163/* %$$
164$subhead rows$$
165$srccode%cpp% */
166        // convert from int to size_t
167        static size_t rows(const matrix& x)
168        {       return size_t( x.rows() ); }
169        static size_t rows(const ad_matrix& x)
170        {       return size_t( x.rows() ); }
171/* %$$
172$subhead cols$$
173$srccode%cpp% */
174        // convert from int to size_t
175        static size_t cols(const matrix& x)
176        {       return size_t( x.cols() ); }
177        static size_t cols(const ad_matrix& x)
178        {       return size_t( x.cols() ); }
179/* %$$
180$subhead forward$$
181$srccode%cpp% */
182        // forward mode routine called by CppAD
183        virtual bool forward(
184                // lowest order Taylor coefficient we are evaluating
185                size_t                          p ,
186                // highest order Taylor coefficient we are evaluating
187                size_t                          q ,
188                // which components of x are variables
189                const CppAD::vector<bool>&      vx ,
190                // which components of y are variables
191                CppAD::vector<bool>&            vy ,
192                // tx [ j * (q+1) + k ] is x_j^k
193                const CppAD::vector<scalar>&    tx ,
194                // ty [ i * (q+1) + k ] is y_i^k
195                CppAD::vector<scalar>&          ty
196        )
197        {       size_t n_order = q + 1;
198                assert( vx.size() == 0 || nx_ == vx.size() );
199                assert( vx.size() == 0 || nx_ == vy.size() );
200                assert( nx_ * n_order == tx.size() );
201                assert( nx_ * n_order == ty.size() );
202                //
203                // -------------------------------------------------------------------
204                // make sure f_arg_ and f_result_ are large enough
205                assert( f_arg_.size() == f_result_.size() );
206                if( f_arg_.size() < n_order )
207                {       f_arg_.resize(n_order);
208                        f_result_.resize(n_order);
209                        //
210                        for(size_t k = 0; k < n_order; k++)
211                        {       f_arg_[k].resize(nr_, nr_);
212                                f_result_[k].resize(nr_, nr_);
213                        }
214                }
215                // -------------------------------------------------------------------
216                // unpack tx into f_arg_
217                for(size_t k = 0; k < n_order; k++)
218                {       // unpack arg values for this order
219                        for(size_t i = 0; i < nx_; i++)
220                                f_arg_[k].data()[i] = tx[ i * n_order + k ];
221                }
222                // -------------------------------------------------------------------
223                // result for each order
224                //
225                // compute LU factorixation of arg_[0]
226                Eigen::FullPivLU<matrix> lu_arg( f_arg_[0] );
227                for(size_t k = 0; k < n_order; k++)
228                {       // initialize sum
229                        if( k == 0 )
230                                f_sum_ = matrix::Identity(nr_, nr_);
231                        else
232                                f_sum_ = matrix::Zero(nr_, nr_);
233                        // compute sum
234                        for(size_t ell = 1; ell <= k; ell++)
235                                f_sum_ -= f_arg_[ell] * f_result_[k-ell];
236                        // solve arg_[0] * result_[k] = sum
237                        // result_[k] = arg_[0]^{-1} * sum
238                        f_result_[k] = lu_arg.solve(f_sum_);
239                }
240                // -------------------------------------------------------------------
241                // pack result_ into ty
242                for(size_t k = 0; k < n_order; k++)
243                {       for(size_t i = 0; i < nx_; i++)
244                                ty[ i * n_order + k ] = f_result_[k].data()[i];
245                }
246
247                // check if we are computing vy
248                if( vx.size() == 0 )
249                        return true;
250                // ------------------------------------------------------------------
251                // compute variable information for y; i.e., vy
252                // (note that the constant zero times a variable is a constant)
253                scalar zero(0.0);
254                assert( n_order == 1 );
255                for(size_t j = 0; j < nr_; j++)
256                {       for(size_t i = 0; i < nr_; i++)
257                        {       // initialize vy as false
258                                size_t index = i * nr_ + j;
259                                vy[index]    = false;
260                        }
261                        for(size_t i = 0; i < nr_; i++)
262                        {       for(size_t ell = 0; ell < nr_; ell++)
263                                {       // arg information
264                                        size_t index   = i * nr_ + ell;
265                                        bool var_arg  = vx[index];
266                                        // identity information
267                                        bool nz_eye   = ell == j;
268                                        // result information
269                                        index = i * nr_ + j;
270                                        vy[index] |= var_arg & nz_eye;
271                                }
272                        }
273                }
274                return true;
275        }
276/* %$$
277$subhead reverse$$
278$srccode%cpp% */
279        // reverse mode routine called by CppAD
280        virtual bool reverse(
281                // highest order Taylor coefficient that we are computing derivative of
282                size_t                     q ,
283                // forward mode Taylor coefficients for x variables
284                const CppAD::vector<double>&     tx ,
285                // forward mode Taylor coefficients for y variables
286                const CppAD::vector<double>&     ty ,
287                // upon return, derivative of G[ F[ {x_j^k} ] ] w.r.t {x_j^k}
288                CppAD::vector<double>&           px ,
289                // derivative of G[ {y_i^k} ] w.r.t. {y_i^k}
290                const CppAD::vector<double>&     py
291        )
292        {       size_t n_order = q + 1;
293                assert( nx_           == nr_ * nr_ );
294                assert( nx_ * n_order == tx.size() );
295                assert( nx_ * n_order == ty.size() );
296                assert( px.size()     == tx.size() );
297                assert( py.size()     == ty.size() );
298                //
299                // -------------------------------------------------------------------
300                // make sure f_arg_ is large enough
301                assert( f_arg_.size() == f_result_.size() );
302                // must have previous run forward with order >= n_order
303                assert( f_arg_.size() >= n_order );
304                // -------------------------------------------------------------------
305                // make sure r_arg_, r_result_ are large enough
306                assert( r_arg_.size() == r_result_.size() );
307                if( r_arg_.size() < n_order )
308                {       r_arg_.resize(n_order);
309                        r_result_.resize(n_order);
310                        //
311                        for(size_t k = 0; k < n_order; k++)
312                        {       r_arg_[k].resize(nr_, nr_);
313                                r_result_[k].resize(nr_, nr_);
314                        }
315                }
316                // -------------------------------------------------------------------
317                // unpack tx into f_arg_
318                for(size_t k = 0; k < n_order; k++)
319                {       // unpack arg values for this order
320                        for(size_t i = 0; i < nx_; i++)
321                                f_arg_[k].data()[i] = tx[ i * n_order + k ];
322                }
323                // -------------------------------------------------------------------
324                // unpack py into r_result_
325                for(size_t k = 0; k < n_order; k++)
326                {       for(size_t i = 0; i < nx_; i++)
327                                r_result_[k].data()[i] = py[ i * n_order + k ];
328                }
329                // -------------------------------------------------------------------
330                // initialize r_arg_ as zero
331                for(size_t k = 0; k < n_order; k++)
332                        r_arg_[k]   = matrix::Zero(nr_, nr_);
333                // -------------------------------------------------------------------
334                // matrix reverse mode calculation
335                //
336                for(size_t k1 = n_order; k1 > 1; k1--)
337                {       size_t k = k1 - 1;
338                        for(size_t ell = 1; ell <= k; ell++)
339                        {       r_arg_[ell]    -=
340                                        f_result_[0] * r_result_[k] * f_result_[k-ell].transpose();
341                                //
342                                r_result_[k-ell] -=
343                                        f_result_[0] * f_arg_[ell].transpose() * r_result_[k];
344                        }
345                        r_result_[0] +=
346                        r_result_[k] * f_result_[k].transpose() * f_arg_[0].transpose();
347                }
348                r_arg_[0] -=
349                f_result_[0].transpose() * r_result_[0] * f_result_[0].transpose();
350                // -------------------------------------------------------------------
351                // pack r_arg into px
352                for(size_t k = 0; k < n_order; k++)
353                {       for(size_t i = 0; i < nx_; i++)
354                                px[ i * n_order + k ] = r_arg_[k].data()[i];
355                }
356                //
357                return true;
358        }
359/* %$$
360$head End Class Definition$$
361$srccode%cpp% */
362}; // End of atomic_eigen_mat_inv class
363
364}  // END_EMPTY_NAMESPACE
365/* %$$
366$$ $comment end nospell$$
367$end
368*/
369
370
371# endif
Note: See TracBrowser for help on using the repository browser.