MoochoPack : Framework for Large-Scale Optimization Algorithms  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
MoochoPack_ReducedHessianSecantUpdateBFGSProjected_Strategy.cpp
1 #if 0
2 
3 // @HEADER
4 // ***********************************************************************
5 //
6 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
7 // Copyright (2003) Sandia Corporation
8 //
9 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
10 // license for use of this work by or on behalf of the U.S. Government.
11 //
12 // Redistribution and use in source and binary forms, with or without
13 // modification, are permitted provided that the following conditions are
14 // met:
15 //
16 // 1. Redistributions of source code must retain the above copyright
17 // notice, this list of conditions and the following disclaimer.
18 //
19 // 2. Redistributions in binary form must reproduce the above copyright
20 // notice, this list of conditions and the following disclaimer in the
21 // documentation and/or other materials provided with the distribution.
22 //
23 // 3. Neither the name of the Corporation nor the names of the
24 // contributors may be used to endorse or promote products derived from
25 // this software without specific prior written permission.
26 //
27 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
28 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
29 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
30 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
31 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
32 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
33 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
34 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
35 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
36 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
37 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
38 //
39 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov)
40 //
41 // ***********************************************************************
42 // @HEADER
43 
44 #include "MoochoPack_ReducedHessianSecantUpdateBFGSProjected_Strategy.hpp"
45 #include "MoochoPack_PBFGS_helpers.hpp"
46 #include "MoochoPack_NLPAlgo.hpp"
47 #include "MoochoPack_NLPAlgoState.hpp"
48 #include "ConstrainedOptPack/src/AbstractLinAlgPack_MatrixSymAddDelUpdateable.hpp"
49 #include "ConstrainedOptPack/src/AbstractLinAlgPack_BFGS_helpers.hpp"
50 #include "AbstractLinAlgPack/src/AbstractLinAlgPack_SpVectorClass.hpp"
51 #include "AbstractLinAlgPack_SpVectorOp.hpp"
52 #include "AbstractLinAlgPack/src/AbstractLinAlgPack_MatrixOpOut.hpp"
53 #include "AbstractLinAlgPack/src/AbstractLinAlgPack_GenPermMatrixSlice.hpp"
54 #include "AbstractLinAlgPack_GenPermMatrixSliceOp.hpp"
55 #include "AbstractLinAlgPack/src/AbstractLinAlgPack_GenPermMatrixSliceOut.hpp"
56 #include "AbstractLinAlgPack/src/AbstractLinAlgPack_MatrixSymInitDiag.hpp"
57 #include "DenseLinAlgPack_LinAlgOpPack.hpp"
58 #include "Midynamic_cast_verbose.h"
59 #include "MiWorkspacePack.h"
60 
61 namespace LinAlgOpPack {
63 }
64 
65 namespace MoochoPack {
66 
68  const bfgs_update_ptr_t& bfgs_update
69  ,value_type act_set_frac_proj_start
70  ,value_type project_error_tol
71  ,value_type super_basic_mult_drop_tol
72  )
73  : bfgs_update_(bfgs_update)
74  , act_set_frac_proj_start_(act_set_frac_proj_start)
75  , project_error_tol_(project_error_tol)
76  , super_basic_mult_drop_tol_(super_basic_mult_drop_tol)
77 {}
78 
80  DVectorSlice* s_bfgs, DVectorSlice* y_bfgs, bool first_update
81  ,std::ostream& out, EJournalOutputLevel olevel, NLPAlgo *algo, NLPAlgoState *s
82  ,MatrixOp *rHL_k
83  )
84 {
85  using std::setw;
86  using std::endl;
87  using std::right;
88  using Teuchos::dyn_cast;
90  using LinAlgOpPack::V_MtV;
91  using AbstractLinAlgPack::norm_inf;
93  using Teuchos::Workspace;
95 
96  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
97  out << "\n*** (PBFGS) Projected BFGS ...\n";
98  }
99 
100 #ifdef _WINDOWS
101  MHSB_t &rHL_super = dynamic_cast<MHSB_t&>(*rHL_k);
102 #else
103  MHSB_t &rHL_super = dyn_cast<MHSB_t>(*rHL_k);
104 #endif
105 
106  const size_type
107  n = algo->nlp().n(),
108  r = algo->nlp().r(),
109  n_pz = n-r;
110  const GenPermMatrixSlice
111  &Q_R = rHL_super.Q_R(),
112  &Q_X = rHL_super.Q_X();
113 
114  bool do_projected_rHL_RR = false;
115 
116  // Check that the current update is sufficiently p.d. before we do anything
117  const value_type
118  sTy = dot(*s_bfgs,*y_bfgs),
119  yTy = dot(*y_bfgs,*y_bfgs);
120  if( !ConstrainedOptPack::BFGS_sTy_suff_p_d(
121  *s_bfgs,*y_bfgs,&sTy
122  , int(olevel) >= int(PRINT_ALGORITHM_STEPS) ? &out : NULL )
123  && !bfgs_update().use_dampening()
124  )
125  {
126  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
127  out << "\nWarning! use_damening == false so there is no way we can perform any kind BFGS update (projected or not) so we will skip it!\n";
128  }
129  quasi_newton_stats_(*s).set_k(0).set_updated_stats(
130  QuasiNewtonStats::INDEF_SKIPED );
131  return true;
132  }
133 
134  // Get matrix scaling
135  value_type
136  rHL_XX_scale = sTy > 0.0 ? yTy/sTy : 1.0;
137 
138  //
139  // Initialize or adjust the active set before the BFGS update
140  //
141  if( !s->nu().updated_k(-1) ) {
142  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
143  out << "\nWarning! nu_k(-1) has not been updated. No adjustment to the set of superbasic variables is possible!\n";
144  }
145  }
146  else if( Q_R.is_identity() ) {
147  // Determine when to start adding and removing rows/cols form rHL_RR
148  if( act_set_stats_(*s).updated_k(-1) ) {
149  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
150  out << "\nDetermining if projected BFGS updating of superbasics should begin ...\n";
151  }
152  // Determine if the active set has calmed down enough
153  const SpVector
154  &nu_km1 = s->nu().get_k(-1);
155  const SpVectorSlice
156  nu_indep = nu_km1(s->var_indep());
157  do_projected_rHL_RR = PBFGSPack::act_set_calmed_down(
158  act_set_stats_(*s).get_k(-1)
159  ,act_set_frac_proj_start()
160  ,olevel,out
161  );
162  if( do_projected_rHL_RR ) {
163  //
164  // Determine the set of initially fixed and free independent variables.
165  //
166  typedef Workspace<size_type> i_x_t;
167  typedef Workspace<ConstrainedOptPack::EBounds> bnd_t;
168  i_x_t i_x_free(wss,n_pz);
169  i_x_t i_x_fixed(wss,n_pz);
170  bnd_t bnd_fixed(wss,n_pz);
171  i_x_t l_x_fixed_sorted(wss,n_pz);
172  size_type n_pz_X = 0, n_pz_R = 0;
173  value_type sRTBRRsR = 0.0, sRTyR = 0.0, sXTBXXsX = 0.0, sXTyX = 0.0;
174  // Get the elements in i_x_free[] for variables that are definitely free
175  // and initialize s_R'*y_R
176  PBFGSPack::init_i_x_free_sRTsR_sRTyR(
177  nu_indep, *s_bfgs, *y_bfgs
178  , &n_pz_R, &i_x_free[0], &sRTBRRsR, &sRTyR ); // We don't really want sRTRBBsR here
179  // rHL_XX = some_scaling * I
180  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
181  out << "\nScaling for diagonal rHL_XX = rHL_XX_scale*I, rHL_XX_scale = " << rHL_XX_scale << std::endl;
182  }
183  Workspace<value_type> rHL_XX_diag_ws(wss,nu_indep.nz());
184  DVectorSlice rHL_XX_diag(&rHL_XX_diag_ws[0],rHL_XX_diag_ws.size());
185  rHL_XX_diag = rHL_XX_scale;
186  // s_R'*B_RR_*s_R
187  Workspace<value_type> Q_R_Q_RT_s_ws(wss,n_pz);
188  DVectorSlice Q_R_Q_RT_s(&Q_R_Q_RT_s_ws[0],Q_R_Q_RT_s_ws.size());
189  Q_R_Q_RT_s = 0.0;
190  {for( size_type k = 0; k < n_pz_R; ++k ) {
191  const size_type i = i_x_free[k];
192  Q_R_Q_RT_s(i) = (*s_bfgs)(i);
193  }}
194  sRTBRRsR = AbstractLinAlgPack::transVtMtV( Q_R_Q_RT_s, *rHL_k, BLAS_Cpp::no_trans, Q_R_Q_RT_s );
195  // Sort fixed variables according to |s_X(i)^2*B_XX(i,i)|/|sRTBRRsR| + |s_X(i)*y_X(i)|/|sRTyR|
196  // and initialize s_X'*B_XX*s_X and s_X*y_X
197  PBFGSPack::sort_fixed_max_cond_viol(
198  nu_indep,*s_bfgs,*y_bfgs,rHL_XX_diag,sRTBRRsR,sRTyR
199  ,&sXTBXXsX,&sXTyX,&l_x_fixed_sorted[0]);
200  // Pick initial set of i_x_free[] and i_x_fixed[] (sorted!)
201  PBFGSPack::choose_fixed_free(
202  project_error_tol(),super_basic_mult_drop_tol(),nu_indep
203  ,*s_bfgs,*y_bfgs,rHL_XX_diag,&l_x_fixed_sorted[0]
204  ,olevel,out,&sRTBRRsR,&sRTyR,&sXTBXXsX,&sXTyX
205  ,&n_pz_X,&n_pz_R,&i_x_free[0],&i_x_fixed[0],&bnd_fixed[0] );
206  //
207  // Delete rows/cols from rHL_RR for fixed variables
208  //
209 #ifdef _WINDOWS
210  MatrixSymAddDelUpdateable
211  &rHL_RR = dynamic_cast<MatrixSymAddDelUpdateable&>(
212  const_cast<MatrixSymWithOpFactorized&>(*rHL_super.B_RR_ptr())
213  );
214 #else
215  MatrixSymAddDelUpdateable
216  &rHL_RR = dyn_cast<MatrixSymAddDelUpdateable>(
217  const_cast<MatrixSymWithOpFactorized&>(*rHL_super.B_RR_ptr())
218  );
219 #endif
220  if( n_pz_R < n_pz ) {
221  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
222  out << "\nDeleting n_pz_X = " << n_pz_X << " rows/columns from rHL_RR for fixed independent variables...\n";
223  }
224  {for( size_type k = n_pz_X; k > 0; --k ) { // Delete from the largest to the smallest index (cheaper!)
225  rHL_RR.delete_update( i_x_fixed[k-1], false );
226  }}
227  TEUCHOS_TEST_FOR_EXCEPT( !( rHL_RR.rows() == n_pz_R ) );
228  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ITERATION_QUANTITIES) ) {
229  out << "\nrHL_RR after rows/columns where removed =\n" << *rHL_super.B_RR_ptr();
230  }
231  // Initialize rHL_XX = rHL_XX_scale*I
232 #ifdef _WINDOWS
233  MatrixSymInitDiag
234  &rHL_XX = dynamic_cast<MatrixSymInitDiag&>(
235  const_cast<MatrixSymOp&>(*rHL_super.B_XX_ptr())
236  );
237 #else
238  MatrixSymInitDiag
239  &rHL_XX = dyn_cast<MatrixSymInitDiag>(
240  const_cast<MatrixSymOp&>(*rHL_super.B_XX_ptr())
241  );
242 #endif
243  rHL_XX.init_identity(n_pz_X,rHL_XX_scale);
244  // Reinitialize rHL for new active set
245  rHL_super.initialize(
246  n_pz, n_pz_R, &i_x_free[0], &i_x_fixed[0], &bnd_fixed[0]
247  ,rHL_super.B_RR_ptr(),NULL,BLAS_Cpp::no_trans,rHL_super.B_XX_ptr()
248  );
249  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ITERATION_QUANTITIES) ) {
250  out << "\nFull rHL after reinitialization but before BFGS update:\n"
251  << "\nrHL =\n" << *rHL_k
252  << "\nQ_R =\n" << rHL_super.Q_R()
253  << "\nQ_X =\n" << rHL_super.Q_X();
254  }
255  }
256  else {
257  do_projected_rHL_RR = false;
258  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
259  out << "\nWith n_pz_X = " << n_pz_X << ", there where no variables to drop from superbasis!\n";
260  }
261  }
262  }
263  }
264  }
265  else {
266  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
267  out << "\nAdjust the set of superbasic variables and the projected reduced Hessian rHL_RR ...\n";
268  }
269  //
270  // Modify rHL_RR by adding and dropping rows/cols for freeded and fixed variables
271  //
272  const SpVectorSlice
273  nu_indep = s->nu().get_k(-1)(s->var_indep());
274  //
275  // Determine new Q_R and Q_X
276  //
277  typedef Workspace<size_type> i_x_t;
278  typedef Workspace<ConstrainedOptPack::EBounds> bnd_t;
279  i_x_t i_x_free(wss,n_pz);
280  i_x_t i_x_fixed(wss,n_pz);
281  bnd_t bnd_fixed(wss,n_pz);
282  i_x_t l_x_fixed_sorted(wss,n_pz);
283  size_type n_pz_X = 0, n_pz_R = 0;
284  value_type sRTBRRsR = 0.0, sRTyR = 0.0, sXTBXXsX = 0.0, sXTyX = 0.0;
285  // Get the elements in i_x_free[] for variables that are definitely free
286  // and initialize s_R'*y_R. This will be the starting point for the new Q_R.
287  PBFGSPack::init_i_x_free_sRTsR_sRTyR(
288  nu_indep, *s_bfgs, *y_bfgs
289  , &n_pz_R, &i_x_free[0], &sRTBRRsR, &sRTyR ); // We don't really want sRTBRRsR here
290  // Initialize rHL_XX_diag = some_scaling * I as though all of the currently fixed variables
291  // will be left out of Q_R. Some of these variables might already be in Q_R and B_RR
292  // and may still be in Q_R and B_RR after we are finished adjusting the sets Q_R and Q_X
293  // and we don't want to delete these rows/cols in B_RR just yet!
294  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
295  out << "\nScaling for diagonal rHL_XX = rHL_XX_scale*I, rHL_XX_scale = " << rHL_XX_scale << std::endl;
296  }
297  Workspace<value_type> rHL_XX_diag_ws(wss,nu_indep.nz());
298  DVectorSlice rHL_XX_diag(&rHL_XX_diag_ws[0],rHL_XX_diag_ws.size());
299  rHL_XX_diag = rHL_XX_scale;
300  // Initialize rHL_XX = rHL_XX_scale * I so that those variables in the current Q_R
301  // not in the estimate i_x_free[] will have their proper value when s_R'*B_RR*s_R computed
302  // for the estimate i_x_free[]. This is needed to change the behavior of *rHL_k which
303  // is used below to compute s_R'*B_RR*s_R
304 #ifdef _WINDOWS
305  MatrixSymInitDiag
306  &rHL_XX = dynamic_cast<MatrixSymInitDiag&>(
307  const_cast<MatrixSymOp&>(*rHL_super.B_XX_ptr())
308  );
309 #else
310  MatrixSymInitDiag
311  &rHL_XX = dyn_cast<MatrixSymInitDiag>(
312  const_cast<MatrixSymOp&>(*rHL_super.B_XX_ptr())
313  );
314 #endif
315  rHL_XX.init_identity(rHL_XX.rows(),rHL_XX_scale); // Don't change its size yet!
316  // s_R'*B_RR_*s_R
317  // This will only include those terms for the variable actually free.
318  Workspace<value_type> Q_R_Q_RT_s_ws(wss,n_pz);
319  DVectorSlice Q_R_Q_RT_s(&Q_R_Q_RT_s_ws[0],Q_R_Q_RT_s_ws.size());
320  Q_R_Q_RT_s = 0.0;
321  {for( size_type k = 0; k < n_pz_R; ++k ) {
322  const size_type i = i_x_free[k];
323  Q_R_Q_RT_s(i) = (*s_bfgs)(i);
324  }}
325  sRTBRRsR = AbstractLinAlgPack::transVtMtV( Q_R_Q_RT_s, *rHL_k, BLAS_Cpp::no_trans, Q_R_Q_RT_s );
326  // Sort fixed variables according to |s_X(i)^2*B_XX(i,i)|/|sRTBRRsR| + |s_X(i)*y_X(i)|/|sRTyR|
327  PBFGSPack::sort_fixed_max_cond_viol(
328  nu_indep,*s_bfgs,*y_bfgs,rHL_XX_diag,sRTBRRsR,sRTyR
329  ,&sXTBXXsX,&sXTyX,&l_x_fixed_sorted[0]);
330  // Pick initial set of i_x_free[] and i_x_fixed[] (sorted!)
331  PBFGSPack::choose_fixed_free(
332  project_error_tol(),super_basic_mult_drop_tol(),nu_indep
333  ,*s_bfgs,*y_bfgs,rHL_XX_diag,&l_x_fixed_sorted[0]
334  ,olevel,out,&sRTBRRsR,&sRTyR,&sXTBXXsX,&sXTyX
335  ,&n_pz_X,&n_pz_R,&i_x_free[0],&i_x_fixed[0],&bnd_fixed[0] );
336  // Get the changes to the set of superbasic variables
337  size_type num_free_to_fixed = 0, num_fixed_to_free = 0;
338  i_x_t i_x_free_to_fixed(wss,Q_R.cols());
339  i_x_t i_x_fixed_to_free(wss,Q_X.cols());
340  i_x_t i_x_free_still(wss,Q_R.cols()); // Will be set to those indices still in Q_R
341  std::fill_n( &i_x_free_still[0], Q_R.cols(), 0 ); // in the same order as in Q_R and B_RR
342  {
343  GenPermMatrixSlice::const_iterator
344  Q_R_begin = Q_R.begin(),
345  Q_R_itr = Q_R_begin,
346  Q_R_end = Q_R.end();
347  const size_type
348  *i_x_free_itr = &i_x_free[0],
349  *i_x_free_end = i_x_free_itr + n_pz_R;
350  for( size_type i = 1; i <= n_pz; ++i ) {
351  if( Q_R_itr == Q_R_end && i_x_free_itr == i_x_free_end ) {
352  break; // The rest of these variables were and still are not superbasic
353  }
354  else if( i_x_free_itr == i_x_free_end ) {
355  // A variable that was in the superbasis now is not
356  i_x_free_to_fixed[num_free_to_fixed] = Q_R_itr->row_i();
357  num_free_to_fixed++;
358  ++Q_R_itr;
359  }
360  else if( Q_R_itr == Q_R_end ) {
361  // A variable that was not in the superbasis now is
362  i_x_fixed_to_free[num_fixed_to_free] = *i_x_free_itr;
363  num_fixed_to_free++;
364  ++i_x_free_itr;
365  }
366  else {
367  if( Q_R_itr->row_i() == *i_x_free_itr ) {
368  // Both still superbasic
369  i_x_free_still[Q_R_itr-Q_R_begin] = Q_R_itr->row_i();
370  ++Q_R_itr;
371  ++i_x_free_itr;
372  }
373  else if( Q_R_itr->row_i() < *i_x_free_itr ) {
374  // A variable that was in the superbasis now is not
375  i_x_free_to_fixed[num_free_to_fixed] = Q_R_itr->row_i();
376  num_free_to_fixed++;
377  ++Q_R_itr;
378  }
379  else {
380  // A variable that was not in the superbasis now is
381  i_x_fixed_to_free[num_fixed_to_free] = *i_x_free_itr;
382  num_fixed_to_free++;
383  ++i_x_free_itr;
384  }
385  }
386  }
387  }
388  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
389  out << "\nThere will be " << num_fixed_to_free << " independent variables added to the superbasis and rHL_RR";
390  if( num_fixed_to_free && int(olevel) >= int(PRINT_ACTIVE_SET) ) {
391  out << " and their indexes are:\n";
392  for(size_type k = 0; k < num_fixed_to_free; ++k)
393  out << " " << i_x_fixed_to_free[k];
394  out << std::endl;
395  }
396  else {
397  out << "\n";
398  }
399  out << "\nThere will be " << num_free_to_fixed << " independent variables removed from the superbasis and rHL_RR";
400  if( num_free_to_fixed && int(olevel) >= int(PRINT_ACTIVE_SET) ) {
401  out << " and their indexes are:\n";
402  for(size_type k = 0; k < num_free_to_fixed; ++k)
403  out << " " << i_x_free_to_fixed[k];
404  out << std::endl;
405  }
406  else {
407  out << "\n";
408  }
409  }
410  // Get reference to rHL_RR = B_RR
411 #ifdef _WINDOWS
412  MatrixSymAddDelUpdateable
413  &rHL_RR = dynamic_cast<MatrixSymAddDelUpdateable&>(
414  const_cast<MatrixSymWithOpFactorized&>(*rHL_super.B_RR_ptr())
415  );
416 #else
417  MatrixSymAddDelUpdateable
418  &rHL_RR = dyn_cast<MatrixSymAddDelUpdateable>(
419  const_cast<MatrixSymWithOpFactorized&>(*rHL_super.B_RR_ptr())
420  );
421 #endif
422  // Remove rows/cols from rHL_RR.
423  if( num_free_to_fixed ) {
424  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
425  out << "\nDeleting " << num_free_to_fixed << " rows/columns from rHL_RR ...\n";
426  }
427  {for( size_type k = i_x_free_still.size(); k > 0; --k ) { // Delete from the largest to the smallest index (cheaper!)
428  if( !i_x_free_still[k-1] )
429  rHL_RR.delete_update( k, false );
430  }}
431  TEUCHOS_TEST_FOR_EXCEPT( !( rHL_RR.rows() == n_pz_R - num_fixed_to_free ) );
432  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ITERATION_QUANTITIES) ) {
433  out << "\nrHL_RR after rows/columns where removed =\n" << *rHL_super.B_RR_ptr();
434  }
435  }
436  // Add new rows/cols to rHL_RR.
437  if( num_fixed_to_free ) {
438  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
439  out << "\nAppending " << num_fixed_to_free << " rows/columns to rHL_RR ...\n";
440  }
441  {for( size_type k = 0; k < num_fixed_to_free; ++k ) {
442  rHL_RR.augment_update( NULL, rHL_XX_scale, false );
443  }}
444  TEUCHOS_TEST_FOR_EXCEPT( !( rHL_RR.rows() == n_pz_R ) );
445  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ITERATION_QUANTITIES) ) {
446  out << "\nrHL_RR after rows/columns where appended =\n" << *rHL_super.B_RR_ptr();
447  }
448  }
449  // Resort i_x_free[] to reflect the actual order of the indices in rHL_RR
450  {
451  size_type tmp_n_pz_R = 0;
452  {for(size_type k = 0; k < i_x_free_still.size(); ++k) {
453  if( i_x_free_still[k] ) {
454  i_x_free[tmp_n_pz_R] = i_x_free_still[k];
455  ++tmp_n_pz_R;
456  }
457  }}
458  {for(size_type k = 0; k < num_fixed_to_free; ++k) {
459  i_x_free[tmp_n_pz_R] = i_x_fixed_to_free[k];
460  ++tmp_n_pz_R;
461  }}
462  TEUCHOS_TEST_FOR_EXCEPT( !( tmp_n_pz_R == n_pz_R ) );
463  }
464  // Initialize rHL_XX = rHL_XX_scale * I resized to the proper dimensions
465  rHL_XX.init_identity(n_pz_X,rHL_XX_scale);
466  // Reinitalize rHL for new active set
467  rHL_super.initialize(
468  n_pz, n_pz_R, &i_x_free[0], &i_x_fixed[0], &bnd_fixed[0]
469  ,rHL_super.B_RR_ptr(),NULL,BLAS_Cpp::no_trans,rHL_super.B_XX_ptr()
470  );
471  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ITERATION_QUANTITIES) ) {
472  out << "\nFull rHL after reinitialization but before BFGS update:\n"
473  << "\nrHL =\n" << *rHL_k
474  << "\nQ_R =\n" << rHL_super.Q_R()
475  << "\nQ_X =\n" << rHL_super.Q_X();
476  }
477  // Now we will do the PBFGS updating from now on!
478  do_projected_rHL_RR = true;
479  }
480  //
481  // Perform the BFGS update
482  //
483  if( do_projected_rHL_RR ) {
484  // Perform BFGS update on smaller rHL_RR.
485  // By the time we get here rHL_RR should be resize and ready to update
486  const GenPermMatrixSlice
487  &Q_R = rHL_super.Q_R(),
488  &Q_X = rHL_super.Q_X();
489  const size_type
490  n_pz_R = Q_R.cols(),
491  n_pz_X = Q_X.cols();
492  TEUCHOS_TEST_FOR_EXCEPT( !( n_pz_R + n_pz_X == n_pz ) );
493  // Get projected BFGS update vectors y_bfgs_R, s_bfgs_R
494  Workspace<value_type>
495  y_bfgs_R_ws(wss,Q_R.cols()),
496  s_bfgs_R_ws(wss,Q_R.cols());
497  DVectorSlice y_bfgs_R(&y_bfgs_R_ws[0],y_bfgs_R_ws.size());
498  DVectorSlice s_bfgs_R(&s_bfgs_R_ws[0],s_bfgs_R_ws.size());
499  V_MtV( &y_bfgs_R, Q_R, BLAS_Cpp::trans, *y_bfgs ); // y_bfgs_R = Q_R'*y_bfgs
500  V_MtV( &s_bfgs_R, Q_R, BLAS_Cpp::trans, *s_bfgs ); // s_bfgs_R = Q_R'*s_bfgs
501  // Update rHL_RR
502  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
503  out << "\nPerform BFGS update on " << n_pz_R << " x " << n_pz_R << " projected reduced Hessian for the superbasic variables where B = rHL_RR...\n";
504  }
505  bfgs_update().perform_update(
506  &s_bfgs_R(),&y_bfgs_R(),first_update,out,olevel,algo->algo_cntr().check_results()
507  ,const_cast<MatrixOp*>(static_cast<const MatrixOp*>(rHL_super.B_RR_ptr().get()))
508  ,&quasi_newton_stats_(*s).set_k(0)
509  );
510  }
511  else {
512  // Update the full reduced Hessain matrix (rHL = rHL_RR)
513  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
514  out << "\nPerform BFGS update on the full reduced Hessian where B = rHL...\n";
515  }
516  bfgs_update().perform_update(
517  s_bfgs,y_bfgs,first_update,out,olevel,algo->algo_cntr().check_results()
518  ,const_cast<MatrixOp*>(static_cast<const MatrixOp*>(rHL_super.B_RR_ptr().get()))
519  ,&quasi_newton_stats_(*s).set_k(0)
520  );
521  }
522 
523  return true;
524 }
525 
526 void ReducedHessianSecantUpdateBFGSProjected_Strategy::print_step( std::ostream& out, const std::string& L ) const
527 {
528  out
529  << L << "*** Perform BFGS update on only free independent (superbasic) variables.\n"
530  << L << "if s'*y < sqrt(macheps) * ||s||2 * ||y||2 and use_dampening == false then"
531  << L << " Skip the update and exit this step!\n"
532  << L << "end\n"
533  << L << "rHL_XX_scale = max(y'*y/s'*y,1.0)\n"
534  << L << "do_projected_rHL_RR = false\n"
535  << L << "if nu_km1 is updated then\n"
536  << L << " if rHL_k.Q_R is the identity matrix then\n"
537  << L << " *** Determine if the active set has calmed down enough\n"
538  << L << " Given (num_active_indep,num_adds_indep,num_drops_indep) from act_set_stats_km1\n"
539  << L << " fact_same =\n"
540  << L << " ( num_adds_indep== NOT_KNOWN || num_drops_indep==NOT_KNOWN\n"
541  << L << " || num_active_indep==0\n"
542  << L << " ? 0.0\n"
543  << L << " : std::_MAX((num_active_indep-num_adds_indep-num_drops_indep)\n"
544  << L << " /num_active_indep,0.0)\n"
545  << L << " )\n"
546  << L << " do_projected_rHL_RR\n"
547  << L << " = fact_same >= act_set_frac_proj_start && num_active_indep > 0\n"
548  << L << " if do_projected_rHL_RR == true then\n"
549  << L << " Determine the sets of superbasic variables given the mapping matrix\n"
550  << L << " Q = [ Q_R, Q_X ] where pz_R = Q_R'*pz <: R^n_pz_R are the superbasic variables and\n"
551  << L << " pz_X = Q_X'*pz <: R^n_pz_X are the nonbasic variables that only contain fixed\n"
552  << L << " variables in nu_km1(indep) where the following condidtions are satisfied:\n"
553  << L << " (s'*Q_X*rHL_XX_scale*Q_X'*s)/(s'*Q_R*Q_R'*rHL_k*Q_R*Q_R'*s) <= project_error_tol\n"
554  << L << " |s'*Q_X*Q_X'*y|/|s'*Q_R*Q_R'*s| <= project_error_tol\n"
555  << L << " |Q_X'*nu_km1(indep)|/||nu_km1(indep)||inf >= super_basic_mult_drop_tol\n"
556  << L << " if n_pz_R < n-r then\n"
557  << L << " Delete rows/columns of rHL_k to form rHL_RR = Q_R'*rHL_k*Q_R\n"
558  << L << " Define new rHL_k = [ Q_R, Q_X ] * [ rHL_RR, 0; 0; rHL_XX_scale*I ] [ Q_R'; Q_X ]\n"
559  << L << " else\n"
560  << L << " do_projected_rHL_RR = false\n"
561  << L << " end\n"
562  << L << " end\n"
563  << L << " else\n"
564  << L << " Determine the new Q_n = [ Q_R_n, Q_X_n ] that satisfies:\n"
565  << L << " (s'*Q_X_n*rHL_XX_scale*Q_X_n'*s)/(s'*Q_R_n*Q_R_n'*rHL_k*Q_R_n*Q_R_n'*s) <= project_error_tol\n"
566  << L << " |s'*Q_X_n*Q_X_n'*y|/|s'*Q_R_n*Q_R_n'*s| <= project_error_tol\n"
567  << L << " |Q_X_n'*nu_km1(indep)|/||nu_km1(indep)||inf >= super_basic_mult_drop_tol\n"
568  << L << " Remove rows/cols from rHL_k.rHL_RR for variables in rHL_k.Q_R that are not in Q_R_n.\n"
569  << L << " Add digonal entries equal to rHL_XX_scale to rHL_k.rHL_RR for variables in Q_R_n\n"
570  << L << " that are not in rHL_k.Q_R\n"
571  << L << " Define new rHL_k = [ Q_R_n, Q_X_n ] * [ rHL_k.rHL_RR, 0; 0; rHL_XX_scale*I ] [ Q_R_n'; Q_X_n ]\n"
572  << L << " do_projected_rHL_RR = true\n"
573  << L << " end\n"
574  << L << "end\n"
575  << L << "if do_projected_rHL_RR == true then\n"
576  << L << " Perform projected BFGS update (see below): (rHL_k.rHL_RR, Q_R_n'*s, Q_R_n'*y) -> rHL_k.rHL_RR\n"
577  << L << "else\n"
578  << L << " Perform full BFGS update: (rHL_k, s, y) -> rHL_k\n"
579  << L << " begin BFGS update where B = rHL_k\n";
580  bfgs_update().print_step( out, L + " " );
581  out
582  << L << " end BFGS update\n"
583  << L << "else\n"
584  ;
585 }
586 
587 } // end namespace MoochoPack
588 
589 #endif // 0
T_To & dyn_cast(T_From &from)
value_type transVtMtV(const Vector &v_rhs1, const MatrixOp &M_rhs2, BLAS_Cpp::Transp trans_rhs2, const Vector &v_rhs3)
bool perform_update(DVectorSlice *s_bfgs, DVectorSlice *y_bfgs, bool first_update, std::ostream &out, EJournalOutputLevel olevel, NLPAlgo *algo, NLPAlgoState *s, MatrixOp *rHL_k)
void print_step(std::ostream &out, const std::string &leading_str) const
size_t size_type
void Vp_StMtV(VectorMutable *v_lhs, value_type alpha, const MatrixOp &M_rhs1, BLAS_Cpp::Transp trans_rhs1, const Vector &v_rhs2, value_type beta=1.0)
value_type dot(const Vector &v_rhs1, const Vector &v_rhs2)
ReducedHessianSecantUpdateBFGSProjected_Strategy(const bfgs_update_ptr_t &bfgs_update=NULL, value_type act_set_frac_proj_start=0.8, value_type project_error_tol=1e-5, value_type super_basic_mult_drop_tol=1e-5)
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
TEUCHOSCORE_LIB_DLL_EXPORT Teuchos::RCP< WorkspaceStore > get_default_workspace_store()