ROL
ROL_TruncatedCG_U.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
44#ifndef ROL_TRUNCATEDCG_U_H
45#define ROL_TRUNCATEDCG_U_H
46
51#include "ROL_TrustRegion_U.hpp"
52#include "ROL_Types.hpp"
53
54namespace ROL {
55
56template<class Real>
57class TruncatedCG_U : public TrustRegion_U<Real> {
58private:
59 Ptr<Vector<Real>> s_, g_, v_, p_, Hp_;
60
61 int maxit_;
62 Real tol1_;
63 Real tol2_;
64
65public:
66
67 // Constructor
68 TruncatedCG_U( ParameterList &parlist ) {
69 // Unravel Parameter List
70 Real em4(1e-4), em2(1e-2);
71 maxit_ = parlist.sublist("General").sublist("Krylov").get("Iteration Limit",20);
72 tol1_ = parlist.sublist("General").sublist("Krylov").get("Absolute Tolerance",em4);
73 tol2_ = parlist.sublist("General").sublist("Krylov").get("Relative Tolerance",em2);
74 }
75
76 void initialize(const Vector<Real> &x, const Vector<Real> &g) {
77 s_ = x.clone();
78 v_ = x.clone();
79 p_ = x.clone();
80 g_ = g.clone();
81 Hp_ = g.clone();
82 }
83
85 Real &snorm,
86 Real &pRed,
87 int &iflag,
88 int &iter,
89 const Real del,
91 Real tol = std::sqrt(ROL_EPSILON<Real>());
92 const Real zero(0), one(1), two(2), half(0.5);
93 // Initialize step
94 s.zero(); s_->zero();
95 snorm = zero;
96 Real snorm2(0), s1norm2(0);
97 // Compute (projected) gradient
98 g_->set(*model.getGradient());
99 Real gnorm = g_->norm(), normg = gnorm;
100 const Real gtol = std::min(tol1_,tol2_*gnorm);
101 // Preconditioned (projected) gradient vector
102 model.precond(*v_,*g_,s,tol);
103 // Initialize basis vector
104 p_->set(*v_); p_->scale(-one);
105 //Real pnorm2 = v_->dot(g_->dual());
106 Real pnorm2 = v_->apply(*g_);
107 if ( pnorm2 <= zero ) {
108 iflag = 4;
109 iter = 0;
110 return;
111 }
112 // Initialize scalar storage
113 iter = 0; iflag = 0;
114 Real kappa(0), beta(0), sigma(0), alpha(0), tmp(0), sMp(0);
115 Real gv = pnorm2; //v_->dot(g_->dual());
116 pRed = zero;
117 // Iterate CG
118 for (iter = 0; iter < maxit_; iter++) {
119 // Apply Hessian to direction p
120 model.hessVec(*Hp_,*p_,s,tol);
121 // Check positivity of Hessian
122 //kappa = p_->dot(Hp_->dual());
123 kappa = p_->apply(*Hp_);
124 if (kappa <= zero) {
125 sigma = (-sMp+sqrt(sMp*sMp+pnorm2*(del*del-snorm2)))/pnorm2;
126 s.axpy(sigma,*p_);
127 snorm = del;
128 iflag = 2;
129 break;
130 }
131 // Update step
132 alpha = gv/kappa;
133 s_->set(s);
134 s_->axpy(alpha,*p_);
135 s1norm2 = snorm2 + two*alpha*sMp + alpha*alpha*pnorm2;
136 // Check if step exceeds trust region radius
137 if (s1norm2 >= del*del) {
138 sigma = (-sMp+sqrt(sMp*sMp+pnorm2*(del*del-snorm2)))/pnorm2;
139 s.axpy(sigma,*p_);
140 snorm = del;
141 iflag = 3;
142 break;
143 }
144 // Update model predicted reduction
145 pRed += half*alpha*gv;
146 // Set step to temporary step and store norm
147 s.set(*s_);
148 snorm2 = s1norm2;
149 // Check for convergence
150 g_->axpy(alpha,*Hp_);
151 normg = g_->norm();
152 if (normg < gtol) {
153 break;
154 }
155 // Preconditioned updated (projected) gradient vector
156 model.precond(*v_,*g_,s,tol);
157 tmp = gv;
158 //gv = v_->dot(g_->dual());
159 gv = v_->apply(*g_);
160 beta = gv/tmp;
161 // Update basis vector
162 p_->scale(beta);
163 p_->axpy(-one,*v_);
164 sMp = beta*(sMp+alpha*pnorm2);
165 pnorm2 = gv + beta*beta*pnorm2;
166 }
167 // Update model predicted reduction
168 if (iflag > 0) {
169 pRed += sigma*(gv-half*sigma*kappa);
170 }
171 else {
172 snorm = std::sqrt(snorm2);
173 }
174 // Check iteration count
175 if (iter == maxit_) {
176 iflag = 1;
177 }
178 if (iflag != 1) {
179 iter++;
180 }
181 }
182};
183
184} // namespace ROL
185
186#endif
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0_ zero()
Contains definitions of custom data types in ROL.
Provides interface for truncated CG trust-region subproblem solver.
Ptr< Vector< Real > > g_
Ptr< Vector< Real > > p_
Ptr< Vector< Real > > s_
Ptr< Vector< Real > > Hp_
Ptr< Vector< Real > > v_
TruncatedCG_U(ParameterList &parlist)
void initialize(const Vector< Real > &x, const Vector< Real > &g)
void solve(Vector< Real > &s, Real &snorm, Real &pRed, int &iflag, int &iter, const Real del, TrustRegionModel_U< Real > &model)
Provides the interface to evaluate trust-region model functions.
virtual void hessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &s, Real &tol) override
Apply Hessian approximation to vector.
virtual void precond(Vector< Real > &Pv, const Vector< Real > &v, const Vector< Real > &s, Real &tol) override
Apply preconditioner to vector.
virtual const Ptr< const Vector< Real > > getGradient(void) const
Provides interface for and implements trust-region subproblem solvers.
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:84
virtual void set(const Vector &x)
Set where .
Definition: ROL_Vector.hpp:209
virtual void zero()
Set to zero vector.
Definition: ROL_Vector.hpp:167
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
virtual void axpy(const Real alpha, const Vector &x)
Compute where .
Definition: ROL_Vector.hpp:153