CppAD: A C++ Algorithmic Differentiation Package  20130918
bender_quad.hpp
Go to the documentation of this file.
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
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines