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"
55#include "ROL_Secant.hpp"
56#include "Teuchos_SerialDenseMatrix.hpp"
57#include "Teuchos_SerialDenseVector.hpp"
58#include "Teuchos_LAPACK.hpp"
59
60namespace 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 ROL::Ptr<Vector<Real> > e = x.clone();
71 ROL::Ptr<Vector<Real> > h = x.dual().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->dual());
79 H(j,i) = e->apply(*h);
80 }
81 }
82
83 return H;
84
85 }
86
87
88 template<class Real>
89 Teuchos::SerialDenseMatrix<int, Real> computeScaledDenseHessian(Objective<Real> &obj, const Vector<Real> &x) {
90
91 Real tol = std::sqrt(ROL_EPSILON<Real>());
92
93 int dim = x.dimension();
94 Teuchos::SerialDenseMatrix<int, Real> H(dim, dim);
95
96 ROL::Ptr<Vector<Real> > ei = x.clone();
97 ROL::Ptr<Vector<Real> > ej = x.dual().clone();
98 ROL::Ptr<Vector<Real> > h = x.dual().clone();
99
100 for (int i=0; i<dim; i++) {
101 ei = ei->basis(i);
102 obj.hessVec(*h, *ei, x, tol);
103 for (int j=0; j<dim; j++) {
104 ej = ej->basis(j);
105 H(j,i) = ej->dot(*h);
106 }
107 }
108
109 return H;
110
111 }
112
113
114 template<class Real>
115 Teuchos::SerialDenseMatrix<int, Real> computeDotMatrix(const Vector<Real> &x) {
116
117 int dim = x.dimension();
118 Teuchos::SerialDenseMatrix<int, Real> M(dim, dim);
119
120 ROL::Ptr<Vector<Real> > ei = x.clone();
121 ROL::Ptr<Vector<Real> > ej = x.clone();
122
123 for (int i=0; i<dim; i++) {
124 ei = x.basis(i);
125 for (int j=0; j<dim; j++) {
126 ej = x.basis(j);
127 M(j,i) = ej->dot(*ei);
128 }
129 }
130
131 return M;
132
133 }
134
135 template<class Real>
136 std::vector<std::vector<Real> > computeEigenvalues(const Teuchos::SerialDenseMatrix<int, Real> & mat) {
137
138 Teuchos::LAPACK<int, Real> lapack;
139
140 Teuchos::SerialDenseMatrix<int, Real> mymat(Teuchos::Copy, mat);
141
142 char jobvl = 'N';
143 char jobvr = 'N';
144
145 int n = mat.numRows();
146
147 std::vector<Real> real(n, 0);
148 std::vector<Real> imag(n, 0);
149 std::vector<std::vector<Real> > eigenvals;
150
151 Real* vl = 0;
152 Real* vr = 0;
153
154 int ldvl = 1;
155 int ldvr = 1;
156
157 int lwork = 4*n;
158
159 std::vector<Real> work(lwork, 0);
160
161 int info = 0;
162
163 lapack.GEEV(jobvl, jobvr, n, &mymat(0,0), n, &real[0], &imag[0], vl, ldvl, vr, ldvr, &work[0], lwork, &info);
164
165 eigenvals.push_back(real);
166 eigenvals.push_back(imag);
167
168 return eigenvals;
169
170 }
171
172
173 template<class Real>
174 std::vector<std::vector<Real> > computeGenEigenvalues(const Teuchos::SerialDenseMatrix<int, Real> & A,
175 const Teuchos::SerialDenseMatrix<int, Real> & B) {
176
177 Teuchos::LAPACK<int, Real> lapack;
178
179 Teuchos::SerialDenseMatrix<int, Real> myA(Teuchos::Copy, A);
180 Teuchos::SerialDenseMatrix<int, Real> myB(Teuchos::Copy, B);
181
182 char jobvl = 'N';
183 char jobvr = 'N';
184
185 int n = A.numRows();
186
187 std::vector<Real> real(n, 0);
188 std::vector<Real> imag(n, 0);
189 std::vector<Real> beta(n, 0);
190 std::vector<std::vector<Real> > eigenvals;
191
192 Real* vl = 0;
193 Real* vr = 0;
194
195 int ldvl = 1;
196 int ldvr = 1;
197
198 int lwork = 10*n;
199
200 std::vector<Real> work(lwork, 0);
201
202 int info = 0;
203
204 lapack.GGEV(jobvl, jobvr, n, &myA(0,0), n, &myB(0,0), n, &real[0], &imag[0], &beta[0],
205 vl, ldvl, vr, ldvr, &work[0], lwork, &info);
206
207 for (int i=0; i<n; i++) {
208 real[i] /= beta[i];
209 imag[i] /= beta[i];
210 }
211
212 eigenvals.push_back(real);
213 eigenvals.push_back(imag);
214
215 return eigenvals;
216
217 }
218
219
220 template<class Real>
221 Teuchos::SerialDenseMatrix<int, Real> computeInverse(const Teuchos::SerialDenseMatrix<int, Real> & mat) {
222
223 Teuchos::LAPACK<int, Real> lapack;
224
225 Teuchos::SerialDenseMatrix<int, Real> mymat(Teuchos::Copy, mat);
226
227 int n = mat.numRows();
228
229 std::vector<int> ipiv(n, 0);
230
231 int lwork = 5*n;
232
233 std::vector<Real> work(lwork, 0);
234
235 int info = 0;
236
237 lapack.GETRF(n, n, &mymat(0,0), n, &ipiv[0], &info);
238 lapack.GETRI(n, &mymat(0,0), n, &ipiv[0], &work[0], lwork, &info);
239
240 return mymat;
241
242 }
243
244
245
246 template<class Real>
247 class ProjectedObjective : public Objective<Real> {
248 private:
249 ROL::Ptr<Objective<Real> > obj_;
250 ROL::Ptr<BoundConstraint<Real> > con_;
251 ROL::Ptr<Secant<Real> > secant_;
252
253 ROL::Ptr<ROL::Vector<Real> > primalV_;
254 ROL::Ptr<ROL::Vector<Real> > dualV_;
256
259 Real eps_;
260
261 public:
263 ROL::Ptr<Secant<Real> > &secant,
264 bool useSecantPrecond = false,
265 bool useSecantHessVec = false,
266 Real eps = 0.0 )
267 : isInitialized_(false), useSecantPrecond_(useSecantPrecond),
268 useSecantHessVec_(useSecantHessVec), eps_(eps) {
269 obj_ = ROL::makePtrFromRef(obj);
270 con_ = ROL::makePtrFromRef(con);
271 secant_ = secant;
272 }
273
274 void update( const Vector<Real> &x, bool flag = true, int iter = -1 ) {
275 obj_->update(x,flag,iter);
276 con_->update(x,flag,iter);
277 }
278
279 Real value( const Vector<Real> &x, Real &tol ) {
280 return obj_->value(x,tol);
281 }
282
283 void gradient( Vector<Real> &g, const Vector<Real> &x, Real &tol ) {
284 obj_->gradient(g,x,tol);
285 }
286
287 Real dirDeriv( const Vector<Real> &x, const Vector<Real> &d, Real &tol ) {
288 return obj_->dirDeriv(x,d,tol);
289 }
290
291 void hessVec( Vector<Real> &Hv, const Vector<Real> &v, const Vector<Real> &x, Real &tol ) {
292 if ( useSecantHessVec_ ) {
293 secant_->applyB( Hv, v );
294 }
295 else {
296 obj_->hessVec( Hv, v, x, tol );
297 }
298 }
299
300 void invHessVec( Vector<Real> &Hv, const Vector<Real> &v, const Vector<Real> &x, Real &tol ) {
301 if ( useSecantHessVec_ ) {
302 secant_->applyH(Hv,v);
303 }
304 else {
305 obj_->invHessVec(Hv,v,x,tol);
306 }
307 }
308
309 void precond( Vector<Real> &Mv, const Vector<Real> &v, const Vector<Real> &x, Real &tol ) {
310 if ( useSecantPrecond_ ) {
311 secant_->applyH( Mv, v );
312 }
313 else {
314 obj_->precond( Mv, v, x, tol );
315 }
316 }
317
330 const Vector<Real> &d, const Vector<Real> &x, Real &tol ) {
331 if ( con_->isActivated() ) {
332 if (!isInitialized_) {
333 primalV_ = x.clone();
334 dualV_ = x.dual().clone();
335 isInitialized_ = true;
336 }
337 // Set vnew to v
338 primalV_->set(v);
339 // Remove elements of vnew corresponding to binding set
340 con_->pruneActive(*primalV_,d,p,eps_);
341 // Apply full Hessian to reduced vector
342 hessVec(Hv,*primalV_,x,tol);
343 // Remove elements of Hv corresponding to binding set
344 con_->pruneActive(Hv,d,p,eps_);
345 // Set vnew to v
346 primalV_->set(v);
347 // Remove Elements of vnew corresponding to complement of binding set
348 con_->pruneInactive(*primalV_,d,p,eps_);
349 dualV_->set(primalV_->dual());
350 con_->pruneInactive(*dualV_,d,p,eps_);
351 // Fill complement of binding set elements in Hp with v
352 Hv.plus(*dualV_);
353 }
354 else {
355 hessVec(Hv,v,x,tol);
356 }
357 }
358
370 const Vector<Real> &x, Real &tol ) {
371 if ( con_->isActivated() ) {
372 if (!isInitialized_) {
373 primalV_ = x.clone();
374 dualV_ = x.dual().clone();
375 isInitialized_ = true;
376 }
377 // Set vnew to v
378 primalV_->set(v);
379 // Remove elements of vnew corresponding to binding set
380 con_->pruneActive(*primalV_,p,eps_);
381 // Apply full Hessian to reduced vector
382 hessVec(Hv,*primalV_,x,tol);
383 // Remove elements of Hv corresponding to binding set
384 con_->pruneActive(Hv,p,eps_);
385 // Set vnew to v
386 primalV_->set(v);
387 // Remove Elements of vnew corresponding to complement of binding set
388 con_->pruneInactive(*primalV_,p,eps_);
389 dualV_->set(primalV_->dual());
390 con_->pruneInactive(*dualV_,p,eps_);
391 // Fill complement of binding set elements in Hp with v
392 Hv.plus(*dualV_);
393 }
394 else {
395 hessVec(Hv,v,x,tol);
396 }
397 }
398
411 const Vector<Real> &d, const Vector<Real> &x, Real &tol ) {
412 if ( con_->isActivated() ) {
413 if (!isInitialized_) {
414 primalV_ = x.clone();
415 dualV_ = x.dual().clone();
416 isInitialized_ = true;
417 }
418 // Set vnew to v
419 dualV_->set(v);
420 // Remove elements of vnew corresponding to binding set
421 con_->pruneActive(*dualV_,d,p,eps_);
422 // Apply full Hessian to reduced vector
423 invHessVec(Hv,*dualV_,x,tol);
424 // Remove elements of Hv corresponding to binding set
425 con_->pruneActive(Hv,d,p,eps_);
426 // Set vnew to v
427 dualV_->set(v);
428 // Remove Elements of vnew corresponding to complement of binding set
429 con_->pruneInactive(*dualV_,d,p,eps_);
430 primalV_->set(dualV_->dual());
431 con_->pruneInactive(*primalV_,d,p,eps_);
432 // Fill complement of binding set elements in Hv with v
433 Hv.plus(*primalV_);
434 }
435 else {
436 invHessVec(Hv,v,x,tol);
437 }
438 }
439
451 const Vector<Real> &x, Real &tol ) {
452 if ( con_->isActivated() ) {
453 if (!isInitialized_) {
454 primalV_ = x.clone();
455 dualV_ = x.dual().clone();
456 isInitialized_ = true;
457 }
458 // Set vnew to v
459 dualV_->set(v);
460 // Remove elements of vnew corresponding to binding set
461 con_->pruneActive(*dualV_,p,eps_);
462 // Apply full Hessian to reduced vector
463 invHessVec(Hv,*dualV_,x,tol);
464 // Remove elements of Hv corresponding to binding set
465 con_->pruneActive(Hv,p,eps_);
466 // Set vnew to v
467 dualV_->set(v);
468 // Remove Elements of vnew corresponding to complement of binding set
469 con_->pruneInactive(*dualV_,p,eps_);
470 primalV_->set(dualV_->dual());
471 con_->pruneInactive(*primalV_,p,eps_);
472 // Fill complement of binding set elements in Hv with v
473 Hv.plus(*primalV_);
474 }
475 else {
476 invHessVec(Hv,v,x,tol);
477 }
478 }
479
491 void reducedPrecond( Vector<Real> &Mv, const Vector<Real> &v, const Vector<Real> &p,
492 const Vector<Real> &d, const Vector<Real> &x, Real &tol ) {
493 if ( con_->isActivated() ) {
494 if (!isInitialized_) {
495 primalV_ = x.clone();
496 dualV_ = x.dual().clone();
497 isInitialized_ = true;
498 }
499 // Set vnew to v
500 dualV_->set(v);
501 // Remove elements of vnew corresponding to binding set
502 con_->pruneActive(*dualV_,d,p,eps_);
503 // Apply full Hessian to reduced vector
504 precond(Mv,*dualV_,x,tol);
505 // Remove elements of Mv corresponding to binding set
506 con_->pruneActive(Mv,d,p,eps_);
507 // Set vnew to v
508 dualV_->set(v);
509 // Remove Elements of vnew corresponding to complement of binding set
510 con_->pruneInactive(*dualV_,d,p,eps_);
511 primalV_->set(dualV_->dual());
512 con_->pruneInactive(*primalV_,d,p,eps_);
513 // Fill complement of binding set elements in Mv with v
514 Mv.plus(*primalV_);
515 }
516 else {
517 precond(Mv,v,x,tol);
518 }
519 }
520
531 void reducedPrecond( Vector<Real> &Mv, const Vector<Real> &v, const Vector<Real> &p,
532 const Vector<Real> &x, Real &tol ) {
533 if ( con_->isActivated() ) {
534 if (!isInitialized_) {
535 primalV_ = x.clone();
536 dualV_ = x.dual().clone();
537 isInitialized_ = true;
538 }
539 // Set vnew to v
540 dualV_->set(v);
541 // Remove elements of vnew corresponding to binding set
542 con_->pruneActive(*dualV_,p,eps_);
543 // Apply full Hessian to reduced vector
544 precond(Mv,*dualV_,x,tol);
545 // Remove elements of Mv corresponding to binding set
546 con_->pruneActive(Mv,p,eps_);
547 // Set vnew to v
548 dualV_->set(v);
549 // Remove Elements of vnew corresponding to complement of binding set
550 con_->pruneInactive(*dualV_,p,eps_);
551 primalV_->set(dualV_->dual());
552 con_->pruneInactive(*primalV_,p,eps_);
553 // Fill complement of binding set elements in Mv with v
554 Mv.plus(*primalV_);
555 }
556 else {
557 precond(Mv,v,x,tol);
558 }
559 }
560
562 con_->project(x);
563 }
564
565 void pruneActive( Vector<Real> &v, const Vector<Real> &g, const Vector<Real> &x ) {
566 con_->pruneActive(v,g,x,eps_);
567 }
568
569 void pruneActive( Vector<Real> &v, const Vector<Real> &x ) {
570 con_->pruneActive(v,x,eps_);
571 }
572
573 void pruneInactive( Vector<Real> &v, const Vector<Real> &g, const Vector<Real> &x ) {
574 con_->pruneInactive(v,g,x,eps_);
575 }
576
578 con_->pruneInactive(v,x,eps_);
579 }
580
581 bool isFeasible( const Vector<Real> &v ) {
582 return con_->isFeasible(v);
583 }
584
585 bool isConActivated(void) {
586 return con_->isActivated();
587 }
588
590 con_->computeProjectedStep(v,x);
591 }
592 };
593
594} // namespace ROL
595
596#endif
Provides the interface to apply upper and lower bound constraints.
Provides the interface to evaluate objective functions.
virtual void hessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply Hessian approximation to vector.
void project(Vector< Real > &x)
ROL::Ptr< ROL::Vector< Real > > dualV_
Real value(const Vector< Real > &x, Real &tol)
Compute value.
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 gradient(Vector< Real > &g, const Vector< Real > &x, Real &tol)
Compute gradient.
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...
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...
ROL::Ptr< ROL::Vector< Real > > primalV_
bool isFeasible(const Vector< Real > &v)
void pruneActive(Vector< Real > &v, const Vector< Real > &g, const Vector< Real > &x)
ROL::Ptr< BoundConstraint< Real > > con_
void pruneActive(Vector< Real > &v, const Vector< Real > &x)
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...
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 ...
void precond(Vector< Real > &Mv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply preconditioner to vector.
ROL::Ptr< Secant< Real > > secant_
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...
void pruneInactive(Vector< Real > &v, const Vector< Real > &x)
void invHessVec(Vector< Real > &Hv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply inverse Hessian approximation to vector.
void computeProjectedStep(Vector< Real > &v, const Vector< Real > &x)
ProjectedObjective(Objective< Real > &obj, BoundConstraint< Real > &con, ROL::Ptr< Secant< Real > > &secant, bool useSecantPrecond=false, bool useSecantHessVec=false, Real eps=0.0)
void update(const Vector< Real > &x, bool flag=true, int iter=-1)
Update objective function.
void pruneInactive(Vector< Real > &v, const Vector< Real > &g, const Vector< Real > &x)
Real dirDeriv(const Vector< Real > &x, const Vector< Real > &d, Real &tol)
Compute directional derivative.
ROL::Ptr< Objective< Real > > obj_
void hessVec(Vector< Real > &Hv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply Hessian approximation to vector.
Provides interface for and implements limited-memory secant operators.
Defines the linear algebra or vector space interface.
virtual const Vector & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis,...
virtual void plus(const Vector &x)=0
Compute , where .
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
virtual int dimension() const
Return dimension of the vector space.
virtual ROL::Ptr< Vector > basis(const int i) const
Return i-th basis vector.
Teuchos::SerialDenseMatrix< int, Real > computeScaledDenseHessian(Objective< Real > &obj, const Vector< Real > &x)
Teuchos::SerialDenseMatrix< int, Real > computeDenseHessian(Objective< Real > &obj, const Vector< Real > &x)
Teuchos::SerialDenseMatrix< int, Real > computeDotMatrix(const Vector< Real > &x)
std::vector< std::vector< Real > > computeGenEigenvalues(const Teuchos::SerialDenseMatrix< int, Real > &A, const Teuchos::SerialDenseMatrix< int, Real > &B)
std::vector< std::vector< Real > > computeEigenvalues(const Teuchos::SerialDenseMatrix< int, Real > &mat)
Teuchos::SerialDenseMatrix< int, Real > computeInverse(const Teuchos::SerialDenseMatrix< int, Real > &mat)
constexpr auto dim