CppAD: A C++ Algorithmic Differentiation Package
20130918
|
00001 /* $Id$ */ 00002 # ifndef CPPAD_BENDER_QUAD_INCLUDED 00003 # define CPPAD_BENDER_QUAD_INCLUDED 00004 00005 /* -------------------------------------------------------------------------- 00006 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-13 Bradley M. Bell 00007 00008 CppAD is distributed under multiple licenses. This distribution is under 00009 the terms of the 00010 Eclipse Public License Version 1.0. 00011 00012 A copy of this license is included in the COPYING file of this distribution. 00013 Please visit http://www.coin-or.org/CppAD/ for information on other licenses. 00014 -------------------------------------------------------------------------- */ 00015 /* 00016 $begin BenderQuad$$ 00017 $spell 00018 cppad.hpp 00019 BAvector 00020 gx 00021 gxx 00022 CppAD 00023 Fy 00024 dy 00025 Jacobian 00026 ADvector 00027 const 00028 dg 00029 ddg 00030 $$ 00031 00032 00033 $index Hessian, Bender$$ 00034 $index Jacobian, Bender$$ 00035 $index BenderQuad$$ 00036 $section Computing Jacobian and Hessian of Bender's Reduced Objective$$ 00037 00038 $head Syntax$$ 00039 $codei% 00040 # include <cppad/cppad.hpp> 00041 BenderQuad(%x%, %y%, %fun%, %g%, %gx%, %gxx%)%$$ 00042 00043 $head See Also$$ 00044 $cref opt_val_hes$$ 00045 00046 $head Problem$$ 00047 The type $cref/ADvector/BenderQuad/ADvector/$$ cannot be determined 00048 form the arguments above 00049 (currently the type $icode ADvector$$ must be 00050 $codei%CPPAD_TESTVECTOR(%Base%)%$$.) 00051 This will be corrected in the future by requiring $icode Fun$$ 00052 to define $icode%Fun%::vector_type%$$ which will specify the 00053 type $icode ADvector$$. 00054 00055 $head Purpose$$ 00056 We are given the optimization problem 00057 $latex \[ 00058 \begin{array}{rcl} 00059 {\rm minimize} & F(x, y) & {\rm w.r.t.} \; (x, y) \in \B{R}^n \times \B{R}^m 00060 \end{array} 00061 \] $$ 00062 that is convex with respect to $latex y$$. 00063 In addition, we are given a set of equations $latex H(x, y)$$ 00064 such that 00065 $latex \[ 00066 H[ x , Y(x) ] = 0 \;\; \Rightarrow \;\; F_y [ x , Y(x) ] = 0 00067 \] $$ 00068 (In fact, it is often the case that $latex H(x, y) = F_y (x, y)$$.) 00069 Furthermore, it is easy to calculate a Newton step for these equations; i.e., 00070 $latex \[ 00071 dy = - [ \partial_y H(x, y)]^{-1} H(x, y) 00072 \] $$ 00073 The purpose of this routine is to compute the 00074 value, Jacobian, and Hessian of the reduced objective function 00075 $latex \[ 00076 G(x) = F[ x , Y(x) ] 00077 \] $$ 00078 Note that if only the value and Jacobian are needed, they can be 00079 computed more quickly using the relations 00080 $latex \[ 00081 G^{(1)} (x) = \partial_x F [x, Y(x) ] 00082 \] $$ 00083 00084 $head x$$ 00085 The $code BenderQuad$$ argument $icode x$$ has prototype 00086 $codei% 00087 const %BAvector% &%x% 00088 %$$ 00089 (see $cref/BAvector/BenderQuad/BAvector/$$ below) 00090 and its size must be equal to $icode n$$. 00091 It specifies the point at which we evaluating 00092 the reduced objective function and its derivatives. 00093 00094 00095 $head y$$ 00096 The $code BenderQuad$$ argument $icode y$$ has prototype 00097 $codei% 00098 const %BAvector% &%y% 00099 %$$ 00100 and its size must be equal to $icode m$$. 00101 It must be equal to $latex Y(x)$$; i.e., 00102 it must solve the problem in $latex y$$ for this given value of $latex x$$ 00103 $latex \[ 00104 \begin{array}{rcl} 00105 {\rm minimize} & F(x, y) & {\rm w.r.t.} \; y \in \B{R}^m 00106 \end{array} 00107 \] $$ 00108 00109 $head fun$$ 00110 The $code BenderQuad$$ object $icode fun$$ 00111 must support the member functions listed below. 00112 The $codei%AD<%Base%>%$$ arguments will be variables for 00113 a tape created by a call to $cref Independent$$ from $code BenderQuad$$ 00114 (hence they can not be combined with variables corresponding to a 00115 different tape). 00116 00117 $subhead fun.f$$ 00118 The $code BenderQuad$$ argument $icode fun$$ supports the syntax 00119 $codei% 00120 %f% = %fun%.f(%x%, %y%) 00121 %$$ 00122 The $icode%fun%.f%$$ argument $icode x$$ has prototype 00123 $codei% 00124 const %ADvector% &%x% 00125 %$$ 00126 (see $cref/ADvector/BenderQuad/ADvector/$$ below) 00127 and its size must be equal to $icode n$$. 00128 The $icode%fun%.f%$$ argument $icode y$$ has prototype 00129 $codei% 00130 const %ADvector% &%y% 00131 %$$ 00132 and its size must be equal to $icode m$$. 00133 The $icode%fun%.f%$$ result $icode f$$ has prototype 00134 $codei% 00135 %ADvector% %f% 00136 %$$ 00137 and its size must be equal to one. 00138 The value of $icode f$$ is 00139 $latex \[ 00140 f = F(x, y) 00141 \] $$. 00142 00143 $subhead fun.h$$ 00144 The $code BenderQuad$$ argument $icode fun$$ supports the syntax 00145 $codei% 00146 %h% = %fun%.h(%x%, %y%) 00147 %$$ 00148 The $icode%fun%.h%$$ argument $icode x$$ has prototype 00149 $codei% 00150 const %ADvector% &%x% 00151 %$$ 00152 and its size must be equal to $icode n$$. 00153 The $icode%fun%.h%$$ argument $icode y$$ has prototype 00154 $codei% 00155 const %BAvector% &%y% 00156 %$$ 00157 and its size must be equal to $icode m$$. 00158 The $icode%fun%.h%$$ result $icode h$$ has prototype 00159 $codei% 00160 %ADvector% %h% 00161 %$$ 00162 and its size must be equal to $icode m$$. 00163 The value of $icode h$$ is 00164 $latex \[ 00165 h = H(x, y) 00166 \] $$. 00167 00168 $subhead fun.dy$$ 00169 The $code BenderQuad$$ argument $icode fun$$ supports the syntax 00170 $codei% 00171 %dy% = %fun%.dy(%x%, %y%, %h%) 00172 00173 %x% 00174 %$$ 00175 The $icode%fun%.dy%$$ argument $icode x$$ has prototype 00176 $codei% 00177 const %BAvector% &%x% 00178 %$$ 00179 and its size must be equal to $icode n$$. 00180 Its value will be exactly equal to the $code BenderQuad$$ argument 00181 $icode x$$ and values depending on it can be stored as private objects 00182 in $icode f$$ and need not be recalculated. 00183 $codei% 00184 00185 %y% 00186 %$$ 00187 The $icode%fun%.dy%$$ argument $icode y$$ has prototype 00188 $codei% 00189 const %BAvector% &%y% 00190 %$$ 00191 and its size must be equal to $icode m$$. 00192 Its value will be exactly equal to the $code BenderQuad$$ argument 00193 $icode y$$ and values depending on it can be stored as private objects 00194 in $icode f$$ and need not be recalculated. 00195 $codei% 00196 00197 %h% 00198 %$$ 00199 The $icode%fun%.dy%$$ argument $icode h$$ has prototype 00200 $codei% 00201 const %ADvector% &%h% 00202 %$$ 00203 and its size must be equal to $icode m$$. 00204 $codei% 00205 00206 %dy% 00207 %$$ 00208 The $icode%fun%.dy%$$ result $icode dy$$ has prototype 00209 $codei% 00210 %ADvector% %dy% 00211 %$$ 00212 and its size must be equal to $icode m$$. 00213 The return value $icode dy$$ is given by 00214 $latex \[ 00215 dy = - [ \partial_y H (x , y) ]^{-1} h 00216 \] $$ 00217 Note that if $icode h$$ is equal to $latex H(x, y)$$, 00218 $latex dy$$ is the Newton step for finding a zero 00219 of $latex H(x, y)$$ with respect to $latex y$$; 00220 i.e., 00221 $latex y + dy$$ is an approximate solution for the equation 00222 $latex H (x, y + dy) = 0$$. 00223 00224 $head g$$ 00225 The argument $icode g$$ has prototype 00226 $codei% 00227 %BAvector% &%g% 00228 %$$ 00229 and has size one. 00230 The input value of its element does not matter. 00231 On output, 00232 it contains the value of $latex G (x)$$; i.e., 00233 $latex \[ 00234 g[0] = G (x) 00235 \] $$ 00236 00237 $head gx$$ 00238 The argument $icode gx$$ has prototype 00239 $codei% 00240 %BAvector% &%gx% 00241 %$$ 00242 and has size $latex n $$. 00243 The input values of its elements do not matter. 00244 On output, 00245 it contains the Jacobian of $latex G (x)$$; i.e., 00246 for $latex j = 0 , \ldots , n-1$$, 00247 $latex \[ 00248 gx[ j ] = G^{(1)} (x)_j 00249 \] $$ 00250 00251 $head gxx$$ 00252 The argument $icode gx$$ has prototype 00253 $codei% 00254 %BAvector% &%gxx% 00255 %$$ 00256 and has size $latex n \times n$$. 00257 The input values of its elements do not matter. 00258 On output, 00259 it contains the Hessian of $latex G (x)$$; i.e., 00260 for $latex i = 0 , \ldots , n-1$$, and 00261 $latex j = 0 , \ldots , n-1$$, 00262 $latex \[ 00263 gxx[ i * n + j ] = G^{(2)} (x)_{i,j} 00264 \] $$ 00265 00266 $head BAvector$$ 00267 The type $icode BAvector$$ must be a 00268 $cref SimpleVector$$ class. 00269 We use $icode Base$$ to refer to the type of the elements of 00270 $icode BAvector$$; i.e., 00271 $codei% 00272 %BAvector%::value_type 00273 %$$ 00274 00275 $head ADvector$$ 00276 The type $icode ADvector$$ must be a 00277 $cref SimpleVector$$ class with elements of type 00278 $codei%AD<%Base%>%$$; i.e., 00279 $codei% 00280 %ADvector%::value_type 00281 %$$ 00282 must be the same type as 00283 $codei% 00284 AD< %BAvector%::value_type > 00285 %$$. 00286 00287 00288 $head Example$$ 00289 $children% 00290 example/bender_quad.cpp 00291 %$$ 00292 The file 00293 $cref bender_quad.cpp$$ 00294 contains an example and test of this operation. 00295 It returns true if it succeeds and false otherwise. 00296 00297 00298 $end 00299 ----------------------------------------------------------------------------- 00300 */ 00301 00302 namespace CppAD { // BEGIN CppAD namespace 00303 00304 template <class BAvector, class Fun> 00305 void BenderQuad( 00306 const BAvector &x , 00307 const BAvector &y , 00308 Fun fun , 00309 BAvector &g , 00310 BAvector &gx , 00311 BAvector &gxx ) 00312 { // determine the base type 00313 typedef typename BAvector::value_type Base; 00314 00315 // check that BAvector is a SimpleVector class 00316 CheckSimpleVector<Base, BAvector>(); 00317 00318 // declare the ADvector type 00319 typedef CPPAD_TESTVECTOR(AD<Base>) ADvector; 00320 00321 // size of the x and y spaces 00322 size_t n = size_t(x.size()); 00323 size_t m = size_t(y.size()); 00324 00325 // check the size of gx and gxx 00326 CPPAD_ASSERT_KNOWN( 00327 g.size() == 1, 00328 "BenderQuad: size of the vector g is not equal to 1" 00329 ); 00330 CPPAD_ASSERT_KNOWN( 00331 size_t(gx.size()) == n, 00332 "BenderQuad: size of the vector gx is not equal to n" 00333 ); 00334 CPPAD_ASSERT_KNOWN( 00335 size_t(gxx.size()) == n * n, 00336 "BenderQuad: size of the vector gxx is not equal to n * n" 00337 ); 00338 00339 // some temporary indices 00340 size_t i, j; 00341 00342 // variable versions x 00343 ADvector vx(n); 00344 for(j = 0; j < n; j++) 00345 vx[j] = x[j]; 00346 00347 // declare the independent variables 00348 Independent(vx); 00349 00350 // evaluate h = H(x, y) 00351 ADvector h(m); 00352 h = fun.h(vx, y); 00353 00354 // evaluate dy (x) = Newton step as a function of x through h only 00355 ADvector dy(m); 00356 dy = fun.dy(x, y, h); 00357 00358 // variable version of y 00359 ADvector vy(m); 00360 for(j = 0; j < m; j++) 00361 vy[j] = y[j] + dy[j]; 00362 00363 // evaluate G~ (x) = F [ x , y + dy(x) ] 00364 ADvector gtilde(1); 00365 gtilde = fun.f(vx, vy); 00366 00367 // AD function object that corresponds to G~ (x) 00368 // We will make heavy use of this tape, so optimize it 00369 ADFun<Base> Gtilde; 00370 Gtilde.Dependent(vx, gtilde); 00371 Gtilde.optimize(); 00372 00373 // value of G(x) 00374 g = Gtilde.Forward(0, x); 00375 00376 // initial forward direction vector as zero 00377 BAvector dx(n); 00378 for(j = 0; j < n; j++) 00379 dx[j] = Base(0); 00380 00381 // weight, first and second order derivative values 00382 BAvector dg(1), w(1), ddw(2 * n); 00383 w[0] = 1.; 00384 00385 00386 // Jacobian and Hessian of G(x) is equal Jacobian and Hessian of Gtilde 00387 for(j = 0; j < n; j++) 00388 { // compute partials in x[j] direction 00389 dx[j] = Base(1); 00390 dg = Gtilde.Forward(1, dx); 00391 gx[j] = dg[0]; 00392 00393 // restore the dx vector to zero 00394 dx[j] = Base(0); 00395 00396 // compute second partials w.r.t x[j] and x[l] for l = 1, n 00397 ddw = Gtilde.Reverse(2, w); 00398 for(i = 0; i < n; i++) 00399 gxx[ i * n + j ] = ddw[ i * 2 + 1 ]; 00400 } 00401 00402 return; 00403 } 00404 00405 } // END CppAD namespace 00406 00407 # endif