KJB
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
diff_hessian_ind.h
Go to the documentation of this file.
1 /* =========================================================================== *
2  |
3  | Copyright (c) 1994-2011 by Kobus Barnard (author)
4  |
5  | Personal and educational use of this code is granted, provided that this
6  | header is kept intact, and that the authorship is not misrepresented, that
7  | its use is acknowledged in publications, and relevant papers are cited.
8  |
9  | For other use contact the author (kobus AT cs DOT arizona DOT edu).
10  |
11  | Please note that the code in this file has not necessarily been adequately
12  | tested. Naturally, there is no guarantee of performance, support, or fitness
13  | for any particular task. Nonetheless, I am interested in hearing about
14  | problems that you encounter.
15  |
16  | Author: Ernesto Brau
17  * =========================================================================== */
18 
19 /* $Id$ */
20 
21 #ifndef DIFF_HESSIAN_IND_H
22 #define DIFF_HESSIAN_IND_H
23 
24 #include <m_cpp/m_matrix.h>
25 #include <diff_cpp/diff_util.h>
26 #include <vector>
27 
28 namespace kjb {
29 
47 template<class Func, class Model, class Adapter>
48 Matrix hessian_ind
49 (
50  const Func& f,
51  const Model& x,
52  const std::vector<double>& dx,
53  const Adapter& adapter
54 )
55 {
56  const size_t D = adapter.size(&x);
57 
58  Matrix H(D, D);
59  Model y = x;
60 
61  for(size_t i = 0; i < D; i++)
62  {
63  for(size_t j = 0; j < D; j++)
64  {
65  if(i == j)
66  {
67  double yi = adapter.get(&y, i);
68 
69  // compute at current spot
70  double fx = f(y, i, i);
71 
72  // move right and compute
73  move_param(y, i, dx[i], adapter);
74  double fxp = f(y, i, i);
75 
76  // move left and compute
77  move_param(y, i, -2*dx[i], adapter);
78  double fxm = f(y, i, i);
79 
80  // hessian
81  H(i, i) = (fxp - 2*fx + fxm) / (dx[i]*dx[i]);
82 
83  // move back to original spot
84  adapter.set(&y, i, yi);
85  }
86  else
87  {
88  double yi = adapter.get(&y, i);
89  double yj = adapter.get(&y, j);
90 
91  // move i and j ++ and compute
92  move_params(y, i, j, dx[i], dx[j], adapter);
93  double fxpp = f(y, i, j);
94 
95  // move i and j +- and compute
96  move_param(y, j, -2*dx[j], adapter);
97  double fxpm = f(y, i, j);
98 
99  // move i and j -- and compute
100  move_param(y, i, -2*dx[i], adapter);
101  double fxmm = f(y, i, j);
102 
103  // move i and j -+ and compute
104  move_param(y, j, 2*dx[j], adapter);
105  double fxmp = f(y, i, j);
106 
107  H(i, j) = (fxpp - fxpm - fxmp + fxmm) / (4*dx[i]*dx[j]);
108 
109  // move back
110  adapter.set(&y, i, j, yi, yj);
111  }
112  }
113  }
114 
115  return H;
116 }
117 
127 template<class Func, class Vec>
128 inline
129 Matrix hessian_ind
130 (
131  const Func& f,
132  const Vec& x,
133  const std::vector<double>& dx
134 )
135 {
136  return hessian_ind(f, x, dx, Vector_adapter<Vec>());
137 }
138 
155 template<class Func, class Model, class Adapter>
157 (
158  const Func& f,
159  const Model& x,
160  const std::vector<double>& dx,
161  const Adapter& adapter
162 )
163 {
164  const size_t D = adapter.size(&x);
165 
166  Matrix H(D, D);
167  Model y = x;
168 
169  for(size_t i = 0; i < D; i++)
170  {
171  for(size_t j = 0; j <= i; j++)
172  {
173  if(i == j)
174  {
175  double yi = adapter.get(&y, i);
176 
177  // compute at current spot
178  double fx = f(y, i, i);
179 
180  // move right and compute
181  move_param(y, i, dx[i], adapter);
182  double fxp = f(y, i, i);
183 
184  // move left and compute
185  move_param(y, i, -2*dx[i], adapter);
186  double fxm = f(y, i, i);
187 
188  // hessian
189  H(i, i) = (fxp - 2*fx + fxm) / (dx[i]*dx[i]);
190 
191  // move back to original spot
192  adapter.set(&y, i, yi);
193  }
194  else
195  {
196  double yi = adapter.get(&y, i);
197  double yj = adapter.get(&y, j);
198 
199  // move i and j ++ and compute
200  move_params(y, i, j, dx[i], dx[j], adapter);
201  double fxpp = f(y, i, j);
202 
203  // move i and j +- and compute
204  move_param(y, j, -2*dx[j], adapter);
205  double fxpm = f(y, i, j);
206 
207  // move i and j -- and compute
208  move_param(y, i, -2*dx[i], adapter);
209  double fxmm = f(y, i, j);
210 
211  // move i and j -+ and compute
212  move_param(y, j, 2*dx[j], adapter);
213  double fxmp = f(y, i, j);
214 
215  H(i, j) = (fxpp - fxpm - fxmp + fxmm) / (4*dx[i]*dx[j]);
216  H(j, i) = H(i, j);
217 
218  // move back
219  adapter.set(&y, i, j, yi, yj);
220  }
221  }
222  }
223 
224  return H;
225 }
226 
236 template<class Func, class Vec>
237 inline
239 (
240  const Func& f,
241  const Vec& x,
242  const std::vector<double>& dx
243 )
244 {
245  return hessian_symmetric_ind(f, x, dx, Vector_adapter<Vec>());
246 }
247 
265 template<class Func, class Model, class Adapter>
267 (
268  const Func& f,
269  const Model& x,
270  const std::vector<double>& dx,
271  const Adapter& adapter,
272  size_t is,
273  size_t ie
274 )
275 {
276  const size_t D = adapter.size(&x);
277 
278  IFT(is < D && ie < D && ie > is, Runtime_error,
279  "Cannot compute Hessian diagonal; bad indices.");
280 
281  Vector H(ie - is + 1);
282  Model y = x;
283 
284  for(size_t i = is; i <= ie; i++)
285  {
286  double yi = adapter.get(&y, i);
287 
288  // compute at current spot
289  double fx = f(y, i);
290 
291  // move right and compute
292  move_param(y, i, dx[i], adapter);
293  double fxp = f(y, i);
294 
295  // move left and compute
296  move_param(y, i, -2*dx[i], adapter);
297  double fxm = f(y, i);
298 
299  // hessian
300  H[i - is] = (fxp - 2*fx + fxm) / (dx[i]*dx[i]);
301 
302  // move back to original spot
303  adapter.set(&y, i, yi);
304  }
305 
306  return H;
307 }
308 
318 template<class Func, class Vec>
319 inline
321 (
322  const Func& f,
323  const Vec& x,
324  const std::vector<double>& dx,
325  size_t is,
326  size_t ie
327 )
328 {
329  return hessian_ind_diagonal(f, x, dx, Vector_adapter<Vec>(), is, ie);
330 }
331 
332 } //namespace kjb
333 
334 #endif /*DIFF_HESSIAN_H_IND */
335 
Definition for the Matrix class, a thin wrapper on the KJB Matrix struct and its related functionalit...
void move_param(Model &x, size_t i, double dv, const Adapter &aptr)
Helper function that moves a parameter by an amount.
Definition: diff_util.h:102
This class implements vectors, in the linear-algebra sense, with real-valued elements.
Definition: m_vector.h:87
#define IFT(a, ex, msg)
Definition: l_exception.h:101
Matrix hessian_symmetric_ind(const Func &f, const Model &x, const std::vector< double > &dx, const Adapter &adapter)
Computes the Hessian of an "independent" function, evaluated at a point, using finite differences...
Definition: diff_hessian_ind.h:157
void move_params(Model &x, size_t i, size_t j, double dv, double dw, const Adapter &aptr)
Helper function that moves a pair of parameters by an amount.
Definition: diff_util.h:111
x
Definition: APPgetLargeConnectedEdges.m:100
Vector & set(Value_type val)
Clone of zero_out(int)
Definition: m_vector.h:707
Matrix hessian_ind(const Func &f, const Model &x, const std::vector< double > &dx, const Adapter &adapter)
Computes the Hessian of a function, evaluated at a point, using finite differences.
Definition: diff_hessian_ind.h:49
Default adapter for the hessian function.
Definition: diff_util.h:42
get the indices of edges in each direction for i
Definition: APPgetLargeConnectedEdges.m:48
This class implements matrices, in the linear-algebra sense, with real-valued elements.
Definition: m_matrix.h:94
D
Definition: APPgetLargeConnectedEdges.m:106
Vector hessian_ind_diagonal(const Func &f, const Model &x, const std::vector< double > &dx, const Adapter &adapter, size_t is, size_t ie)
Computes the Hessian diagonal of a function, evaluated at a point, using finite differences.
Definition: diff_hessian_ind.h:267
Object thrown when computation fails somehow during execution.
Definition: l_exception.h:321