00001 # ifndef CPPAD_REV_SPARSE_HES_INCLUDED 00002 # define CPPAD_REV_SPARSE_HES_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 /* 00016 $begin RevSparseHes$$ 00017 $spell 00018 VecAD 00019 Jacobian 00020 Jac 00021 Hessian 00022 Hes 00023 const 00024 Bool 00025 Dep 00026 proportional 00027 var 00028 $$ 00029 00030 $section Hessian Sparsity Pattern: Reverse Mode$$ 00031 00032 $index RevSparseHes$$ 00033 $index reverse, sparse Hessian$$ 00034 $index sparse, reverse Hessian$$ 00035 $index pattern, reverse Hessian$$ 00036 00037 $head Syntax$$ 00038 $syntax%%h% = %f%.RevSparseHes(%q%, %s%)%$$ 00039 $pre 00040 $$ 00041 00042 00043 $head Purpose$$ 00044 We use $latex F : B^n \rightarrow B^m$$ to denote the 00045 $xref/glossary/AD Function/AD function/$$ corresponding to $italic f$$. 00046 For a fixed $latex n \times q$$ matrix $latex R$$ 00047 and a fixed $latex m \times 1$$ matrix $latex S$$, 00048 the second partial of $latex S * F[ x + R * u ]$$ 00049 with respect to $latex u$$ at $latex u = 0$$ and with respect to x 00050 $latex \[ 00051 H(x) = R^T * (S * F)^{(2)} ( x ) 00052 \] $$ 00053 where $latex (S * F)^{(2)} (x)$$ is the Hessian of the scalar 00054 valued function $latex S * F (x)$$. 00055 Given a 00056 $xref/glossary/Sparsity Pattern/sparsity pattern/$$ 00057 for $latex R$$ and $latex S$$, 00058 $code RevSparseHes$$ returns a sparsity pattern for the $latex H(x)$$. 00059 00060 $head f$$ 00061 The object $italic f$$ has prototype 00062 $syntax% 00063 const ADFun<%Base%> %f% 00064 %$$ 00065 00066 $head x$$ 00067 If no $xref/VecAD/$$ objects are used by the 00068 AD of $italic Base$$ 00069 $xref/glossary/Operation/Sequence/operation sequence/1/$$ 00070 stored in $italic f$$, 00071 the sparsity pattern is valid for all values $latex x \in B^n$$. 00072 $pre 00073 00074 $$ 00075 If $xref/SeqProperty/use_VecAD/f.use_VecAD/$$ is true, 00076 the sparsity patter is only valid for the value of $italic x$$ 00077 in the previous $xref/ForwardZero//zero order forward mode/$$ call 00078 $syntax% 00079 %f%.Forward(0, %x%) 00080 %$$ 00081 If there is no previous zero order forward mode call using $italic f$$, 00082 the value of the $xref/Independent//independent/$$ variables 00083 during the recording of the AD sequence of operations is used 00084 for $italic x$$. 00085 00086 $head q$$ 00087 The argument $italic q$$ has prototype 00088 $syntax% 00089 size_t %q% 00090 %$$ 00091 It specifies the number of columns in the Jacobian $latex J(x)$$. 00092 It must be the same value as in the previous $xref/ForSparseJac/$$ call 00093 $syntax% 00094 %f%.ForSparseJac(%q%, %r%) 00095 %$$ 00096 Note that the memory required for the calculation is proportional 00097 to $latex q$$ times the total number of variables 00098 in the AD operation sequence corresponding to $italic f$$ 00099 ($xref/SeqProperty/size_var/f.size_var/$$). 00100 00101 $head r$$ 00102 The argument $italic r$$ in the previous call 00103 $syntax% 00104 %f%.ForSparseJac(%q%, %r%) 00105 %$$ 00106 is a sparsity pattern for the matrix $latex R$$ above; i.e., 00107 for $latex i = 0 , \ldots , n-1$$ and $latex j = 0 , \ldots , q-1$$. 00108 $latex \[ 00109 R_{i,j} \neq 0 ; \Rightarrow \; r [ i * q + j ] = {\rm true} 00110 \] $$ 00111 00112 $head s$$ 00113 The argument $italic s$$ has prototype 00114 $syntax% 00115 const %Vector% &%s% 00116 %$$ 00117 (see $xref/RevSparseHes/Vector/Vector/$$ below) 00118 and its size is $latex m$$. 00119 It specifies a 00120 $xref/glossary/Sparsity Pattern/sparsity pattern/$$ 00121 for the matrix $italic S$$ as follows: 00122 for $latex j = 0 , \ldots , m-1$$. 00123 $latex \[ 00124 S_{0,j} \neq 0 ; \Rightarrow \; s [ j ] = {\rm true} 00125 \] $$ 00126 00127 $head t$$ 00128 The result $italic h$$ has prototype 00129 $syntax% 00130 %Vector% &%h% 00131 %$$ 00132 (see $xref/RevSparseHes/Vector/Vector/$$ below) 00133 and its size is $latex q * n$$, 00134 It specifies a 00135 $xref/glossary/Sparsity Pattern/sparsity pattern/$$ 00136 for the matrix $latex H(x)$$ as follows: 00137 for $latex x \in B^n$$, 00138 for $latex i = 0 , \ldots , q-1$$, 00139 and $latex j = 0 , \ldots , n-1$$ 00140 $latex \[ 00141 H(x)_{i,j} \neq 0 ; \Rightarrow \; h [ i * n + j ] = {\rm true} 00142 \] $$ 00143 00144 $head Vector$$ 00145 The type $italic Vector$$ must be a $xref/SimpleVector/$$ class with 00146 $xref/SimpleVector/Elements of Specified Type/elements of type bool/$$. 00147 The routine $xref/CheckSimpleVector/$$ will generate an error message 00148 if this is not the case. 00149 In order to save memory, 00150 you may want to use a class that packs multiple elements into one 00151 storage location; for example, 00152 $xref/CppAD_vector/vectorBool/vectorBool/$$. 00153 00154 $head Entire Sparsity Pattern$$ 00155 Suppose that $latex q = n$$ and 00156 $latex R$$ is the $latex n \times n$$ identity matrix, 00157 If follows that 00158 $latex \[ 00159 r [ i * q + j ] = \left\{ \begin{array}{ll} 00160 {\rm true} & {\rm if} \; i = j \\ 00161 {\rm false} & {\rm otherwise} 00162 \end{array} \right. 00163 \] $$ 00164 is an efficient sparsity pattern for $latex R$$; 00165 i.e., the choice for $italic r$$ has as few true values as possible. 00166 Further suppose that the $latex S$$ is the $th k$$ 00167 $xref/glossary/Elementary Vector/elementary vector/$$ 00168 If follows that 00169 $latex \[ 00170 s [ j ] = \left\{ \begin{array}{ll} 00171 {\rm true} & {\rm if} \; j = k \\ 00172 {\rm false} & {\rm otherwise} 00173 \end{array} \right. 00174 \] $$ 00175 is an efficient sparsity pattern for $latex S$$. 00176 In the case defined above, 00177 the result $italic h$$ corresponds to a 00178 sparsity pattern for the Hessian $latex F_k^{(2)} (x)$$. 00179 00180 $head Example$$ 00181 $children% 00182 example/rev_sparse_hes.cpp 00183 %$$ 00184 The file 00185 $xref/RevSparseHes.cpp/$$ 00186 contains an example and test of this operation. 00187 It returns true if it succeeds and false otherwise. 00188 00189 $end 00190 ----------------------------------------------------------------------------- 00191 */ 00192 # include <algorithm> 00193 00194 // BEGIN CppAD namespace 00195 namespace CppAD { 00196 00197 template <class Base> 00198 template <class Vector> 00199 Vector ADFun<Base>::RevSparseHes(size_t q, const Vector &s) const 00200 { 00201 // type used to pack bits (must support standard bit operations) 00202 typedef size_t Pack; 00203 00204 // temporary indices 00205 size_t i, j, k, p; 00206 00207 // check Vector is Simple Vector class with bool elements 00208 CheckSimpleVector<bool, Vector>(); 00209 00210 // range and domain dimensions for F 00211 size_t m = dep_taddr.size(); 00212 size_t n = ind_taddr.size(); 00213 00214 CPPAD_ASSERT_KNOWN( 00215 q > 0, 00216 "RevSparseHes: q (first argument) is not greater than zero" 00217 ); 00218 CPPAD_ASSERT_KNOWN( 00219 q == ForJacBitDim, 00220 "RevSparseHes: q (first argument) is not equal to its value" 00221 " in the previous call to ForSparseJac with this ADFun object." 00222 ); 00223 CPPAD_ASSERT_KNOWN( 00224 s.size() == m, 00225 "RevSparseHes: s (third argument) length is not equal to\n" 00226 "range dimension for ADFun object." 00227 ); 00228 00229 // number of packed values per variable on the tape 00230 size_t npv = 1 + (q - 1) / sizeof(Pack); 00231 CPPAD_ASSERT_UNKNOWN( npv <= ForJacColDim ); 00232 00233 // array that will hold packed reverse Jacobian values 00234 Pack *RevJac = CPPAD_NULL; 00235 RevJac = CPPAD_TRACK_NEW_VEC(totalNumVar, RevJac); 00236 00237 // array that will hold packed reverse Hessain values 00238 Pack *RevHes = CPPAD_NULL; 00239 RevHes = CPPAD_TRACK_NEW_VEC(totalNumVar * npv, RevHes); 00240 00241 // update maximum memory requirement 00242 // memoryMax = std::max( memoryMax, 00243 // Memory() + totalNumVar * (npv + 1) * sizeof(Pack) 00244 // ); 00245 00246 // initialize entire RevHes and RevJac matrix to false 00247 for(i = 0; i < totalNumVar; i++) 00248 { RevJac[i] = 0; 00249 for(k = 0; k < npv; k++) 00250 RevHes[ i * npv + k ] = 0; 00251 } 00252 for(i = 0; i < m; i++) 00253 { CPPAD_ASSERT_UNKNOWN( dep_taddr[i] < totalNumVar ); 00254 if( s[i] ) 00255 { // set all the bits to true 00256 RevJac[ dep_taddr[i] ] = ~0; 00257 } 00258 } 00259 00260 // comput the reverse mode Jacobian sparsity 00261 RevJacSweep(1, totalNumVar, &Rec, TaylorColDim, Taylor, RevJac); 00262 00263 // compute the Hessian sparsity patterns 00264 RevHesSweep(npv, totalNumVar, &Rec, TaylorColDim, Taylor, ForJac, 00265 RevJac, RevHes 00266 ); 00267 00268 // return values corresponding to independent variables 00269 Pack mask; 00270 Vector h(n * q); 00271 00272 // j is index corresponding to reverse mode martial 00273 for(j = 0; j < n; j++) 00274 { CPPAD_ASSERT_UNKNOWN( ind_taddr[j] < totalNumVar ); 00275 00276 // ind_taddr[j] is operator taddr for j-th independent variable 00277 CPPAD_ASSERT_UNKNOWN( Rec.GetOp( ind_taddr[j] ) == InvOp ); 00278 00279 // i is index corresponding to forward mode partial 00280 for(i = 0; i < q; i++) 00281 { k = i / sizeof(Pack); 00282 p = i - k * sizeof(Pack); 00283 mask = Pack(1) << p; 00284 mask &= RevHes[ ind_taddr[j] * npv + k ]; 00285 h[ i * n + j ] = (mask != 0); 00286 } 00287 } 00288 00289 // free memory used for the calculation 00290 CPPAD_TRACK_DEL_VEC(RevJac); 00291 CPPAD_TRACK_DEL_VEC(RevHes); 00292 00293 return h; 00294 } 00295 00296 } // END CppAD namespace 00297 00298 # endif