CppAD: A C++ Algorithmic Differentiation Package  20171217
Go to the documentation of this file.
3
4 /* --------------------------------------------------------------------------
6
8 the terms of the
9  Eclipse Public License Version 1.0.
10
11 A copy of this license is included in the COPYING file of this distribution.
13 -------------------------------------------------------------------------- */
14 /*
15 $begin BenderQuad$$16 spell 17 cppad.hpp 18 BAvector 19 gx 20 gxx 21 CppAD 22 Fy 23 dy 24 Jacobian 25 ADvector 26 const 27 dg 28 ddg 29$$ 30 31 32$section Computing Jacobian and Hessian of Bender's Reduced Objective$$33 mindex BenderQuad$$
34
35 $head Syntax$$36 codei% 37 # include <cppad/cppad.hpp> 38 BenderQuad(%x%, %y%, %fun%, %g%, %gx%, %gxx%)%$$ 39 40$head See Also$$41 cref opt_val_hes$$
42
43 $head Problem$$44 The type cref/ADvector/BenderQuad/ADvector/$$ cannot be determined 45 form the arguments above 46 (currently the type$icode ADvector$$must be 47 codei%CPPAD_TESTVECTOR(%Base%)%$$.)
48 This will be corrected in the future by requiring $icode Fun$$49 to define icode%Fun%::vector_type%$$ which will specify the 50 type$icode ADvector$$. 51 52 head Purpose$$
53 We are given the optimization problem
54 $latex $55 \begin{array}{rcl} 56 {\rm minimize} & F(x, y) & {\rm w.r.t.} \; (x, y) \in \B{R}^n \times \B{R}^m 57 \end{array} 58$ $$59 that is convex with respect to latex y$$. 60 In addition, we are given a set of equations$latex H(x, y)$$61 such that 62 latex $63 H[ x , Y(x) ] = 0 \;\; \Rightarrow \;\; F_y [ x , Y(x) ] = 0 64$$$
65 (In fact, it is often the case that $latex H(x, y) = F_y (x, y)$$.) 66 Furthermore, it is easy to calculate a Newton step for these equations; i.e., 67 latex $68 dy = - [ \partial_y H(x, y)]^{-1} H(x, y) 69$$$ 70 The purpose of this routine is to compute the 71 value, Jacobian, and Hessian of the reduced objective function 72$latex $73 G(x) = F[ x , Y(x) ] 74$ $$75 Note that if only the value and Jacobian are needed, they can be 76 computed more quickly using the relations 77 latex $78 G^{(1)} (x) = \partial_x F [x, Y(x) ] 79$$$
80
81 $head x$$82 The code BenderQuad$$ argument$icode x$$has prototype 83 codei% 84 const %BAvector% &%x% 85 %$$
86 (see $cref/BAvector/BenderQuad/BAvector/$$below) 87 and its size must be equal to icode n$$. 88 It specifies the point at which we evaluating 89 the reduced objective function and its derivatives. 90 91 92$head y$$93 The code BenderQuad$$ argument $icode y$$has prototype 94 codei% 95 const %BAvector% &%y% 96 %$$ 97 and its size must be equal to$icode m$$. 98 It must be equal to latex Y(x)$$; i.e.,
99 it must solve the problem in $latex y$$for this given value of latex x$$ 100$latex $101 \begin{array}{rcl} 102 {\rm minimize} & F(x, y) & {\rm w.r.t.} \; y \in \B{R}^m 103 \end{array} 104$ $$105 106 head fun$$
107 The $code BenderQuad$$object icode fun$$ 108 must support the member functions listed below. 109 The$codei%AD<%Base%>%$$arguments will be variables for 110 a tape created by a call to cref Independent$$ from $code BenderQuad$$111 (hence they can not be combined with variables corresponding to a 112 different tape). 113 114 subhead fun.f$$ 115 The$code BenderQuad$$argument icode fun$$ supports the syntax
116 $codei% 117 %f% = %fun%.f(%x%, %y%) 118 %$$119 The icode%fun%.f%$$ argument$icode x$$has prototype 120 codei% 121 const %ADvector% &%x% 122 %$$
123 (see $cref/ADvector/BenderQuad/ADvector/$$below) 124 and its size must be equal to icode n$$. 125 The$icode%fun%.f%$$argument icode y$$ has prototype
126 $codei% 127 const %ADvector% &%y% 128 %$$129 and its size must be equal to icode m$$. 130 The$icode%fun%.f%$$result icode f$$ has prototype
131 $codei% 132 %ADvector% %f% 133 %$$134 and its size must be equal to one. 135 The value of icode f$$ is 136$latex $137 f = F(x, y) 138$ $$. 139 140 subhead fun.h$$
141 The $code BenderQuad$$argument icode fun$$ supports the syntax 142$codei%
143  %h% = %fun%.h(%x%, %y%)
144 %$$145 The icode%fun%.h%$$ argument $icode x$$has prototype 146 codei% 147 const %ADvector% &%x% 148 %$$ 149 and its size must be equal to$icode n$$. 150 The icode%fun%.h%$$ argument $icode y$$has prototype 151 codei% 152 const %BAvector% &%y% 153 %$$ 154 and its size must be equal to$icode m$$. 155 The icode%fun%.h%$$ result $icode h$$has prototype 156 codei% 157 %ADvector% %h% 158 %$$ 159 and its size must be equal to$icode m$$. 160 The value of icode h$$ is
161 $latex $162 h = H(x, y) 163$ $$. 164 165 subhead fun.dy$$ 166 The$code BenderQuad$$argument icode fun$$ supports the syntax
167 $codei% 168 %dy% = %fun%.dy(%x%, %y%, %h%) 169 170 %x% 171 %$$172 The icode%fun%.dy%$$ argument$icode x$$has prototype 173 codei% 174 const %BAvector% &%x% 175 %$$
176 and its size must be equal to $icode n$$. 177 Its value will be exactly equal to the code BenderQuad$$ argument 178$icode x$$and values depending on it can be stored as private objects 179 in icode f$$ and need not be recalculated.
180 $codei% 181 182 %y% 183 %$$184 The icode%fun%.dy%$$ argument$icode y$$has prototype 185 codei% 186 const %BAvector% &%y% 187 %$$
188 and its size must be equal to $icode m$$. 189 Its value will be exactly equal to the code BenderQuad$$ argument 190$icode y$$and values depending on it can be stored as private objects 191 in icode f$$ and need not be recalculated.
192 $codei% 193 194 %h% 195 %$$196 The icode%fun%.dy%$$ argument$icode h$$has prototype 197 codei% 198 const %ADvector% &%h% 199 %$$
200 and its size must be equal to $icode m$$. 201 codei% 202 203 %dy% 204 %$$ 205 The$icode%fun%.dy%$$result icode dy$$ has prototype
206 $codei% 207 %ADvector% %dy% 208 %$$209 and its size must be equal to icode m$$. 210 The return value$icode dy$$is given by 211 latex $212 dy = - [ \partial_y H (x , y) ]^{-1} h 213$$$
214 Note that if $icode h$$is equal to latex H(x, y)$$, 215$latex dy$$is the Newton step for finding a zero 216 of latex H(x, y)$$ with respect to $latex y$$; 217 i.e., 218 latex y + dy$$ is an approximate solution for the equation 219$latex H (x, y + dy) = 0$$. 220 221 head g$$
222 The argument $icode g$$has prototype 223 codei% 224 %BAvector% &%g% 225 %$$ 226 and has size one. 227 The input value of its element does not matter. 228 On output, 229 it contains the value of$latex G (x)$$; i.e., 230 latex $231 g[0] = G (x) 232$$$
233
234 $head gx$$235 The argument icode gx$$ has prototype 236$codei%
237  %BAvector% &%gx%
238 %$$239 and has size latex n$$.
240 The input values of its elements do not matter.
241 On output,
242 it contains the Jacobian of $latex G (x)$$; i.e., 243 for latex j = 0 , \ldots , n-1$$, 244$latex $245 gx[ j ] = G^{(1)} (x)_j 246$ $$247 248 head gxx$$
249 The argument $icode gx$$has prototype 250 codei% 251 %BAvector% &%gxx% 252 %$$ 253 and has size$latex n \times n$$. 254 The input values of its elements do not matter. 255 On output, 256 it contains the Hessian of latex G (x)$$; i.e.,
257 for $latex i = 0 , \ldots , n-1$$, and 258 latex j = 0 , \ldots , n-1$$, 259$latex $260 gxx[ i * n + j ] = G^{(2)} (x)_{i,j} 261$ $$262 263 head BAvector$$
264 The type $icode BAvector$$must be a 265 cref SimpleVector$$ class. 266 We use$icode Base$$to refer to the type of the elements of 267 icode BAvector$$; i.e.,
268 $codei% 269 %BAvector%::value_type 270 %$$271 272 head ADvector$$ 273 The type$icode ADvector$$must be a 274 cref SimpleVector$$ class with elements of type
275 $codei%AD<%Base%>%$$; i.e., 276 codei% 277 %ADvector%::value_type 278 %$$ 279 must be the same type as 280$codei%
282 %$$. 283 284 285 head Example$$
286 $children% 287 example/general/bender_quad.cpp 288 %$$289 The file 290 cref bender_quad.cpp$$ 291 contains an example and test of this operation. 292 It returns true if it succeeds and false otherwise. 293 294 295$end
296 -----------------------------------------------------------------------------
297 */
298
300
301 template <class BAvector, class Fun>
303  const BAvector &x ,
304  const BAvector &y ,
305  Fun fun ,
306  BAvector &g ,
307  BAvector &gx ,
308  BAvector &gxx )
309 { // determine the base type
310  typedef typename BAvector::value_type Base;
311
312  // check that BAvector is a SimpleVector class
313  CheckSimpleVector<Base, BAvector>();
314
315  // declare the ADvector type
317
318  // size of the x and y spaces
319  size_t n = size_t(x.size());
320  size_t m = size_t(y.size());
321
322  // check the size of gx and gxx
324  g.size() == 1,
325  "BenderQuad: size of the vector g is not equal to 1"
326  );
328  size_t(gx.size()) == n,
329  "BenderQuad: size of the vector gx is not equal to n"
330  );
332  size_t(gxx.size()) == n * n,
333  "BenderQuad: size of the vector gxx is not equal to n * n"
334  );
335
336  // some temporary indices
337  size_t i, j;
338
339  // variable versions x
341  for(j = 0; j < n; j++)
342  vx[j] = x[j];
343
344  // declare the independent variables
345  Independent(vx);
346
347  // evaluate h = H(x, y)
349  h = fun.h(vx, y);
350
351  // evaluate dy (x) = Newton step as a function of x through h only
353  dy = fun.dy(x, y, h);
354
355  // variable version of y
357  for(j = 0; j < m; j++)
358  vy[j] = y[j] + dy[j];
359
360  // evaluate G~ (x) = F [ x , y + dy(x) ]
362  gtilde = fun.f(vx, vy);
363
364  // AD function object that corresponds to G~ (x)
365  // We will make heavy use of this tape, so optimize it
367  Gtilde.Dependent(vx, gtilde);
368  Gtilde.optimize();
369
370  // value of G(x)
371  g = Gtilde.Forward(0, x);
372
373  // initial forward direction vector as zero
374  BAvector dx(n);
375  for(j = 0; j < n; j++)
376  dx[j] = Base(0.0);
377
378  // weight, first and second order derivative values
379  BAvector dg(1), w(1), ddw(2 * n);
380  w[0] = 1.;
381
382
383  // Jacobian and Hessian of G(x) is equal Jacobian and Hessian of Gtilde
384  for(j = 0; j < n; j++)
385  { // compute partials in x[j] direction
386  dx[j] = Base(1.0);
387  dg = Gtilde.Forward(1, dx);
388  gx[j] = dg[0];
389
390  // restore the dx vector to zero
391  dx[j] = Base(0.0);
392
393  // compute second partials w.r.t x[j] and x[l] for l = 1, n
394  ddw = Gtilde.Reverse(2, w);
395  for(i = 0; i < n; i++)
396  gxx[ i * n + j ] = ddw[ i * 2 + 1 ];
397  }
398
399  return;
400 }
401
402 } // END CppAD namespace
403
404 # endif
Check that exp is true, if not print msg and terminate execution.
Class used to hold function objects.
Declaration of independent variables.
void BenderQuad(const BAvector &x, const BAvector &y, Fun fun, BAvector &g, BAvector &gx, BAvector &gxx)