CppAD: A C++ Algorithmic Differentiation Package 20110419
sparse_hessian.hpp
Go to the documentation of this file.
00001 /* $Id$ */
00002 # ifndef CPPAD_SPARSE_HESSIAN_INCLUDED
00003 # define CPPAD_SPARSE_HESSIAN_INCLUDED
00004 
00005 /* --------------------------------------------------------------------------
00006 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-10 Bradley M. Bell
00007 
00008 CppAD is distributed under multiple licenses. This distribution is under
00009 the terms of the 
00010                     Common 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 /*
00017 $begin sparse_hessian$$
00018 $spell
00019         CppAD
00020         valarray
00021         std
00022         Bool
00023         hes
00024         const
00025         Taylor
00026 $$
00027 
00028 $section Sparse Hessian: Easy Driver$$
00029 $index SparseHessian$$
00030 $index hessian, sparse$$
00031 
00032 $head Syntax$$
00033 $codei%%hes% = %f%.SparseHessian(%x%, %w%)
00034 %$$
00035 $codei%%hes% = %f%.SparseHessian(%x%, %w%, %p%)%$$
00036 
00037 $head Purpose$$
00038 We use $latex F : \R^n \rightarrow \R^m$$ do denote the
00039 $cref/AD function/glossary/AD Function/$$
00040 corresponding to $icode f$$. 
00041 The syntax above sets $icode hes$$ to the Hessian 
00042 $latex \[
00043         hes = \dpow{2}{x} \sum_{i=1}^m w_i F_i (x) 
00044 \] $$
00045 This routine assumes  that the matrix $icode hes$$ is sparse
00046 and uses this assumption sparse to reduce the amount of computation necessary.
00047 One should use speed tests to verify that results are computed faster
00048 than when using the routine $cref/Hessian/$$.
00049 
00050 $head f$$
00051 The object $icode f$$ has prototype
00052 $codei%
00053         ADFun<%Base%> %f%
00054 %$$
00055 Note that the $cref/ADFun/$$ object $icode f$$ is not $code const$$
00056 (see $cref/Uses Forward/sparse_hessian/Uses Forward/$$ below).
00057 
00058 $head x$$
00059 The argument $icode x$$ has prototype
00060 $codei%
00061         const %VectorBase%& %x%
00062 %$$
00063 (see $cref/VectorBase/sparse_hessian/VectorBase/$$ below)
00064 and its size 
00065 must be equal to $icode n$$, the dimension of the
00066 $cref/domain/seq_property/Domain/$$ space for $icode f$$.
00067 It specifies
00068 that point at which to evaluate the Hessian.
00069 
00070 $head w$$
00071 The argument $icode w$$ has prototype
00072 $codei%
00073         const %VectorBase%& %w%
00074 %$$
00075 and size $latex m$$.
00076 It specifies the value of $latex w_i$$ in the expression 
00077 for $icode hes$$.
00078 The more components of $latex w$$ that are identically zero,
00079 the more sparse the resulting Hessian may be (and hence the more efficient
00080 the calculation of $icode hes$$ may be).
00081 
00082 $head p$$
00083 The argument $icode p$$ is optional and has prototype
00084 $codei%
00085         const %VectorSet%& %p%
00086 %$$
00087 (see $cref/VectorSet/sparse_hessian/VectorSet/$$ below)
00088 If it has elements of type $code bool$$,
00089 its size is $latex n * n$$.
00090 If it has elements of type $code std::set<size_t>$$,
00091 its size is $latex n$$ and all its set elements are between
00092 zero and $latex n - 1$$.
00093 It specifies a
00094 $cref/sparsity pattern/glossary/Sparsity Pattern/$$ 
00095 for the Hessian $icode hes$$.
00096 $pre
00097 
00098 $$
00099 If this sparsity pattern does not change between calls to 
00100 $codei SparseHessian$$, it should be faster to calculate $icode p$$ once and
00101 pass this argument to $codei SparseHessian$$.
00102 In addition,
00103 if you specify $icode p$$, CppAD will use the same
00104 type of sparsity representation
00105 (vectors of $code bool$$ or vectors of $code std::set<size_t>$$)
00106 for its internal calculations.
00107 Otherwise, the representation
00108 for the internal calculations is unspecified.
00109 
00110 $head hes$$
00111 The result $icode hes$$ has prototype
00112 $codei%
00113         %VectorBase% %hes%
00114 %$$
00115 and its size is $latex n * n$$.
00116 For $latex j = 0 , \ldots , n - 1 $$ 
00117 and $latex \ell = 0 , \ldots , n - 1$$
00118 $latex \[
00119         hes [ j * n + \ell ] = \DD{ w^{\rm T} F }{ x_j }{ x_\ell } ( x )
00120 \] $$
00121 
00122 $head VectorBase$$
00123 The type $icode VectorBase$$ must be a $cref/SimpleVector/$$ class with
00124 $cref/elements of type/SimpleVector/Elements of Specified Type/$$
00125 $icode Base$$.
00126 The routine $cref/CheckSimpleVector/$$ will generate an error message
00127 if this is not the case.
00128 
00129 $head VectorSet$$
00130 The type $icode VectorSet$$ must be a $xref/SimpleVector/$$ class with
00131 $xref/SimpleVector/Elements of Specified Type/elements of type/$$
00132 $code bool$$ or $code std::set<size_t>$$;
00133 see $cref/sparsity pattern/glossary/Sparsity Pattern/$$ for a discussion
00134 of the difference.
00135 The routine $xref/CheckSimpleVector/$$ will generate an error message
00136 if this is not the case.
00137 
00138 $subhead Restrictions$$
00139 If $icode VectorSet$$ has elements of $code std::set<size_t>$$,
00140 then $icode%p%[%i%]%$$ must return a reference (not a copy) to the 
00141 corresponding set.
00142 According to section 26.3.2.3 of the 1998 C++ standard,
00143 $code std::valarray< std::set<size_t> >$$ does not satisfy
00144 this condition. 
00145 
00146 $head Uses Forward$$
00147 After each call to $cref/Forward/$$,
00148 the object $icode f$$ contains the corresponding 
00149 $cref/Taylor coefficients/glossary/Taylor Coefficient/$$.
00150 After $code SparseHessian$$,
00151 the previous calls to $xref/Forward/$$ are undefined.
00152 
00153 $head Example$$
00154 $children%
00155         example/sparse_hessian.cpp
00156 %$$
00157 The routine
00158 $cref/sparse_hessian.cpp/$$ 
00159 is examples and tests of $code sparse_hessian$$.
00160 It return $code true$$, if it succeeds and $code false$$ otherwise.
00161 
00162 
00163 $end
00164 -----------------------------------------------------------------------------
00165 */
00166 
00167 CPPAD_BEGIN_NAMESPACE
00168 /*!
00169 \file sparse_hessian.hpp
00170 Sparse Hessian driver routine and helper functions.
00171 */
00172 
00173 /*!
00174 Private helper function for SparseHessian(x, w, p).
00175 
00176 All of the description in the public member function SparseHessian(x, w, p)
00177 applies.
00178 
00179 \param set_type
00180 is a \c bool value. This argument is used to dispatch to the proper source
00181 code depending on the value of \c VectorSet::value_type.
00182 
00183 \param x
00184 See \c SparseHessian(x, w, p).
00185 
00186 \param w
00187 See \c SparseHessian(x, w, p).
00188 
00189 \param p
00190 See \c SparseHessian(x, w, p).
00191 
00192 \param hes
00193 is the return value for the corresponding call to \c SparseHessian(x, w, p).
00194 On input, it must have size equal to the domain times range dimension
00195 for this ADFun<Base> object.
00196 On return, it will contain the Hessian.
00197 */
00198 
00199 template <class Base>
00200 template <class VectorBase, class VectorSet>
00201 void ADFun<Base>::SparseHessianCase(
00202         bool               set_type         ,
00203         const VectorBase&  x                ,
00204         const VectorBase&  w                ,
00205         const VectorSet&   p                ,
00206         VectorBase&        hes              )
00207 {
00208         typedef CppAD::vector<size_t> SizeVector;
00209         typedef CppAD::vectorBool     VectorBool;
00210         size_t j, k, l;
00211 
00212         size_t n = Domain();
00213 
00214         // check Vector is Simple Vector class with bool elements
00215         CheckSimpleVector<bool, VectorSet>();
00216 
00217         // check Vector is Simple Vector class with Base type elements
00218         CheckSimpleVector<Base, VectorBase>();
00219 
00220         CPPAD_ASSERT_KNOWN(
00221                 x.size() == n,
00222                 "SparseHessian: size of x not equal domain dimension for f"
00223         ); 
00224         CPPAD_ASSERT_KNOWN(
00225                 w.size() == Range(),
00226                 "SparseHessian: size of w not equal range dimension for f"
00227         );
00228         CPPAD_ASSERT_KNOWN(
00229                 p.size() == n * n,
00230                 "SparseHessian: using bool values and size of p "
00231                 "not equal square of domain dimension for f"
00232         );
00233         
00234         // initial coloring
00235         SizeVector color(n);
00236         for(j = 0; j < n; j++)
00237                 color[j] = j;
00238 
00239         // See GreedyPartialD2Coloring Algorithm Section 3.6.2 of
00240         // Graph Coloring in Optimization Revisited by
00241         // Assefaw Gebremedhin, Fredrik Maane, Alex Pothen
00242         VectorBool forbidden(n);
00243         for(j = 0; j < n; j++)
00244         {       // initial all colors as ok for this column
00245                 for(k = 0; k < n; k++)
00246                         forbidden[k] = false;
00247                 // for each row that is connected to column j
00248                 for(l = 0; l < n; l++) if( p[l * n + j] )
00249                 {       // for each column that is connected to row l
00250                         for(k = 0; k < n; k++) if( p[l * n + k] & (j != k) )    
00251                                 forbidden[ color[k] ] = true;
00252                 }
00253                 k = 0;
00254                 while( forbidden[k] && k < n )
00255                 {       k++;
00256                         CPPAD_ASSERT_UNKNOWN( k < n );
00257                 }
00258                 color[j] = k;
00259         }
00260         size_t n_color = 1;
00261         for(k = 0; k < n; k++) n_color = std::max(n_color, color[k] + 1);
00262 
00263         // some values
00264         const Base zero(0);
00265         const Base one(1);
00266 
00267         // point at which we are evaluating the Hessian
00268         Forward(0, x);
00269 
00270         // initialize the return value
00271         for(j = 0; j < n; j++)
00272                 for(k = 0; k < n; k++)
00273                         hes[j * n + k] = zero;
00274 
00275         // direction vector for calls to forward
00276         VectorBase u(n);
00277 
00278         // location for return values from Reverse
00279         VectorBase ddw(n * 2);
00280 
00281         // loop over colors
00282         size_t c;
00283         for(c = 0; c < n_color; c++)
00284         {       // determine all the colums with this color
00285                 for(j = 0; j < n; j++)
00286                 {       if( color[j] == c )
00287                                 u[j] = one;
00288                         else    u[j] = zero;
00289                 }
00290                 // call forward mode for all these columns at once
00291                 Forward(1, u);
00292 
00293                 // evaluate derivative of w^T * F'(x) * u
00294                 ddw = Reverse(2, w);
00295 
00296                 // set the corresponding components of the result
00297                 for(j = 0; j < n; j++) if( color[j] == c )
00298                 {       for(l = 0; l < n; l++) 
00299                                 if( p[ l * n + j ] )
00300                                         hes[l * n + j] = ddw[l * 2 + 1];
00301                 }
00302         }
00303 }
00304 /*!
00305 Private helper function for SparseHessian(x, w, p).
00306 
00307 All of the description in the public member function SparseHessian(x, w, p)
00308 applies.
00309 
00310 \param set_type
00311 is a \c std::set<size_t> value.
00312 This argument is used to dispatch to the proper source
00313 code depending on the value of \c VectorSet::value_type.
00314 
00315 \param x
00316 See \c SparseHessian(x, w, p).
00317 
00318 \param w
00319 See \c SparseHessian(x, w, p).
00320 
00321 \param p
00322 See \c SparseHessian(x, w, p).
00323 
00324 \param hes
00325 is the return value for the corresponding call to \c SparseHessian(x, w, p).
00326 On input, it must have size equal to the domain times range dimension
00327 for this ADFun<Base> object.
00328 On return, it will contain the Hessian.
00329 */
00330 
00331 template <class Base>
00332 template <class VectorBase, class VectorSet>
00333 void ADFun<Base>::SparseHessianCase(
00334         const std::set<size_t>&  set_type         ,
00335         const VectorBase&        x                ,
00336         const VectorBase&        w                ,
00337         const VectorSet&         p                ,
00338         VectorBase&              hes              )
00339 {
00340         typedef CppAD::vector<size_t> SizeVector;
00341         typedef CppAD::vectorBool     VectorBool;
00342         size_t j, k, l;
00343 
00344         size_t n = Domain();
00345 
00346         // check VectorSet is Simple Vector class with sets for elements
00347         static std::set<size_t> two, three;
00348         if( two.empty() )
00349         {       two.insert(2);
00350                 three.insert(3);
00351         }
00352         CPPAD_ASSERT_UNKNOWN( two.size() == 1 );
00353         CPPAD_ASSERT_UNKNOWN( three.size() == 1 );
00354         CheckSimpleVector<std::set<size_t>, VectorSet>(two, three);
00355 
00356         // check Vector is Simple Vector class with Base type elements
00357         CheckSimpleVector<Base, VectorBase>();
00358 
00359         CPPAD_ASSERT_KNOWN(
00360                 x.size() == n,
00361                 "SparseHessian: size of x not equal domain dimension for f"
00362         ); 
00363         CPPAD_ASSERT_KNOWN(
00364                 w.size() == Range(),
00365                 "SparseHessian: size of w not equal range dimension for f"
00366         );
00367         CPPAD_ASSERT_KNOWN(
00368                 p.size() == n,
00369                 "SparseHessian: using set values and size of p "
00370                 "not equal domain dimension for f"
00371         );
00372         
00373         // initial coloring
00374         SizeVector color(n);
00375         for(j = 0; j < n; j++)
00376                 color[j] = j;
00377 
00378         // See GreedyPartialD2Coloring Algorithm Section 3.6.2 of
00379         // Graph Coloring in Optimization Revisited by
00380         // Assefaw Gebremedhin, Fredrik Maane, Alex Pothen
00381         VectorBool forbidden(n);
00382         std::set<size_t>::const_iterator itr_k, itr_l;
00383         for(j = 0; j < n; j++)
00384         {       // initial all colors as ok for this column
00385                 for(k = 0; k < n; k++)
00386                         forbidden[k] = false;
00387 
00388                 // for each row that is connected to column j
00389                 itr_l = p[j].begin();
00390                 while( itr_l != p[j].end() )
00391                 {       l = *itr_l++;
00392                         // for each column that is connected to row l
00393                         // (same as each row that is connect to column l)
00394                         itr_k = p[l].begin();
00395                         while( itr_k != p[l].end() )
00396                         {       k = *itr_k++;
00397                                 forbidden[ color[k] ] = (j != k);
00398                         }
00399                 }
00400                 k = 0;
00401                 while( forbidden[k] && k < n )
00402                 {       k++;
00403                         CPPAD_ASSERT_UNKNOWN( k < n );
00404                 }
00405                 color[j] = k;
00406         }
00407         size_t n_color = 1;
00408         for(k = 0; k < n; k++) n_color = std::max(n_color, color[k] + 1);
00409 
00410         // some values
00411         const Base zero(0);
00412         const Base one(1);
00413 
00414         // point at which we are evaluating the Hessian
00415         Forward(0, x);
00416 
00417         // initialize the return value
00418         for(j = 0; j < n; j++)
00419                 for(k = 0; k < n; k++)
00420                         hes[j * n + k] = zero;
00421 
00422         // direction vector for calls to forward
00423         VectorBase u(n);
00424 
00425         // location for return values from Reverse
00426         VectorBase ddw(n * 2);
00427 
00428         // loop over colors
00429         size_t c;
00430         for(c = 0; c < n_color; c++)
00431         {       // determine all the colums with this color
00432                 for(j = 0; j < n; j++)
00433                 {       if( color[j] == c )
00434                                 u[j] = one;
00435                         else    u[j] = zero;
00436                 }
00437                 // call forward mode for all these columns at once
00438                 Forward(1, u);
00439 
00440                 // evaluate derivative of w^T * F'(x) * u
00441                 ddw = Reverse(2, w);
00442 
00443                 // set the corresponding components of the result
00444                 for(j = 0; j < n; j++) if( color[j] == c )
00445                 {       itr_l = p[j].begin();
00446                         while( itr_l != p[j].end() )
00447                         {       l = *itr_l++;
00448                                 hes[l * n + j] = ddw[l * 2 + 1];
00449                         }
00450                 }
00451         }
00452 }
00453 
00454 /*!
00455 Compute a sparse Hessian.
00456 
00457 The C++ source code coresponding to this operation is
00458 \verbatim
00459         hes = SparseHessian(x, w)
00460 \endverbatim
00461 
00462 
00463 \tparam Base
00464 is the base type for the recording that is stored in this
00465 ADFun<Base object.
00466 
00467 \tparam VectorBase
00468 is a simple vector class with elements of the \a Base.
00469 
00470 \tparam VectorSet
00471 is a simple vector class with elements of type
00472 \c bool or \c std::set<size_t>.
00473 
00474 \param x
00475 is a vector specifing the point at which to compute the Hessian.
00476 
00477 \param w
00478 The Hessian is computed for a weighted sum of the components
00479 of the function corresponding to this ADFun<Base> object.
00480 The argument \a w specifies the weights for each component.
00481 It must have size equal to the range dimension for this ADFun<Base> object.
00482 
00483 \param p
00484 is a sparsity pattern for the Hessian.
00485 
00486 \return
00487 Will be a vector of size \c n * n containing the Hessian of 
00488 at the point specified by \a x
00489 (where \c n is the domain dimension for this ADFun<Base> object).
00490 */
00491 template <class Base>
00492 template <class VectorBase, class VectorSet>
00493 VectorBase ADFun<Base>::SparseHessian(
00494         const VectorBase& x, const VectorBase& w, const VectorSet& p
00495 )
00496 {       size_t n = Domain();
00497         VectorBase hes(n * n);
00498 
00499         typedef typename VectorSet::value_type Set_type;
00500         SparseHessianCase(Set_type(), x, w, p, hes);
00501 
00502         return hes;
00503 }
00504 
00505 /*!
00506 Compute a sparse Hessian
00507 
00508 The C++ source code coresponding to this operation is
00509 \verbatim
00510         hes = SparseHessian(x, w)
00511 \endverbatim
00512 
00513 
00514 \tparam Base
00515 is the base type for the recording that is stored in this
00516 ADFun<Base object.
00517 
00518 \tparam VectorBase
00519 is a simple vector class with elements of the \a Base.
00520 
00521 \param x
00522 is a vector specifing the point at which to compute the Hessian.
00523 
00524 \param w
00525 The Hessian is computed for a weighted sum of the components
00526 of the function corresponding to this ADFun<Base> object.
00527 The argument \a w specifies the weights for each component.
00528 It must have size equal to the range dimension for this ADFun<Base> object.
00529 
00530 \return
00531 Will be a vector of size \c n * n containing the Hessian of 
00532 at the point specified by \a x
00533 (where \c n is the domain dimension for this ADFun<Base> object).
00534 */
00535 template <class Base>
00536 template <class VectorBase>
00537 VectorBase ADFun<Base>::SparseHessian(const VectorBase &x, const VectorBase &w)
00538 {       size_t i, j, k;
00539         typedef CppAD::vectorBool VectorBool;
00540 
00541         size_t m = Range();
00542         size_t n = Domain();
00543 
00544         // determine the sparsity pattern p for Hessian of w^T F
00545         VectorBool r(n * n);
00546         for(j = 0; j < n; j++)
00547         {       for(k = 0; k < n; k++)
00548                         r[j * n + k] = false;
00549                 r[j * n + j] = true;
00550         }
00551         ForSparseJac(n, r);
00552         //
00553         VectorBool s(m);
00554         for(i = 0; i < m; i++)
00555                 s[i] = w[i] != 0;
00556         VectorBool p = RevSparseHes(n, s);
00557 
00558         // compute sparse Hessian
00559         VectorBase hes(n * n);
00560         bool set_type = true; // only used to dispatch compiler to proper case
00561         SparseHessianCase(set_type, x, w, p, hes);
00562 
00563         return hes;
00564 }
00565 
00566 CPPAD_END_NAMESPACE
00567 # endif