source: trunk/cppad/example/eigen_cholesky.hpp @ 3832

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

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

commit d9bb35cfcbbf4a69640d4e20aa59422061e65891
Author: Brad Bell <bradbell@…>
Date: Tue Sep 27 07:50:07 2016 -0700

  1. Advance version to cppad-20160927.
  2. Change cholesky -> cholesky_theory.


eigen_cholesky.hpp: make computation of M_k agree with theory.

File size: 12.0 KB
Line 
1// $Id$
2# ifndef CPPAD_EXAMPLE_EIGEN_CHOLESKY_HPP
3# define CPPAD_EXAMPLE_EIGEN_CHOLESKY_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_cholesky.hpp$$
18$spell
19        Eigen
20        Taylor
21        Cholesky
22        op
23$$
24
25$section Atomic Eigen Cholesky Factorization Class$$
26
27$head Purpose$$
28Construct an atomic operation that computes a lower triangular matrix
29$latex L $$ such that $latex L L^\R{T} = A$$
30for any positive integer $latex p$$
31and symmetric positive definite matrix $latex A \in \B{R}^{p \times p}$$.
32
33$head Start Class Definition$$
34$srccode%cpp% */
35# include <cppad/cppad.hpp>
36# include <Eigen/Dense>
37
38
39/* %$$
40$head Public$$
41
42$subhead Types$$
43$srccode%cpp% */
44namespace { // BEGIN_EMPTY_NAMESPACE
45
46template <class Base>
47class atomic_eigen_cholesky : public CppAD::atomic_base<Base> {
48public:
49        // -----------------------------------------------------------
50        // type of elements during calculation of derivatives
51        typedef Base              scalar;
52        // type of elements during taping
53        typedef CppAD::AD<scalar> ad_scalar;
54        //
55        // type of matrix during calculation of derivatives
56        typedef Eigen::Matrix<
57                scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>        matrix;
58        // type of matrix during taping
59        typedef Eigen::Matrix<
60                ad_scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > ad_matrix;
61        //
62        // lower triangular scalar matrix
63        typedef Eigen::TriangularView<matrix, Eigen::Lower>             lower_view;
64/* %$$
65$subhead Constructor$$
66$srccode%cpp% */
67        // constructor
68        atomic_eigen_cholesky(void) : CppAD::atomic_base<Base>(
69                "atom_eigen_cholesky"                             ,
70                CppAD::atomic_base<Base>::set_sparsity_enum
71        )
72        { }
73/* %$$
74$subhead op$$
75$srccode%cpp% */
76        // use atomic operation to invert an AD matrix
77        ad_matrix op(const ad_matrix& arg)
78        {       size_t nr = size_t( arg.rows() );
79                size_t ny = ( (nr + 1 ) * nr ) / 2;
80                size_t nx = 1 + ny;
81                assert( nr == size_t( arg.cols() ) );
82                // -------------------------------------------------------------------
83                // packed version of arg
84                CPPAD_TESTVECTOR(ad_scalar) packed_arg(nx);
85                size_t index = 0;
86                packed_arg[index++] = ad_scalar( nr );
87                // lower triangle of symmetric matrix A
88                for(size_t i = 0; i < nr; i++)
89                {       for(size_t j = 0; j <= i; j++)
90                                packed_arg[index++] = arg(i, j);
91                }
92                assert( index == nx );
93                // -------------------------------------------------------------------
94                // packed version of result = arg^{-1}.
95                // This is an atomic_base function call that CppAD uses to
96                // store the atomic operation on the tape.
97                CPPAD_TESTVECTOR(ad_scalar) packed_result(ny);
98                (*this)(packed_arg, packed_result);
99                // -------------------------------------------------------------------
100                // unpack result matrix L
101                ad_matrix result = ad_matrix::Zero(nr, nr);
102                index = 0;
103                for(size_t i = 0; i < nr; i++)
104                {       for(size_t j = 0; j <= i; j++)
105                                result(i, j) = packed_result[index++];
106                }
107                return result;
108        }
109        /* %$$
110$head Private$$
111
112$subhead Variables$$
113$srccode%cpp% */
114private:
115        // -------------------------------------------------------------
116        // one forward mode vector of matrices for argument and result
117        CppAD::vector<matrix> f_arg_, f_result_;
118        // one reverse mode vector of matrices for argument and result
119        CppAD::vector<matrix> r_arg_, r_result_;
120        // -------------------------------------------------------------
121/* %$$
122$subhead forward$$
123$srccode%cpp% */
124        // forward mode routine called by CppAD
125        virtual bool forward(
126                // lowest order Taylor coefficient we are evaluating
127                size_t                          p ,
128                // highest order Taylor coefficient we are evaluating
129                size_t                          q ,
130                // which components of x are variables
131                const CppAD::vector<bool>&      vx ,
132                // which components of y are variables
133                CppAD::vector<bool>&            vy ,
134                // tx [ j * (q+1) + k ] is x_j^k
135                const CppAD::vector<scalar>&    tx ,
136                // ty [ i * (q+1) + k ] is y_i^k
137                CppAD::vector<scalar>&          ty
138        )
139        {       size_t n_order = q + 1;
140                size_t nr      = size_t( CppAD::Integer( tx[ 0 * n_order + 0 ] ) );
141                size_t ny      = ((nr + 1) * nr) / 2;
142                size_t nx      = 1 + ny;
143                assert( vx.size() == 0 || nx == vx.size() );
144                assert( vx.size() == 0 || ny == vy.size() );
145                assert( nx * n_order == tx.size() );
146                assert( ny * n_order == ty.size() );
147                //
148                // -------------------------------------------------------------------
149                // make sure f_arg_ and f_result_ are large enough
150                assert( f_arg_.size() == f_result_.size() );
151                if( f_arg_.size() < n_order )
152                {       f_arg_.resize(n_order);
153                        f_result_.resize(n_order);
154                        //
155                        for(size_t k = 0; k < n_order; k++)
156                        {       f_arg_[k].resize(nr, nr);
157                                f_result_[k].resize(nr, nr);
158                        }
159                }
160                // -------------------------------------------------------------------
161                // unpack tx into f_arg_
162                for(size_t k = 0; k < n_order; k++)
163                {       size_t index = 1;
164                        // unpack arg values for this order
165                        for(size_t i = 0; i < nr; i++)
166                        {       for(size_t j = 0; j <= i; j++)
167                                {       f_arg_[k](i, j) = tx[ index * n_order + k ];
168                                        f_arg_[k](j, i) = f_arg_[k](i, j);
169                                        index++;
170                                }
171                        }
172                }
173                // -------------------------------------------------------------------
174                // result for each order
175                // (we could avoid recalculting f_result_[k] for k=0,...,p-1)
176                //
177                Eigen::LLT<matrix> cholesky(f_arg_[0]);
178                f_result_[0]   = cholesky.matrixL();
179                lower_view L_0 = f_result_[0].template triangularView<Eigen::Lower>();
180                for(size_t k = 1; k < n_order; k++)
181                {       // initialize sum as A_k
182                        matrix f_sum = f_arg_[k];
183                        // compute A_k - B_k
184                        for(size_t ell = 1; ell < k; ell++)
185                                f_sum -= f_result_[ell] * f_result_[k-ell].transpose();
186                        // compute L_0^{-1} * (A_k - B_k) * L_0^{-T}
187                        matrix temp = L_0.template solve<Eigen::OnTheLeft>(f_sum);
188                        temp   = L_0.transpose().template solve<Eigen::OnTheRight>(temp);
189                        // divide the diagonal by 2
190                        for(size_t i = 0; i < nr; i++)
191                                temp(i, i) /= scalar(2.0);
192                        // L_k = L_0 * low[ L_0^{-1} * (A_k - B_k) * L_0^{-T} ]
193                        lower_view view = temp.template triangularView<Eigen::Lower>();
194                        f_result_[k] = f_result_[0] * view;
195                }
196                // -------------------------------------------------------------------
197                // pack result_ into ty
198                for(size_t k = 0; k < n_order; k++)
199                {       size_t index = 0;
200                        for(size_t i = 0; i < nr; i++)
201                        {       for(size_t j = 0; j <= i; j++)
202                                {       ty[ index * n_order + k ] = f_result_[k](i, j);
203                                        index++;
204                                }
205                        }
206                }
207                // -------------------------------------------------------------------
208                // check if we are computing vy
209                if( vx.size() == 0 )
210                        return true;
211                // ------------------------------------------------------------------
212                // This is a very dumb algorithm that over estimates which
213                // elements of the inverse are variables (which is not efficient).
214                bool var = false;
215                for(size_t i = 0; i < ny; i++)
216                        var |= vx[1 + i];
217                for(size_t i = 0; i < ny; i++)
218                        vy[i] = var;
219                //
220                return true;
221        }
222/* %$$
223$subhead reverse$$
224$srccode%cpp% */
225        // reverse mode routine called by CppAD
226        virtual bool reverse(
227                // highest order Taylor coefficient that we are computing derivative of
228                size_t                     q ,
229                // forward mode Taylor coefficients for x variables
230                const CppAD::vector<double>&     tx ,
231                // forward mode Taylor coefficients for y variables
232                const CppAD::vector<double>&     ty ,
233                // upon return, derivative of G[ F[ {x_j^k} ] ] w.r.t {x_j^k}
234                CppAD::vector<double>&           px ,
235                // derivative of G[ {y_i^k} ] w.r.t. {y_i^k}
236                const CppAD::vector<double>&     py
237        )
238        {       size_t n_order = q + 1;
239                size_t nr = size_t( CppAD::Integer( tx[ 0 * n_order + 0 ] ) );
240                size_t ny = ( (nr + 1 ) * nr ) / 2;
241                size_t nx = 1 + ny;
242                //
243                assert( nx * n_order == tx.size() );
244                assert( ny * n_order == ty.size() );
245                assert( px.size()    == tx.size() );
246                assert( py.size()    == ty.size() );
247                // -------------------------------------------------------------------
248                // make sure f_arg_ is large enough
249                assert( f_arg_.size() == f_result_.size() );
250                // must have previous run forward with order >= n_order
251                assert( f_arg_.size() >= n_order );
252                // -------------------------------------------------------------------
253                // make sure r_arg_, r_result_ are large enough
254                assert( r_arg_.size() == r_result_.size() );
255                if( r_arg_.size() < n_order )
256                {       r_arg_.resize(n_order);
257                        r_result_.resize(n_order);
258                        //
259                        for(size_t k = 0; k < n_order; k++)
260                        {       r_arg_[k].resize(nr, nr);
261                                r_result_[k].resize(nr, nr);
262                        }
263                }
264                // -------------------------------------------------------------------
265                // unpack tx into f_arg_
266                for(size_t k = 0; k < n_order; k++)
267                {       size_t index = 1;
268                        // unpack arg values for this order
269                        for(size_t i = 0; i < nr; i++)
270                        {       for(size_t j = 0; j <= i; j++)
271                                {       f_arg_[k](i, j) = tx[ index * n_order + k ];
272                                        f_arg_[k](j, i) = f_arg_[k](i, j);
273                                        index++;
274                                }
275                        }
276                }
277                // -------------------------------------------------------------------
278                // unpack py into r_result_
279                for(size_t k = 0; k < n_order; k++)
280                {       r_result_[k] = matrix::Zero(nr, nr);
281                        size_t index = 0;
282                        for(size_t i = 0; i < nr; i++)
283                        {       for(size_t j = 0; j <= i; j++)
284                                {       r_result_[k](i, j) = py[ index * n_order + k ];
285                                        index++;
286                                }
287                        }
288                }
289                // -------------------------------------------------------------------
290                // initialize r_arg_ as zero
291                for(size_t k = 0; k < n_order; k++)
292                        r_arg_[k]   = matrix::Zero(nr, nr);
293                // -------------------------------------------------------------------
294                // matrix reverse mode calculation
295                lower_view L_0 = f_result_[0].template triangularView<Eigen::Lower>();
296                //
297                for(size_t k1 = n_order; k1 > 1; k1--)
298                {       size_t k = k1 - 1;
299                        //
300                        // L_0^T * bar{L}_k
301                        matrix tmp1 = L_0.transpose() * r_result_[k];
302                        //
303                        //low[ L_0^T * bar{L}_k ]
304                        for(size_t i = 0; i < nr; i++)
305                                tmp1(i, i) /= scalar(2.0);
306                        matrix tmp2 = tmp1.template triangularView<Eigen::Lower>();
307                        //
308                        // L_0^{-T} low[ L_0^T * bar{L}_k ]
309                        tmp1 = L_0.transpose().template solve<Eigen::OnTheLeft>( tmp2 );
310                        //
311                        // M_k = L_0^{-T} * low[ L_0^T * bar{L}_k ]^{T} L_0^{-1}
312                        matrix M_k = L_0.transpose().template
313                                solve<Eigen::OnTheLeft>( tmp1.transpose() );
314                        //
315                        // remove L_k and compute bar{B}_k
316                        matrix barB_k = scalar(0.5) * ( M_k + M_k.transpose() );
317                        r_arg_[k]    += barB_k;
318                        barB_k        = scalar(-1.0) * barB_k;
319                        //
320                        // 2.0 * lower( bar{B}_k L_k )
321                        matrix temp = scalar(2.0) * barB_k * f_result_[k];
322                        temp        = temp.template triangularView<Eigen::Lower>();
323                        //
324                        // remove C_k
325                        r_result_[0] += temp;
326                        //
327                        // remove B_k
328                        for(size_t ell = 1; ell < k; ell++)
329                        {       // bar{L}_ell = 2 * lower( \bar{B}_k * L_{k-ell} )
330                                temp = scalar(2.0) * barB_k * f_result_[k-ell];
331                                r_result_[ell] += temp.template triangularView<Eigen::Lower>();
332                        }
333                }
334                // M_0 = L_0^{-T} * low[ L_0^T * bar{L}_0 ]^{T} L_0^{-1}
335                matrix M_0 = L_0.transpose() * r_result_[0];
336                for(size_t i = 0; i < nr; i++)
337                        M_0(i, i) /= scalar(2.0);
338                M_0 = M_0.template triangularView<Eigen::Lower>();
339                M_0 = L_0.template solve<Eigen::OnTheRight>( M_0 );
340                M_0 = L_0.transpose().template solve<Eigen::OnTheLeft>( M_0 );
341                // remove L_0
342                r_arg_[0] += scalar(0.5) * ( M_0 + M_0.transpose() );
343                // -------------------------------------------------------------------
344                // pack r_arg into px
345                // note that only the lower triangle of barA_k is stored in px
346                for(size_t k = 0; k < n_order; k++)
347                {       size_t index = 0;
348                        px[ index * n_order + k ] = 0.0;
349                        index++;
350                        for(size_t i = 0; i < nr; i++)
351                        {       for(size_t j = 0; j < i; j++)
352                                {       px[ index * n_order + k ] = 2.0 * r_arg_[k](i, j);
353                                        index++;
354                                }
355                                px[ index * n_order + k] = r_arg_[k](i, i);
356                                index++;
357                        }
358                }
359                // -------------------------------------------------------------------
360                return true;
361        }
362/* %$$
363$head End Class Definition$$
364$srccode%cpp% */
365}; // End of atomic_eigen_cholesky class
366
367}  // END_EMPTY_NAMESPACE
368/* %$$
369$end
370*/
371
372
373# endif
Note: See TracBrowser for help on using the repository browser.