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

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

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

commit 33d1f5b8d837d83c8baf656045d85d8b45f0b297
Author: Brad Bell <bradbell@…>
Date: Sun Mar 27 04:01:14 2016 -0700

  1. Test non-symmetric reverse mode in eigen_mat_inv (and fix this example).
  2. Advance version to cppad-20160327.

commit 238e214ff36ca835efca615c609576dc8bf5038d
Author: Brad Bell <bradbell@…>
Date: Sat Mar 26 18:10:59 2016 -0700

eigen_mat_inv.hpp: use inverse matrix (since we are computing it).
eigen_mat_mul.hpp: comment seperator.
eigen_mat_inv.cpp: matrix during recording is non-singular.
eigen_mat_inv.cpp: test first and second order forward.

commit fc918b0476cc8ea66abdf2904a71fc93472d279d
Author: Brad Bell <bradbell@…>
Date: Sat Mar 26 16:26:51 2016 -0700

Add test_more/eigen_mat_inv.cpp.


eigen_mat_inv.hpp: fix calculation of vy.
eigen_mat_inv.cpp: fix section title.
CMakeLists.txt: no longer necessary to have special library for eigen tests.
eigen_mat_inv.cpp: Test with non-symmetric matrix.

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