CppAD: A C++ Algorithmic Differentiation Package
20130918
|
00001 /* $Id$ */ 00002 /* -------------------------------------------------------------------------- 00003 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-14 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 \{ 00146 \file link_sparse_hessian.cpp 00147 Defines and implement sparse Hessian speed link to package specific code. 00148 */ 00149 namespace { 00150 using CppAD::vector; 00151 00152 /*! 00153 Class used by choose_row_col to determin order of row and column indices 00154 */ 00155 class Key { 00156 public: 00157 /// row index 00158 size_t row_; 00159 /// column index 00160 size_t col_; 00161 /// default constructor 00162 Key(void) 00163 { } 00164 /*! 00165 Construct from a value for row and col 00166 00167 \param row 00168 row value for this key 00169 00170 \param col 00171 column value for this key 00172 */ 00173 Key(size_t row, size_t col) 00174 : row_(row), col_(col) 00175 { } 00176 /*! 00177 Compare this key with another key using < operator 00178 00179 \param other 00180 the other key. 00181 */ 00182 bool operator<(const Key& other) const 00183 { if( row_ == other.row_ ) 00184 return col_ < other.col_; 00185 return row_ < other.row_; 00186 } 00187 }; 00188 00189 /*! 00190 Function that randomly choose the row and column indices 00191 00192 \param n [in] 00193 is the dimension of the argument space for the function f(x). 00194 00195 \param row [out] 00196 the input size and elements of \c row do not matter. 00197 Upon return it is the chosen row indices. 00198 00199 \param col [out] 00200 the input size and elements of \c col do not matter. 00201 Upon return it is the chosen column indices. 00202 */ 00203 void choose_row_col( 00204 size_t n , 00205 vector<size_t>& row , 00206 vector<size_t>& col ) 00207 { size_t r, c, k, K = 5 * n; 00208 00209 // get the random indices 00210 vector<double> random(2 * K); 00211 CppAD::uniform_01(2 * K, random); 00212 00213 // sort the temporary row and colunn choices 00214 vector<Key> keys(K + n); 00215 vector<size_t> ind(K + n); 00216 for(k = 0; k < K; k++) 00217 { r = size_t( n * random[k] ); 00218 r = std::min(n-1, r); 00219 // 00220 c = size_t( n * random[k + K] ); 00221 c = std::min(n-1, c); 00222 // 00223 // force to lower triangle 00224 if( c > r ) 00225 std::swap(r, c); 00226 // 00227 keys[k] = Key(r, c); 00228 } 00229 // include the diagonal 00230 for(k = 0; k < n; k++) 00231 keys[k + K] = Key(k, k); 00232 CppAD::index_sort(keys, ind); 00233 00234 // remove duplicates while setting the return value for row and col 00235 row.resize(0); 00236 col.resize(0); 00237 size_t r_previous = keys[ ind[0] ].row_; 00238 size_t c_previous = keys[ ind[0] ].col_; 00239 row.push_back(r_previous); 00240 col.push_back(c_previous); 00241 for(k = 1; k < K; k++) 00242 { r = keys[ ind[k] ].row_; 00243 c = keys[ ind[k] ].row_; 00244 if( r != r_previous || c != c_previous) 00245 { row.push_back(r); 00246 col.push_back(c); 00247 } 00248 r_previous = r; 00249 c_previous = c; 00250 } 00251 } 00252 } 00253 00254 /*! 00255 Package specific implementation of a sparse Hessian claculation. 00256 00257 \param size [in] 00258 is the size of the domain space; i.e. specifies \c n. 00259 00260 \param repeat [in] 00261 number of times tha the test is repeated. 00262 00263 \param x [out] 00264 is a vector of size \c n containing 00265 the argument at which the Hessian was evaluated during the last repetition. 00266 00267 \param row [in] 00268 is the row indices correpsonding to non-zero Hessian entries. 00269 00270 \param col [in] 00271 is the column indices corresponding to non-zero Hessian entries. 00272 00273 \param hessian [out] 00274 is a vector with size <code>n * n</code> 00275 containing the value of the Hessian of f(x) 00276 corresponding to the last repetition. 00277 00278 \return 00279 is true, if the sparse Hessian speed test is implemented for this package, 00280 and false otherwise. 00281 */ 00282 extern bool link_sparse_hessian( 00283 size_t size , 00284 size_t repeat , 00285 const CppAD::vector<size_t>& row , 00286 const CppAD::vector<size_t>& col , 00287 CppAD::vector<double>& x , 00288 CppAD::vector<double>& hessian 00289 ); 00290 00291 /*! 00292 Is sparse Hessian test avaialable. 00293 00294 \return 00295 true, if spare Hessian available for this package, and false otherwise. 00296 */ 00297 bool available_sparse_hessian(void) 00298 { 00299 size_t n = 2; 00300 size_t repeat = 1; 00301 vector<double> x(n); 00302 vector<double> hessian(n * n); 00303 vector<size_t> row, col; 00304 choose_row_col(n, row, col); 00305 00306 return link_sparse_hessian(n, repeat, row, col, x, hessian); 00307 exit(0); 00308 } 00309 /*! 00310 Does final sparse Hessian value pass correctness test. 00311 00312 \param is_package_double 00313 if true, we are checking function values instead of derivatives. 00314 00315 \return 00316 true, if correctness test passes, and false otherwise. 00317 */ 00318 bool correct_sparse_hessian(bool is_package_double) 00319 { 00320 size_t n = 10; 00321 size_t repeat = 1; 00322 vector<double> x(n); 00323 vector<double> hessian(n * n); 00324 vector<size_t> row, col; 00325 choose_row_col(n, row, col); 00326 00327 link_sparse_hessian(n, repeat, row, col, x, hessian); 00328 00329 size_t order, size; 00330 if( is_package_double) 00331 { order = 0; // check function value 00332 size = 1; 00333 } 00334 else 00335 { order = 2; // check hessian value 00336 size = n * n; 00337 } 00338 CppAD::vector<double> check(size); 00339 CppAD::sparse_hes_fun<double>(n, x, row, col, order, check); 00340 bool ok = true; 00341 size_t k; 00342 for( k = 0; k < size; k++) 00343 ok &= CppAD::NearEqual(check[k], hessian[k], 1e-10, 1e-10); 00344 00345 return ok; 00346 } 00347 00348 /*! 00349 Sparse Hessian speed test. 00350 00351 \param size 00352 is the dimension of the argument space for this speed test. 00353 00354 \param repeat 00355 is the number of times to repeate the speed test. 00356 */ 00357 void speed_sparse_hessian(size_t size, size_t repeat) 00358 { size_t n = size; 00359 vector<double> x(n); 00360 vector<double> hessian(n * n); 00361 vector<size_t> row, col; 00362 choose_row_col(n, row, col); 00363 00364 // note that cppad/sparse_hessian.cpp assumes that x.size() == size 00365 link_sparse_hessian(n, repeat, row, col, x, hessian); 00366 return; 00367 } 00368