ROL
ROL_NonlinearCG.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_NONLINEARCG_H
45 #define ROL_NONLINEARCG_H
46 
72 #include "ROL_Types.hpp"
73 
74 namespace ROL {
75 
76 template<class Real>
77 struct NonlinearCGState {
78  std::vector<Teuchos::RCP<Vector<Real> > > grad; // Gradient Storage
79  std::vector<Teuchos::RCP<Vector<Real> > > pstep; // Step Storage
80  int iter; // Nonlinear-CG Iteration Counter
81  int restart; // Reinitialize every 'restart' iterations
82  ENonlinearCG nlcg_type; // Nonlinear-CG Type
83 };
84 
85 template<class Real>
86 class NonlinearCG {
87 private:
88 
89  Teuchos::RCP<NonlinearCGState<Real> > state_; // Nonlinear-CG State
90 
91  Teuchos::RCP<Vector<Real> > y_;
92  Teuchos::RCP<Vector<Real> > yd_;
93 
94 public:
95 
96  virtual ~NonlinearCG() {}
97 
98  // Constructor
99  NonlinearCG(ENonlinearCG type, int restart = 100) {
100  state_ = Teuchos::rcp( new NonlinearCGState<Real> );
101  state_->iter = 0;
102  state_->grad.resize(1);
103  state_->pstep.resize(1);
104  TEUCHOS_TEST_FOR_EXCEPTION(!(isValidNonlinearCG(type)),
105  std::invalid_argument,
106  ">>> ERROR (ROL_NonlinearCG.hpp): Invalid nonlinear CG type in constructor!");
107  state_->nlcg_type = type;
108  TEUCHOS_TEST_FOR_EXCEPTION((restart < 1),
109  std::invalid_argument,
110  ">>> ERROR (ROL_NonlinearCG.hpp): Non-positive restart integer in constructor!");
111  state_->restart = restart;
112  }
113 
114  Teuchos::RCP<NonlinearCGState<Real> >& get_state() { return this->state_; }
115 
116  // Run one step of nonlinear CG.
117  virtual void run( Vector<Real> &s , const Vector<Real> &g, const Vector<Real> &x, Objective<Real> &obj ) {
118  // Initialize vector storage
119  if ( state_->iter == 0 ) {
120  if ( state_->nlcg_type != NONLINEARCG_FLETCHER_REEVES &&
121  state_->nlcg_type != NONLINEARCG_FLETCHER_CONJDESC ) {
122  y_ = g.clone();
123  }
124  if ( state_->nlcg_type == NONLINEARCG_HAGER_ZHANG ||
125  state_->nlcg_type == NONLINEARCG_OREN_LUENBERGER ) {
126  yd_ = g.clone();
127  }
128  }
129 
130  s.set(g.dual());
131 
132  if ((state_->iter % state_->restart) != 0) {
133  Real beta = 0.0;
134  switch(state_->nlcg_type) {
135 
137  y_->set(g);
138  y_->axpy(-1.0, *(state_->grad[0]));
139  beta = - g.dot(*y_) / (state_->pstep[0]->dot(y_->dual()));
140  beta = std::max(beta, 0.0);
141  break;
142  }
143 
145  beta = g.dot(g) / (state_->grad[0])->dot(*(state_->grad[0]));
146  break;
147  }
148 
149  case NONLINEARCG_DANIEL: {
150  Real htol = 0.0;
151  obj.hessVec( *y_, *(state_->pstep[0]), x, htol );
152  beta = - g.dot(*y_) / (state_->pstep[0])->dot(y_->dual());
153  beta = std::max(beta, 0.0);
154  break;
155  }
156 
158  y_->set(g);
159  y_->axpy(-1.0, *(state_->grad[0]));
160  beta = g.dot(*y_) / (state_->grad[0])->dot(*(state_->grad[0]));
161  beta = std::max(beta, 0.0);
162  break;
163  }
164 
166  beta = g.dot(g) / (state_->pstep[0])->dot((state_->grad[0])->dual());
167  break;
168  }
169 
170  case NONLINEARCG_LIU_STOREY: {
171  y_->set(g);
172  y_->axpy(-1.0, *(state_->grad[0]));
173  beta = g.dot(*y_) / (state_->pstep[0])->dot((state_->grad[0])->dual());
174  //beta = std::max(beta, 0.0); // Is this needed? May need research.
175  break;
176  }
177 
178  case NONLINEARCG_DAI_YUAN: {
179  y_->set(g);
180  y_->axpy(-1.0, *(state_->grad[0]));
181  beta = - g.dot(g) / (state_->pstep[0])->dot(y_->dual());
182  break;
183  }
184 
186  Real eta_0 = 1e-2;
187  y_->set(g);
188  y_->axpy(-1.0, *(state_->grad[0]));
189  yd_->set(*y_);
190  Real mult = 2.0 * ( y_->dot(*y_) / (state_->pstep[0])->dot(y_->dual()) );
191  yd_->axpy(-mult, (state_->pstep[0])->dual());
192  beta = - yd_->dot(g) / (state_->pstep[0])->dot(y_->dual());
193  Real eta = -1.0 / ((state_->pstep[0])->norm()*std::min(eta_0,(state_->grad[0])->norm()));
194  beta = std::max(beta, eta);
195  break;
196  }
197 
199  Real eta_0 = 1e-2;
200  y_->set(g);
201  y_->axpy(-1.0, *(state_->grad[0]));
202  yd_->set(*y_);
203  Real mult = ( y_->dot(*y_) / (state_->pstep[0])->dot(y_->dual()) );
204  yd_->axpy(-mult, (state_->pstep[0])->dual());
205  beta = - yd_->dot(g) / (state_->pstep[0])->dot(y_->dual());
206  Real eta = -1.0 / ((state_->pstep[0])->norm()*std::min(eta_0,(state_->grad[0])->norm()));
207  beta = std::max(beta, eta);
208  break;
209  }
210 
211  default:
212  TEUCHOS_TEST_FOR_EXCEPTION(!(isValidNonlinearCG(state_->nlcg_type)),
213  std::invalid_argument,
214  ">>> ERROR (ROL_NonlinearCG.hpp): Invalid nonlinear CG type in the 'run' method!");
215  }
216 
217  s.axpy(beta, *(state_->pstep[0]));
218  }
219 
220  // Update storage.
221  if (state_->iter == 0) {
222  (state_->grad[0]) = g.clone();
223  (state_->pstep[0]) = s.clone();
224  }
225  (state_->grad[0])->set(g);
226  (state_->pstep[0])->set(s);
227  state_->iter++;
228  }
229 
230 };
231 
232 }
233 
234 #endif
Contains definitions of custom data types in ROL.
ENonlinearCG
Enumeration of nonlinear CG algorithms.
Definition: ROL_Types.hpp:536
int isValidNonlinearCG(ENonlinearCG s)
Verifies validity of a NonlinearCG enum.
Definition: ROL_Types.hpp:572