MoochoPack : Framework for Large-Scale Optimization Algorithms  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
MoochoPack_ReducedHessianSecantUpdateStd_Step.cpp
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
5 // Copyright (2003) 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 Roscoe A. Bartlett (rabartl@sandia.gov)
38 //
39 // ***********************************************************************
40 // @HEADER
41 
42 #include <ostream>
43 
44 #include "MoochoPack_ReducedHessianSecantUpdateStd_Step.hpp"
45 #include "MoochoPack_moocho_algo_conversion.hpp"
46 #include "MoochoPack_Exceptions.hpp"
47 #include "IterationPack_print_algorithm_step.hpp"
48 #include "AbstractLinAlgPack_VectorSpace.hpp"
49 #include "AbstractLinAlgPack_VectorStdOps.hpp"
50 #include "AbstractLinAlgPack_VectorMutable.hpp"
51 #include "AbstractLinAlgPack_VectorOut.hpp"
52 #include "AbstractLinAlgPack_MatrixSymOpNonsing.hpp"
53 #include "AbstractLinAlgPack_MatrixSymInitDiag.hpp"
54 #include "AbstractLinAlgPack_MatrixOpOut.hpp"
55 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
56 #include "Teuchos_dyn_cast.hpp"
57 
59  const secant_update_ptr_t& secant_update
60  )
61  :secant_update_(secant_update)
62  ,num_basis_(NO_BASIS_UPDATED_YET)
63  ,iter_k_rHL_init_ident_(-1) // not a valid iteration?
64 {}
65 
67  Algorithm& _algo, poss_type step_poss, IterationPack::EDoStepType type
68  ,poss_type assoc_step_poss
69  )
70 {
71  using Teuchos::dyn_cast;
75  using LinAlgOpPack::Vp_V;
76  using LinAlgOpPack::V_VpV;
77  using LinAlgOpPack::V_VmV;
78  using LinAlgOpPack::V_StV;
79  using LinAlgOpPack::V_MtV;
80 
81  NLPAlgo &algo = rsqp_algo(_algo);
82  NLPAlgoState &s = algo.rsqp_state();
83 
84  EJournalOutputLevel olevel = algo.algo_cntr().journal_output_level();
85  EJournalOutputLevel ns_olevel = algo.algo_cntr().null_space_journal_output_level();
86  std::ostream& out = algo.track().journal_out();
87 
88  // print step header.
89  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
90  using IterationPack::print_algorithm_step;
91  print_algorithm_step( algo, step_poss, type, assoc_step_poss, out );
92  }
93 
94  bool return_val = true;
95 
96  // Get iteration quantities
97  IterQuantityAccess<index_type>
98  &num_basis_iq = s.num_basis();
99  IterQuantityAccess<VectorMutable>
100  &pz_iq = s.pz(),
101  &rGf_iq = s.rGf(),
102  &w_iq = s.w();
103  IterQuantityAccess<MatrixOp>
104  &Z_iq = s.Z();
105  IterQuantityAccess<MatrixSymOp>
106  &rHL_iq = s.rHL();
107 
108  // problem size
109  const NLP &nlp = algo.nlp();
110  const size_type
111  //n = nlp.n(),
112  m = nlp.m(),
113  nind = m ? Z_iq.get_k(Z_iq.last_updated()).cols() : 0;
114  //r = m - nind;
115 
116  // See if a new basis has been selected
117  bool new_basis = false;
118  {
119  const int last_updated_k = num_basis_iq.last_updated();
120  if( last_updated_k != IterQuantity::NONE_UPDATED ) {
121  const index_type num_basis_last = num_basis_iq.get_k(last_updated_k);
122  if( num_basis_ == NO_BASIS_UPDATED_YET )
123  num_basis_ = num_basis_last;
124  else if( num_basis_ != num_basis_last )
125  new_basis = true;
126  num_basis_ = num_basis_last;
127  }
128  }
129 
130  // If rHL has already been updated for this iteration then just leave it.
131  if( !rHL_iq.updated_k(0) ) {
132 
133  // If a new basis has been selected, reinitialize
134  if( new_basis ) {
135 
136  if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
137  out << "\nBasis changed. Reinitializing rHL_k = eye(n-r) ...\n";
138  }
139  dyn_cast<MatrixSymInitDiag>(rHL_iq.set_k(0)).init_identity(
140  Z_iq.get_k(0).space_rows()
141  );
142  if( static_cast<int>(ns_olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) )
143  if( algo.algo_cntr().calc_matrix_norms() )
144  out << "\n||rHL_k||inf = " << rHL_iq.get_k(0).calc_norm(MatrixOp::MAT_NORM_INF).value << std::endl;
145  if( (int)ns_olevel >= (int)PRINT_ITERATION_QUANTITIES )
146  out << "\nrHL_k = \n" << rHL_iq.get_k(0);
147  quasi_newton_stats_(s).set_k(0).set_updated_stats(
148  QuasiNewtonStats::REINITIALIZED );
149  iter_k_rHL_init_ident_ = s.k(); // remember what iteration this was
150 
151  }
152  else {
153 
154  // Determine if rHL has been initialized and if we
155  // can perform the update. To perform the BFGS update
156  // rHL_km1 and rGf_km1 must have been computed.
157  if( rHL_iq.updated_k(-1) && rGf_iq.updated_k(-1) ) {
158 
159  // /////////////////////////////////////////////////////
160  // Perform the Secant update
161 
162  if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS )
163  {
164  out
165  << "\nPerforming Secant update ...\n";
166  }
167 
168  const Vector
169  &rGf_k = rGf_iq.get_k(0),
170  &rGf_km1 = rGf_iq.get_k(-1),
171  &pz_km1 = pz_iq.get_k(-1);
172  const value_type
173  alpha_km1 = s.alpha().get_k(-1);
174  VectorSpace::vec_mut_ptr_t
175  y_bfgs = rGf_k.space().create_member(),
176  s_bfgs = pz_km1.space().create_member();
177 
178  // /////////////////////////////////////////////////////
179  // y_bfgs = rGf_k - rGf_km1 - alpha_km1 * w_km1
180 
181  // y_bfgs = rGf_k - rGf_km1
182  V_VmV( y_bfgs.get(), rGf_k, rGf_km1 );
183 
184  if( w_iq.updated_k(-1) )
185  // y_bfgs += - alpha_km1 * w_km1
186  Vp_StV( y_bfgs.get(), - alpha_km1, w_iq.get_k(-1) );
187 
188  // /////////////////////////////////////////////////////
189  // s_bfgs = alpha_km1 * pz_km1
190  V_StV( s_bfgs.get(), alpha_km1, pz_iq.get_k(-1) );
191 
192  if( static_cast<int>(ns_olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
193  out << "\n||y_bfgs||inf = " << y_bfgs->norm_inf() << std::endl;
194  out << "\n||s_bfgs||inf = " << s_bfgs->norm_inf() << std::endl;
195  }
196 
197  if( static_cast<int>(ns_olevel) >= static_cast<int>(PRINT_VECTORS) ) {
198  out << "\ny_bfgs =\n" << *y_bfgs;
199  out << "\ns_bfgs =\n" << *s_bfgs;
200  }
201 
202  // Update from last
203  MatrixSymOp
204  &rHL_k = rHL_iq.set_k(0,-1);
205 
206  // Perform the secant update
207  if(!secant_update().perform_update(
208  s_bfgs.get(), y_bfgs.get()
209  ,iter_k_rHL_init_ident_ == s.k() - 1
210  ,out, ns_olevel, &algo, &s, &rHL_k
211  ))
212  {
213  return_val = false; // redirect control of algorithm!
214  }
215 
216  }
217  else {
218  // We do not have the info to perform the update
219 
220  int k_last_offset = rHL_iq.last_updated();
221  bool set_current = false;
222  if( k_last_offset != IterQuantity::NONE_UPDATED && k_last_offset < 0 ) {
223  const MatrixSymOp &rHL_k_last = rHL_iq.get_k(k_last_offset);
224  const size_type nind_last = rHL_k_last.rows();
225  if( nind_last != nind) {
226  if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
227  out
228  << "No new basis was selected.\n"
229  << "The previous matrix rHL_k(" << k_last_offset << ") was found but its dimmension\n"
230  << "rHL_k(" << k_last_offset << ").rows() = " << nind_last << " != n-r = " << nind << std::endl;
231  }
232  }
233  else {
234  if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
235  out
236  << "No new basis was selected so using previously updated...\n "
237  << "rHL_k = rHL_k(" << k_last_offset << ")\n";
238  }
239  rHL_iq.set_k(0) = rHL_k_last;
240  quasi_newton_stats_(s).set_k(0).set_updated_stats(
241  QuasiNewtonStats::SKIPED );
242 
243  if( (int)ns_olevel >= (int)PRINT_ITERATION_QUANTITIES ) {
244  rHL_iq.get_k(0).output( out << "\nrHL_k = \n" );
245  }
246  set_current = true;
247  }
248  }
249  if( !set_current ) {
250  if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
251  out
252  << "\nInitializing rHL = eye(n-r) "
253  << "(k = " << algo.state().k() << ")...\n";
254  }
255 
256  // Now I will assume that since I can't perform the BFGS update and rHL has
257  // not been set for this iteration yet, that it is up to me to initialize rHL_k = 0
258  dyn_cast<MatrixSymInitDiag>(rHL_iq.set_k(0)).init_identity(
259  Z_iq.get_k(0).space_rows() );
260  iter_k_rHL_init_ident_ = s.k(); // remember what iteration this was
261  quasi_newton_stats_(s).set_k(0).set_updated_stats(
262  QuasiNewtonStats::REINITIALIZED );
263 
264  }
265  }
266 
267  }
268 
269  // Print rHL_k
270 
271  MatrixOp::EMatNormType mat_nrm_inf = MatrixOp::MAT_NORM_INF;
272 
273  if( (int)ns_olevel >= (int)PRINT_ALGORITHM_STEPS ) {
274  if(algo.algo_cntr().calc_matrix_norms())
275  out << "\n||rHL_k||inf = " << rHL_iq.get_k(0).calc_norm(mat_nrm_inf).value << std::endl;
276  if(algo.algo_cntr().calc_conditioning()) {
277  const MatrixSymOpNonsing
278  *rHL_ns_k = dynamic_cast<const MatrixSymOpNonsing*>(&rHL_iq.get_k(0));
279  if(rHL_ns_k)
280  out << "\ncond_inf(rHL_k) = " << rHL_ns_k->calc_cond_num(mat_nrm_inf).value << std::endl;
281  }
282  }
283  if( (int)ns_olevel >= (int)PRINT_ITERATION_QUANTITIES ) {
284  out << "\nrHL_k = \n" << rHL_iq.get_k(0);
285  }
286 
287  }
288  else {
289  if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
290  out << "\nThe matrix rHL_k has already been updated so leave it\n";
291  }
292  }
293 
294  return return_val;
295 }
296 
298  const Algorithm& algo, poss_type step_poss, IterationPack::EDoStepType type, poss_type assoc_step_poss
299  ,std::ostream& out, const std::string& L
300  ) const
301 {
302  out
303  << L << "*** Calculate the reduced hessian of the Lagrangian rHL = Z' * HL * Z\n"
304  << L << "default: num_basis_remembered = NO_BASIS_UPDATED_YET\n"
305  << L << " iter_k_rHL_init_ident = -1\n"
306  << L << "if num_basis_remembered = NO_BASIS_UPDATED_YET then\n"
307  << L << " num_basis_remembered = num_basis\n"
308  << L << "end\n"
309  << L << "if num_basis_remembered != num_basis then\n"
310  << L << " num_basis_remembered = num_basis\n"
311  << L << " new_basis = true\n"
312  << L << "end\n"
313  << L << "if rHL_k is not updated then\n"
314  << L << " if new_basis == true then\n"
315  << L << " *** Transition rHL to the new basis by just starting over.\n"
316  << L << " rHL_k = eye(n-r) *** must support MatrixSymInitDiag interface\n"
317  << L << " iter_k_rHL_init_ident = k\n"
318  << L << " goto next step\n"
319  << L << " end\n"
320  << L << " if rHL_km1 and rGf_km1 are updated then\n"
321  << L << " *** We should have the information to perform a BFGS update\n"
322  << L << " y = rGf_k - rGf_km1\n"
323  << L << " s = alpha_km1 * pz_km1\n"
324  << L << " if k - 1 == iter_k_rHL_init_ident then\n"
325  << L << " first_update = true\n"
326  << L << " else\n"
327  << L << " first_update = false\n"
328  << L << " end\n"
329  << L << " rHL_k = rHL_km1\n"
330  << L << " begin secant update\n"
331  << L << " (" << typeName(secant_update()) << ")\n"
332  ;
333  secant_update().print_step( out, L+" " );
334  out
335  << L << " end secant update\n"
336  << L << " else\n"
337  << L << " *** We have no information for which to preform a BFGS update.\n"
338  << L << " k_last_offset = last iteration rHL was updated for\n"
339  << L << " if k_last_offset does not exist then\n"
340  << L << " *** We are left with no choise but to initialize rHL\n"
341  << L << " rHL_k = eye(n-r) *** must support MatrixSymInitDiag interface\n"
342  << L << " iter_k_rHL_init_ident = k\n"
343  << L << " else\n"
344  << L << " *** No new basis has been selected so we may as well\n"
345  << L << " *** just use the last rHL that was updated\n"
346  << L << " rHL_k = rHL_k(k_last_offset)\n"
347  << L << " end\n"
348  << L << " end\n"
349  << L << "end\n";
350 }
ReducedHessianSecantUpdateStd_Step(const secant_update_ptr_t &secant_update=Teuchos::null)
void Vt_S(VectorMutable *v_lhs, const value_type &alpha)
void Vp_StV(VectorMutable *v_lhs, const value_type &alpha, const Vector &v_rhs)
T_To & dyn_cast(T_From &from)
rSQP Algorithm control class.
void print_step(const Algorithm &algo, poss_type step_poss, IterationPack::EDoStepType type, poss_type assoc_step_poss, std::ostream &out, const std::string &leading_str) const
void V_StV(VectorMutable *v_lhs, value_type alpha, const V &V_rhs)
bool do_step(Algorithm &algo, poss_type step_poss, IterationPack::EDoStepType type, poss_type assoc_step_poss)
virtual std::ostream & journal_out() const
void V_VpV(VectorMutable *v_lhs, const V1 &V1_rhs1, const V2 &V2_rhs2)
Reduced space SQP state encapsulation interface.
size_t size_type
AlgorithmTracker & track()
value_type dot(const Vector &v_rhs1, const Vector &v_rhs2)
AlgorithmState & state()
NLPAlgoState & rsqp_state()
<<std aggr>="">> members for algo_cntr
size_type cols(size_type rows, size_type cols, BLAS_Cpp::Transp _trans)
void V_VmV(VectorMutable *v_lhs, const V1 &V1_rhs1, const V2 &V2_rhs2)
std::string typeName(const T &t)
void Vp_V(VectorMutable *v_lhs, const V &V_rhs)