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