ROL
ROL_HelperFunctions.hpp
Go to the documentation of this file.
1 // @HEADER
2 // ************************************************************************
3 //
4 // Rapid Optimization Library (ROL) Package
5 // Copyright (2014) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact lead developers:
38 // Drew Kouri (dpkouri@sandia.gov) and
39 // Denis Ridzal (dridzal@sandia.gov)
40 //
41 // ************************************************************************
42 // @HEADER
43 
49 #ifndef ROL_HELPERFUNCTIONS_HPP
50 #define ROL_HELPERFUNCTIONS_HPP
51 
52 #include "ROL_Vector.hpp"
53 #include "ROL_Objective.hpp"
54 #include "ROL_BoundConstraint.hpp"
55 #include "ROL_Secant.hpp"
56 #include "Teuchos_SerialDenseMatrix.hpp"
57 #include "Teuchos_SerialDenseVector.hpp"
58 #include "Teuchos_LAPACK.hpp"
59 
60 namespace ROL {
61 
62  template<class Real>
63  Teuchos::SerialDenseMatrix<int, Real> computeDenseHessian(Objective<Real> &obj, const Vector<Real> &x) {
64 
65  Real tol = std::sqrt(ROL_EPSILON<Real>());
66 
67  int dim = x.dimension();
68  Teuchos::SerialDenseMatrix<int, Real> H(dim, dim);
69 
70  Teuchos::RCP<Vector<Real> > e = x.clone();
71  Teuchos::RCP<Vector<Real> > h = x.clone();
72 
73  for (int i=0; i<dim; i++) {
74  e = x.basis(i);
75  obj.hessVec(*h, *e, x, tol);
76  for (int j=0; j<dim; j++) {
77  e = x.basis(j);
78  H(j,i) = e->dot(*h);
79  }
80  }
81 
82  return H;
83 
84  }
85 
86 
87  template<class Real>
88  Teuchos::SerialDenseMatrix<int, Real> computeDotMatrix(const Vector<Real> &x) {
89 
90  int dim = x.dimension();
91  Teuchos::SerialDenseMatrix<int, Real> M(dim, dim);
92 
93  Teuchos::RCP<Vector<Real> > ei = x.clone();
94  Teuchos::RCP<Vector<Real> > ej = x.clone();
95 
96  for (int i=0; i<dim; i++) {
97  ei = x.basis(i);
98  for (int j=0; j<dim; j++) {
99  ej = x.basis(j);
100  M(j,i) = ej->dot(*ei);
101  }
102  }
103 
104  return M;
105 
106  }
107 
108  template<class Real>
109  std::vector<std::vector<Real> > computeEigenvalues(const Teuchos::SerialDenseMatrix<int, Real> & mat) {
110 
111  Teuchos::LAPACK<int, Real> lapack;
112 
113  Teuchos::SerialDenseMatrix<int, Real> mymat(Teuchos::Copy, mat);
114 
115  char jobvl = 'N';
116  char jobvr = 'N';
117 
118  int n = mat.numRows();
119 
120  std::vector<Real> real(n, 0);
121  std::vector<Real> imag(n, 0);
122  std::vector<std::vector<Real> > eigenvals;
123 
124  Real* vl = 0;
125  Real* vr = 0;
126 
127  int ldvl = 1;
128  int ldvr = 1;
129 
130  int lwork = 4*n;
131 
132  std::vector<Real> work(lwork, 0);
133 
134  int info = 0;
135 
136  lapack.GEEV(jobvl, jobvr, n, &mymat(0,0), n, &real[0], &imag[0], vl, ldvl, vr, ldvr, &work[0], lwork, &info);
137 
138  eigenvals.push_back(real);
139  eigenvals.push_back(imag);
140 
141  return eigenvals;
142 
143  }
144 
145 
146  template<class Real>
147  std::vector<std::vector<Real> > computeGenEigenvalues(const Teuchos::SerialDenseMatrix<int, Real> & A,
148  const Teuchos::SerialDenseMatrix<int, Real> & B) {
149 
150  Teuchos::LAPACK<int, Real> lapack;
151 
152  Teuchos::SerialDenseMatrix<int, Real> myA(Teuchos::Copy, A);
153  Teuchos::SerialDenseMatrix<int, Real> myB(Teuchos::Copy, B);
154 
155  char jobvl = 'N';
156  char jobvr = 'N';
157 
158  int n = A.numRows();
159 
160  std::vector<Real> real(n, 0);
161  std::vector<Real> imag(n, 0);
162  std::vector<Real> beta(n, 0);
163  std::vector<std::vector<Real> > eigenvals;
164 
165  Real* vl = 0;
166  Real* vr = 0;
167 
168  int ldvl = 1;
169  int ldvr = 1;
170 
171  int lwork = 10*n;
172 
173  std::vector<Real> work(lwork, 0);
174 
175  int info = 0;
176 
177  lapack.GGEV(jobvl, jobvr, n, &myA(0,0), n, &myB(0,0), n, &real[0], &imag[0], &beta[0],
178  vl, ldvl, vr, ldvr, &work[0], lwork, &info);
179 
180  for (int i=0; i<n; i++) {
181  real[i] /= beta[i];
182  imag[i] /= beta[i];
183  }
184 
185  eigenvals.push_back(real);
186  eigenvals.push_back(imag);
187 
188  return eigenvals;
189 
190  }
191 
192 
193  template<class Real>
194  Teuchos::SerialDenseMatrix<int, Real> computeInverse(const Teuchos::SerialDenseMatrix<int, Real> & mat) {
195 
196  Teuchos::LAPACK<int, Real> lapack;
197 
198  Teuchos::SerialDenseMatrix<int, Real> mymat(Teuchos::Copy, mat);
199 
200  int n = mat.numRows();
201 
202  std::vector<int> ipiv(n, 0);
203 
204  int lwork = 5*n;
205 
206  std::vector<Real> work(lwork, 0);
207 
208  int info = 0;
209 
210  lapack.GETRF(n, n, &mymat(0,0), n, &ipiv[0], &info);
211  lapack.GETRI(n, &mymat(0,0), n, &ipiv[0], &work[0], lwork, &info);
212 
213  return mymat;
214 
215  }
216 
217 
218 
219  template<class Real>
220  class ProjectedObjective : public Objective<Real> {
221  private:
222  Teuchos::RCP<Objective<Real> > obj_;
223  Teuchos::RCP<BoundConstraint<Real> > con_;
224  Teuchos::RCP<Secant<Real> > secant_;
225 
226  Teuchos::RCP<ROL::Vector<Real> > primalV_;
227  Teuchos::RCP<ROL::Vector<Real> > dualV_;
229 
232  Real eps_;
233 
234  public:
236  Teuchos::RCP<Secant<Real> > &secant,
237  bool useSecantPrecond = false,
238  bool useSecantHessVec = false,
239  Real eps = 0.0 )
240  : isInitialized_(false), useSecantPrecond_(useSecantPrecond),
241  useSecantHessVec_(useSecantHessVec), eps_(eps) {
242  obj_ = Teuchos::rcpFromRef(obj);
243  con_ = Teuchos::rcpFromRef(con);
244  secant_ = secant;
245  }
246 
247  void update( const Vector<Real> &x, bool flag = true, int iter = -1 ) {
248  obj_->update(x,flag,iter);
249  con_->update(x,flag,iter);
250  }
251 
252  Real value( const Vector<Real> &x, Real &tol ) {
253  return obj_->value(x,tol);
254  }
255 
256  void gradient( Vector<Real> &g, const Vector<Real> &x, Real &tol ) {
257  obj_->gradient(g,x,tol);
258  }
259 
260  Real dirDeriv( const Vector<Real> &x, const Vector<Real> &d, Real &tol ) {
261  return obj_->dirDeriv(x,d,tol);
262  }
263 
264  void hessVec( Vector<Real> &Hv, const Vector<Real> &v, const Vector<Real> &x, Real &tol ) {
265  if ( useSecantHessVec_ ) {
266  secant_->applyB( Hv, v );
267  }
268  else {
269  obj_->hessVec( Hv, v, x, tol );
270  }
271  }
272 
273  void invHessVec( Vector<Real> &Hv, const Vector<Real> &v, const Vector<Real> &x, Real &tol ) {
274  if ( useSecantHessVec_ ) {
275  secant_->applyH(Hv,v);
276  }
277  else {
278  obj_->invHessVec(Hv,v,x,tol);
279  }
280  }
281 
282  void precond( Vector<Real> &Mv, const Vector<Real> &v, const Vector<Real> &x, Real &tol ) {
283  if ( useSecantPrecond_ ) {
284  secant_->applyH( Mv, v );
285  }
286  else {
287  obj_->precond( Mv, v, x, tol );
288  }
289  }
290 
302  void reducedHessVec( Vector<Real> &Hv, const Vector<Real> &v, const Vector<Real> &p,
303  const Vector<Real> &d, const Vector<Real> &x, Real &tol ) {
304  if ( con_->isActivated() ) {
305  if (!isInitialized_) {
306  primalV_ = x.clone();
307  dualV_ = x.dual().clone();
308  isInitialized_ = true;
309  }
310  // Set vnew to v
311  primalV_->set(v);
312  // Remove elements of vnew corresponding to binding set
313  con_->pruneActive(*primalV_,d,p,eps_);
314  // Apply full Hessian to reduced vector
315  hessVec(Hv,*primalV_,x,tol);
316  // Remove elements of Hv corresponding to binding set
317  con_->pruneActive(Hv,d,p,eps_);
318  // Set vnew to v
319  primalV_->set(v);
320  // Remove Elements of vnew corresponding to complement of binding set
321  con_->pruneInactive(*primalV_,d,p,eps_);
322  dualV_->set(primalV_->dual());
323  con_->pruneInactive(*dualV_,d,p,eps_);
324  // Fill complement of binding set elements in Hp with v
325  Hv.plus(*dualV_);
326  }
327  else {
328  hessVec(Hv,v,x,tol);
329  }
330  }
331 
342  void reducedHessVec( Vector<Real> &Hv, const Vector<Real> &v, const Vector<Real> &p,
343  const Vector<Real> &x, Real &tol ) {
344  if ( con_->isActivated() ) {
345  if (!isInitialized_) {
346  primalV_ = x.clone();
347  dualV_ = x.dual().clone();
348  isInitialized_ = true;
349  }
350  // Set vnew to v
351  primalV_->set(v);
352  // Remove elements of vnew corresponding to binding set
353  con_->pruneActive(*primalV_,p,eps_);
354  // Apply full Hessian to reduced vector
355  hessVec(Hv,*primalV_,x,tol);
356  // Remove elements of Hv corresponding to binding set
357  con_->pruneActive(Hv,p,eps_);
358  // Set vnew to v
359  primalV_->set(v);
360  // Remove Elements of vnew corresponding to complement of binding set
361  con_->pruneInactive(*primalV_,p,eps_);
362  dualV_->set(primalV_->dual());
363  con_->pruneInactive(*dualV_,p,eps_);
364  // Fill complement of binding set elements in Hp with v
365  Hv.plus(*dualV_);
366  }
367  else {
368  hessVec(Hv,v,x,tol);
369  }
370  }
371 
383  void reducedInvHessVec( Vector<Real> &Hv, const Vector<Real> &v, const Vector<Real> &p,
384  const Vector<Real> &d, const Vector<Real> &x, Real &tol ) {
385  if ( con_->isActivated() ) {
386  if (!isInitialized_) {
387  primalV_ = x.clone();
388  dualV_ = x.dual().clone();
389  isInitialized_ = true;
390  }
391  // Set vnew to v
392  dualV_->set(v);
393  // Remove elements of vnew corresponding to binding set
394  con_->pruneActive(*dualV_,d,p,eps_);
395  // Apply full Hessian to reduced vector
396  invHessVec(Hv,*dualV_,x,tol);
397  // Remove elements of Hv corresponding to binding set
398  con_->pruneActive(Hv,d,p,eps_);
399  // Set vnew to v
400  dualV_->set(v);
401  // Remove Elements of vnew corresponding to complement of binding set
402  con_->pruneInactive(*dualV_,d,p,eps_);
403  primalV_->set(dualV_->dual());
404  con_->pruneInactive(*primalV_,d,p,eps_);
405  // Fill complement of binding set elements in Hv with v
406  Hv.plus(*primalV_);
407  }
408  else {
409  invHessVec(Hv,v,x,tol);
410  }
411  }
412 
423  void reducedInvHessVec( Vector<Real> &Hv, const Vector<Real> &v, const Vector<Real> &p,
424  const Vector<Real> &x, Real &tol ) {
425  if ( con_->isActivated() ) {
426  if (!isInitialized_) {
427  primalV_ = x.clone();
428  dualV_ = x.dual().clone();
429  isInitialized_ = true;
430  }
431  // Set vnew to v
432  dualV_->set(v);
433  // Remove elements of vnew corresponding to binding set
434  con_->pruneActive(*dualV_,p,eps_);
435  // Apply full Hessian to reduced vector
436  invHessVec(Hv,*dualV_,x,tol);
437  // Remove elements of Hv corresponding to binding set
438  con_->pruneActive(Hv,p,eps_);
439  // Set vnew to v
440  dualV_->set(v);
441  // Remove Elements of vnew corresponding to complement of binding set
442  con_->pruneInactive(*dualV_,p,eps_);
443  primalV_->set(dualV_->dual());
444  con_->pruneInactive(*primalV_,p,eps_);
445  // Fill complement of binding set elements in Hv with v
446  Hv.plus(*primalV_);
447  }
448  else {
449  invHessVec(Hv,v,x,tol);
450  }
451  }
452 
464  void reducedPrecond( Vector<Real> &Mv, const Vector<Real> &v, const Vector<Real> &p,
465  const Vector<Real> &d, const Vector<Real> &x, Real &tol ) {
466  if ( con_->isActivated() ) {
467  if (!isInitialized_) {
468  primalV_ = x.clone();
469  dualV_ = x.dual().clone();
470  isInitialized_ = true;
471  }
472  // Set vnew to v
473  dualV_->set(v);
474  // Remove elements of vnew corresponding to binding set
475  con_->pruneActive(*dualV_,d,p,eps_);
476  // Apply full Hessian to reduced vector
477  precond(Mv,*dualV_,x,tol);
478  // Remove elements of Mv corresponding to binding set
479  con_->pruneActive(Mv,d,p,eps_);
480  // Set vnew to v
481  dualV_->set(v);
482  // Remove Elements of vnew corresponding to complement of binding set
483  con_->pruneInactive(*dualV_,d,p,eps_);
484  primalV_->set(dualV_->dual());
485  con_->pruneInactive(*primalV_,d,p,eps_);
486  // Fill complement of binding set elements in Mv with v
487  Mv.plus(*primalV_);
488  }
489  else {
490  precond(Mv,v,x,tol);
491  }
492  }
493 
504  void reducedPrecond( Vector<Real> &Mv, const Vector<Real> &v, const Vector<Real> &p,
505  const Vector<Real> &x, Real &tol ) {
506  if ( con_->isActivated() ) {
507  if (!isInitialized_) {
508  primalV_ = x.clone();
509  dualV_ = x.dual().clone();
510  isInitialized_ = true;
511  }
512  // Set vnew to v
513  dualV_->set(v);
514  // Remove elements of vnew corresponding to binding set
515  con_->pruneActive(*dualV_,p,eps_);
516  // Apply full Hessian to reduced vector
517  precond(Mv,*dualV_,x,tol);
518  // Remove elements of Mv corresponding to binding set
519  con_->pruneActive(Mv,p,eps_);
520  // Set vnew to v
521  dualV_->set(v);
522  // Remove Elements of vnew corresponding to complement of binding set
523  con_->pruneInactive(*dualV_,p,eps_);
524  primalV_->set(dualV_->dual());
525  con_->pruneInactive(*primalV_,p,eps_);
526  // Fill complement of binding set elements in Mv with v
527  Mv.plus(*primalV_);
528  }
529  else {
530  precond(Mv,v,x,tol);
531  }
532  }
533 
534  void project( Vector<Real> &x ) {
535  con_->project(x);
536  }
537 
538  void pruneActive( Vector<Real> &v, const Vector<Real> &g, const Vector<Real> &x ) {
539  con_->pruneActive(v,g,x,eps_);
540  }
541 
542  void pruneActive( Vector<Real> &v, const Vector<Real> &x ) {
543  con_->pruneActive(v,x,eps_);
544  }
545 
546  void pruneInactive( Vector<Real> &v, const Vector<Real> &g, const Vector<Real> &x ) {
547  con_->pruneInactive(v,g,x,eps_);
548  }
549 
550  void pruneInactive( Vector<Real> &v, const Vector<Real> &x ) {
551  con_->pruneInactive(v,x,eps_);
552  }
553 
554  bool isFeasible( const Vector<Real> &v ) {
555  return con_->isFeasible(v);
556  }
557 
558  bool isConActivated(void) {
559  return con_->isActivated();
560  }
561 
563  con_->computeProjectedStep(v,x);
564  }
565  };
566 
567 } // namespace ROL
568 
569 #endif
Provides the interface to evaluate objective functions.
Real value(const Vector< Real > &x, Real &tol)
Compute value.
Teuchos::RCP< ROL::Vector< Real > > primalV_
virtual void plus(const Vector &x)=0
Compute , where .
void reducedHessVec(Vector< Real > &Hv, const Vector< Real > &v, const Vector< Real > &p, const Vector< Real > &d, const Vector< Real > &x, Real &tol)
Apply the reduced Hessian to a vector, v. The reduced Hessian first removes elements of v correspondi...
virtual void hessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply Hessian approximation to vector.
Real dirDeriv(const Vector< Real > &x, const Vector< Real > &d, Real &tol)
Compute directional derivative.
void pruneInactive(Vector< Real > &v, const Vector< Real > &g, const Vector< Real > &x)
virtual Teuchos::RCP< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
Teuchos::RCP< BoundConstraint< Real > > con_
bool isFeasible(const Vector< Real > &v)
void pruneActive(Vector< Real > &v, const Vector< Real > &g, const Vector< Real > &x)
void reducedPrecond(Vector< Real > &Mv, const Vector< Real > &v, const Vector< Real > &p, const Vector< Real > &x, Real &tol)
Apply the reduced preconditioner to a vector, v. The reduced preconditioner first removes elements ...
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:74
void hessVec(Vector< Real > &Hv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply Hessian approximation to vector.
Teuchos::SerialDenseMatrix< int, Real > computeDenseHessian(Objective< Real > &obj, const Vector< Real > &x)
Teuchos::RCP< ROL::Vector< Real > > dualV_
std::vector< std::vector< Real > > computeGenEigenvalues(const Teuchos::SerialDenseMatrix< int, Real > &A, const Teuchos::SerialDenseMatrix< int, Real > &B)
virtual Teuchos::RCP< Vector > basis(const int i) const
Return i-th basis vector.
Definition: ROL_Vector.hpp:172
void reducedInvHessVec(Vector< Real > &Hv, const Vector< Real > &v, const Vector< Real > &p, const Vector< Real > &x, Real &tol)
Apply the reduced inverse Hessian to a vector, v. The reduced inverse Hessian first removes element...
virtual const Vector & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis...
Definition: ROL_Vector.hpp:213
ProjectedObjective(Objective< Real > &obj, BoundConstraint< Real > &con, Teuchos::RCP< Secant< Real > > &secant, bool useSecantPrecond=false, bool useSecantHessVec=false, Real eps=0.0)
virtual int dimension() const
Return dimension of the vector space.
Definition: ROL_Vector.hpp:183
Provides interface for and implements limited-memory secant operators.
Definition: ROL_Secant.hpp:70
Teuchos::RCP< Objective< Real > > obj_
std::vector< std::vector< Real > > computeEigenvalues(const Teuchos::SerialDenseMatrix< int, Real > &mat)
void reducedInvHessVec(Vector< Real > &Hv, const Vector< Real > &v, const Vector< Real > &p, const Vector< Real > &d, const Vector< Real > &x, Real &tol)
Apply the reduced inverse Hessian to a vector, v. The reduced inverse Hessian first removes element...
void pruneInactive(Vector< Real > &v, const Vector< Real > &x)
void reducedHessVec(Vector< Real > &Hv, const Vector< Real > &v, const Vector< Real > &p, const Vector< Real > &x, Real &tol)
Apply the reduced Hessian to a vector, v. The reduced Hessian first removes elements of v correspondi...
Provides the interface to apply upper and lower bound constraints.
void gradient(Vector< Real > &g, const Vector< Real > &x, Real &tol)
Compute gradient.
void update(const Vector< Real > &x, bool flag=true, int iter=-1)
Update objective function.
Teuchos::SerialDenseMatrix< int, Real > computeDotMatrix(const Vector< Real > &x)
void computeProjectedStep(Vector< Real > &v, const Vector< Real > &x)
void reducedPrecond(Vector< Real > &Mv, const Vector< Real > &v, const Vector< Real > &p, const Vector< Real > &d, const Vector< Real > &x, Real &tol)
Apply the reduced preconditioner to a vector, v. The reduced preconditioner first removes elements ...
void project(Vector< Real > &x)
void invHessVec(Vector< Real > &Hv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply inverse Hessian approximation to vector.
Teuchos::SerialDenseMatrix< int, Real > computeInverse(const Teuchos::SerialDenseMatrix< int, Real > &mat)
void precond(Vector< Real > &Mv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply preconditioner to vector.
void pruneActive(Vector< Real > &v, const Vector< Real > &x)
Teuchos::RCP< Secant< Real > > secant_