CppAD: A C++ Algorithmic Differentiation Package
20130102
|
00001 /* $Id$ */ 00002 /* -------------------------------------------------------------------------- 00003 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-12 Bradley M. Bell 00004 00005 CppAD is distributed under multiple licenses. This distribution is under 00006 the terms of the 00007 Eclipse Public License Version 1.0. 00008 00009 A copy of this license is included in the COPYING file of this distribution. 00010 Please visit http://www.coin-or.org/CppAD/ for information on other licenses. 00011 -------------------------------------------------------------------------- */ 00012 00013 /* 00014 $begin link_sparse_hessian$$ 00015 $spell 00016 const 00017 bool 00018 CppAD 00019 $$ 00020 00021 $index link_sparse_hessian$$ 00022 $index sparse, speed test$$ 00023 $index speed, test sparse$$ 00024 $index test, sparse speed$$ 00025 00026 $section Speed Testing Sparse Hessian$$ 00027 00028 $head Prototype$$ 00029 $codei%extern bool link_sparse_hessian( 00030 size_t %size% , 00031 size_t %repeat% , 00032 CppAD::vector<double>& %x% , 00033 const CppAD::vector<size_t>& %row% , 00034 const CppAD::vector<size_t>& %col% , 00035 CppAD::vector<double>& %hessian% 00036 ); 00037 %$$ 00038 00039 $head Method$$ 00040 Given a row index vector $latex row$$ 00041 and a second column vector $latex col$$, 00042 the corresponding function 00043 $latex f : \B{R}^n \rightarrow \B{R} $$ 00044 is defined by $cref sparse_hes_fun$$. 00045 The non-zero entries in the Hessian of this function have 00046 one of the following forms: 00047 $latex \[ 00048 \DD{f}{x[row[k]]}{x[row[k]]} 00049 \; , \; 00050 \DD{f}{x[row[k]]}{x[col[k]]} 00051 \; , \; 00052 \DD{f}{x[col[k]]}{x[row[k]]} 00053 \; , \; 00054 \DD{f}{x[col[k]]}{x[col[k]]} 00055 \] $$ 00056 for some $latex k $$ between zero and $latex K-1 $$. 00057 All the other terms of the Hessian are zero. 00058 00059 $head size$$ 00060 The argument $icode size$$, referred to as $latex n$$ below, 00061 is the dimension of the domain space for $latex f(x)$$. 00062 00063 $head repeat$$ 00064 The argument $icode repeat$$ is the number of times 00065 to repeat the test 00066 (with a different value for $icode x$$ corresponding to 00067 each repetition). 00068 00069 $head x$$ 00070 The argument $icode x$$ has prototype 00071 $codei% 00072 CppAD::vector<double> &%x% 00073 %$$ 00074 and its size is $latex n$$; i.e., $icode%x%.size() == %size%$$. 00075 The input value of the elements of $icode x$$ does not matter. 00076 On output, it has been set to the 00077 argument value for which the function, 00078 or its derivative, is being evaluated. 00079 The value of this vector need not change with each repetition. 00080 00081 $head row$$ 00082 The argument $icode row$$ has prototype 00083 $codei% 00084 const CppAD::vector<size_t> %row% 00085 %$$ 00086 Its size defines the value $latex K$$. 00087 It contains the row indices for the corresponding function $latex f(x)$$. 00088 All the elements of $icode row$$ are between zero and $latex n-1$$. 00089 00090 $head col$$ 00091 The argument $icode col$$ has prototype 00092 $codei% 00093 const CppAD::vector<size_t> %col% 00094 %$$ 00095 Its size must be the same as $icode row$$; i.e., $latex K$$. 00096 It contains the column indices for the corresponding function 00097 $latex f(x)$$. 00098 All the elements of $icode col$$ are between zero and $latex n-1$$. 00099 $pre 00100 00101 $$ 00102 There are no duplicate row and column entires; i.e., if $icode%j% != %k%$$, 00103 $codei% 00104 %row%[%j%] != %row%[%k%] || %col%[%j%] != %col%[%k%] 00105 %$$ 00106 Only the lower triangle of the Hessian is included in the indices; i.e., 00107 $codei%row%[%k%] >= %col%[%k%]%$$. 00108 Furthermore, for all the non-zero entries in the lower triangle 00109 are included; i.e., if $latex i \geq j$$ and 00110 $latex \DD{f}{x[i]}{x[j]} \neq 0$$, 00111 there is an index $icode k$$ such that 00112 $icode%i% = %row%[%k%]%$$ and 00113 $icode%j% = %col%[%k%]%$$. 00114 00115 $head hessian$$ 00116 The argument $icode hessian$$ has prototype 00117 $codei% 00118 CppAD::vector<double> &hessian 00119 %$$ 00120 and its size is $latex n \times n$$. 00121 The input value of its elements does not matter. 00122 The output value of its elements is the Hessian of the function $latex f(x)$$. 00123 To be more specific, for 00124 $latex i = 0 , \ldots , n-1$$, 00125 $latex j = 0 , \ldots , n-1$$, 00126 $latex \[ 00127 \DD{f}{x[i]}{x[j]} (x) = hessian [ i * n + j ] 00128 \] $$ 00129 00130 $subhead double$$ 00131 In the case where $icode package$$ is $code double$$, 00132 only the first element of $icode hessian$$ is used and it is actually 00133 the value of $latex f(x)$$ (derivatives are not computed). 00134 00135 $end 00136 ----------------------------------------------------------------------------- 00137 */ 00138 # include <cppad/vector.hpp> 00139 # include <cppad/near_equal.hpp> 00140 # include <cppad/speed/sparse_hes_fun.hpp> 00141 # include <cppad/speed/uniform_01.hpp> 00142 # include <cppad/index_sort.hpp> 00143 00144 /*! 00145 \defgroup link_sparse_hessian_cpp link_sparse_hessian.cpp 00146 \{ 00147 \file link_sparse_hessian.cpp 00148 Defines and implement sparse Hessian speed link to package specific code. 00149 */ 00150 namespace { 00151 using CppAD::vector; 00152 00153 /*! 00154 Class used by choose_row_col to determin order of row and column indices 00155 */ 00156 class Key { 00157 public: 00158 /// row index 00159 size_t row_; 00160 /// column index 00161 size_t col_; 00162 /// default constructor 00163 Key(void) 00164 { } 00165 /*! 00166 Construct from a value for row and col 00167 00168 \param row 00169 row value for this key 00170 00171 \param col 00172 column value for this key 00173 */ 00174 Key(size_t row, size_t col) 00175 : row_(row), col_(col) 00176 { } 00177 /*! 00178 Compare this key with another key using < operator 00179 00180 \param other 00181 the other key. 00182 */ 00183 bool operator<(const Key& other) const 00184 { if( row_ == other.row_ ) 00185 return col_ < other.col_; 00186 return row_ < other.row_; 00187 } 00188 }; 00189 00190 /*! 00191 Function that randomly choose the row and column indices 00192 00193 \param n [in] 00194 is the dimension of the argument space for the function f(x). 00195 00196 \param row [out] 00197 the input size and elements of \c row do not matter. 00198 Upon return it is the chosen row indices. 00199 00200 \param col [out] 00201 the input size and elements of \c col do not matter. 00202 Upon return it is the chosen column indices. 00203 */ 00204 void choose_row_col( 00205 size_t n , 00206 vector<size_t>& row , 00207 vector<size_t>& col ) 00208 { size_t r, c, k, K = 5 * n; 00209 00210 // get the random indices 00211 vector<double> random(2 * K); 00212 CppAD::uniform_01(2 * K, random); 00213 00214 // sort the temporary row and colunn choices 00215 vector<Key> keys(K + n); 00216 vector<size_t> ind(K + n); 00217 for(k = 0; k < K; k++) 00218 { r = size_t( n * random[k] ); 00219 r = std::min(n-1, r); 00220 // 00221 c = size_t( n * random[k + K] ); 00222 c = std::min(n-1, c); 00223 // 00224 // force to lower triangle 00225 if( c > r ) 00226 std::swap(r, c); 00227 // 00228 keys[k] = Key(r, c); 00229 } 00230 // include the diagonal 00231 for(k = 0; k < n; k++) 00232 keys[k + K] = Key(k, k); 00233 CppAD::index_sort(keys, ind); 00234 00235 // remove duplicates while setting the return value for row and col 00236 row.resize(0); 00237 col.resize(0); 00238 size_t r_previous = keys[ ind[0] ].row_; 00239 size_t c_previous = keys[ ind[0] ].col_; 00240 row.push_back(r_previous); 00241 col.push_back(c_previous); 00242 for(k = 1; k < K; k++) 00243 { r = keys[ ind[k] ].row_; 00244 c = keys[ ind[k] ].row_; 00245 if( r != r_previous || c != c_previous) 00246 { row.push_back(r); 00247 col.push_back(c); 00248 } 00249 r_previous = r; 00250 c_previous = c; 00251 } 00252 } 00253 } 00254 00255 /*! 00256 Package specific implementation of a sparse Hessian claculation. 00257 00258 \param size [in] 00259 is the size of the domain space; i.e. specifies \c n. 00260 00261 \param repeat [in] 00262 number of times tha the test is repeated. 00263 00264 \param x [out] 00265 is a vector of size \c n containing 00266 the argument at which the Hessian was evaluated during the last repetition. 00267 00268 \param row [in] 00269 is the row indices correpsonding to non-zero Hessian entries. 00270 00271 \param col [in] 00272 is the column indices corresponding to non-zero Hessian entries. 00273 00274 \param hessian [out] 00275 is a vector with size <code>n * n</code> 00276 containing the value of the Hessian of f(x) 00277 corresponding to the last repetition. 00278 00279 \return 00280 is true, if the sparse Hessian speed test is implemented for this package, 00281 and false otherwise. 00282 */ 00283 extern bool link_sparse_hessian( 00284 size_t size , 00285 size_t repeat , 00286 CppAD::vector<double> &x , 00287 const CppAD::vector<size_t> &row , 00288 const CppAD::vector<size_t> &col , 00289 CppAD::vector<double> &hessian 00290 ); 00291 00292 /*! 00293 Is sparse Hessian test avaialable. 00294 00295 \return 00296 true, if spare Hessian available for this package, and false otherwise. 00297 */ 00298 bool available_sparse_hessian(void) 00299 { 00300 size_t n = 2; 00301 size_t repeat = 1; 00302 vector<double> x(n); 00303 vector<double> hessian(n * n); 00304 vector<size_t> row, col; 00305 choose_row_col(n, row, col); 00306 00307 return link_sparse_hessian(n, repeat, x, row, col, hessian); 00308 exit(0); 00309 } 00310 /*! 00311 Does final sparse Hessian value pass correctness test. 00312 00313 \return 00314 true, if correctness test passes, and false otherwise. 00315 */ 00316 bool correct_sparse_hessian(bool is_package_double) 00317 { 00318 size_t n = 10; 00319 size_t repeat = 1; 00320 vector<double> x(n); 00321 vector<double> hessian(n * n); 00322 vector<size_t> row, col; 00323 choose_row_col(n, row, col); 00324 00325 link_sparse_hessian(n, repeat, x, row, col, hessian); 00326 00327 size_t order, size; 00328 if( is_package_double) 00329 { order = 0; // check function value 00330 size = 1; 00331 } 00332 else 00333 { order = 2; // check hessian value 00334 size = n * n; 00335 } 00336 CppAD::vector<double> check(size); 00337 CppAD::sparse_hes_fun<double>(n, x, row, col, order, check); 00338 bool ok = true; 00339 size_t k; 00340 for( k = 0; k < size; k++) 00341 ok &= CppAD::NearEqual(check[k], hessian[k], 1e-10, 1e-10); 00342 00343 return ok; 00344 } 00345 00346 /*! 00347 Sparse Hessian speed test. 00348 00349 \param size 00350 is the dimension of the argument space for this speed test. 00351 00352 \param repeat 00353 is the number of times to repeate the speed test. 00354 */ 00355 void speed_sparse_hessian(size_t size, size_t repeat) 00356 { size_t n = size; 00357 vector<double> x(n); 00358 vector<double> hessian(n * n); 00359 vector<size_t> row, col; 00360 choose_row_col(n, row, col); 00361 00362 // note that cppad/sparse_hessian.cpp assumes that x.size() == size 00363 link_sparse_hessian(n, repeat, x, row, col, hessian); 00364 return; 00365 } 00366 00367 /*! \} */