ConstrainedOptPack: C++ Tools for Constrained (and Unconstrained) Optimization  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
ConstrainedOptPack_MatrixSymPosDefLBFGS.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 // Let's define a compact representation for the matrix B^{k} and
43 // its inverse H^{k} = inv(B^{k}).
44 //
45 // Bk = (1/gk)*I - [ (1/gk)*S Y ] * inv( [ (1/gk)*S'S L ] [ (1/gk)*S' ]
46 // [ L' -D ] ) * [ Y' ]
47 // \________________/
48 // Q
49 //
50 // Hk = gk*I + [ S gk*Y ] * [ inv(R')*(D+gk*Y'Y)*inv(R) -inv(R') ] * [ S' ]
51 // [ -inv(R) 0 ] [ gk*Y' ]
52 //
53 // where:
54 //
55 // gk = gamma_k <: R
56 //
57 // [ s^{1}'*y^{1} s^{1}'*y^{2} ... s^{1}'*y^{m} ]
58 // S'Y = [ s^{2}'*y^{1} s^{2}'*y^{2} ... s^{2}'*y^{m} ] <: R^(m x m)
59 // [ . . . ]
60 // [ s^{m}'*y^{1} s^{m}'*y^{2} ... s^{m}'*y^{m} ]
61 //
62 // [ s^{1}'*y^{1} 0 ... 0 ]
63 // D = [ 0 s^{2}'*y^{2} ... 0 ] <: R^(m x m)
64 // [ . . . ]
65 // [ 0 0 ... s^{m}'*y^{m} ]
66 //
67 // R = upper triangular part of S'Y
68 //
69 // L = lower tirangular part of S'Y with zeros on the diagonal
70 //
71 
72 #include <assert.h>
73 
74 #include <typeinfo>
75 
76 #include "ConstrainedOptPack_MatrixSymPosDefLBFGS.hpp"
77 #include "AbstractLinAlgPack_BFGS_helpers.hpp"
78 #include "AbstractLinAlgPack_VectorDenseEncap.hpp"
79 #include "AbstractLinAlgPack_MatrixOpOut.hpp"
80 #include "AbstractLinAlgPack_VectorStdOps.hpp"
81 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
82 #include "DenseLinAlgPack_LinAlgOpPack.hpp"
83 #include "DenseLinAlgPack_DMatrixOut.hpp"
84 #include "DenseLinAlgLAPack.hpp"
85 #include "Teuchos_Workspace.hpp"
86 #include "Teuchos_Assert.hpp"
87 
88 namespace {
89 
90  using DenseLinAlgPack::DVectorSlice;
92 
100  void comp_Cb( const DMatrixSlice& Lb, const DVectorSlice& Db_diag
101  , DMatrixSlice* Cb );
102 
103 } // end namespace
104 
105 namespace ConstrainedOptPack {
106 
107 // /////////////////////////////////
108 // Inline private member functions
109 
110 inline
111 const DMatrixSliceTri MatrixSymPosDefLBFGS::R() const
112 {
113  return DenseLinAlgPack::tri( STY_(1,m_bar_,1,m_bar_), BLAS_Cpp::upper, BLAS_Cpp::nonunit );
114 }
115 
116 inline
117 const DMatrixSliceTri MatrixSymPosDefLBFGS::Lb() const
118 {
119  return DenseLinAlgPack::tri( STY_(2,m_bar_,1,m_bar_-1), BLAS_Cpp::lower, BLAS_Cpp::nonunit );
120 }
121 
122 inline
123 DMatrixSlice MatrixSymPosDefLBFGS::STY()
124 {
125  return STY_(1,m_bar_,1,m_bar_);
126 }
127 
128 inline
129 const DMatrixSlice MatrixSymPosDefLBFGS::STY() const
130 {
131  return STY_(1,m_bar_,1,m_bar_);
132 }
133 
134 inline
135 DMatrixSliceSym MatrixSymPosDefLBFGS::STS()
136 {
137  return DenseLinAlgPack::nonconst_sym( STSYTY_(2,m_bar_+1,1,m_bar_),BLAS_Cpp::lower );
138 }
139 
140 inline
141 const DMatrixSliceSym MatrixSymPosDefLBFGS::STS() const
142 {
143  return DenseLinAlgPack::sym( STSYTY_(2,m_bar_+1,1,m_bar_),BLAS_Cpp::lower );
144 }
145 
146 inline
147 DMatrixSliceSym MatrixSymPosDefLBFGS::YTY()
148 {
149  return DenseLinAlgPack::nonconst_sym( STSYTY_(1,m_bar_,2,m_bar_+1),BLAS_Cpp::upper );
150 }
151 
152 inline
153 const DMatrixSliceSym MatrixSymPosDefLBFGS::YTY() const
154 {
155  return DenseLinAlgPack::sym( STSYTY_(1,m_bar_,2,m_bar_+1),BLAS_Cpp::upper );
156 }
157 
158 // ///////////////////////
159 // Nonlinined functions
160 
162  size_type m
163  ,bool maintain_original
164  ,bool maintain_inverse
165  ,bool auto_rescaling
166  )
167 {
168  initial_setup(m,maintain_original,maintain_inverse,auto_rescaling);
169 }
170 
172  size_type m,
173  bool maintain_original,
174  bool maintain_inverse,
175  bool auto_rescaling
176  )
177 {
178  // Validate input
180  !maintain_original && !maintain_inverse, std::invalid_argument
181  ,"MatrixSymPosDefLBFGS::initial_setup(...) : "
182  "Error, both maintain_original and maintain_inverse can not both be false!" );
184  m < 1, std::invalid_argument
185  ,"MatrixSymPosDefLBFGS::set_num_updates_stored(m) : "
186  "Error, the number of storage locations must be > 0" );
187  vec_spc_ = Teuchos::null;
188  maintain_original_ = maintain_original;
189  maintain_inverse_ = maintain_inverse;
190  m_ = m;
191  n_ = 0; // make uninitialized
192  num_secant_updates_= 0;
193  auto_rescaling_ = auto_rescaling;
194 }
195 
196 // Overridden from MatrixOp
197 
198 const VectorSpace& MatrixSymPosDefLBFGS::space_cols() const
199 {
200  return *vec_spc_;
201 }
202 
203 std::ostream& MatrixSymPosDefLBFGS::output(std::ostream& out) const
204 {
205  assert_initialized();
206  out << "\n*** Limited Memory BFGS matrix\n";
207  out << "\nConversion to dense =\n";
208  MatrixOp::output(out);
209  out << "\nStored quantities\n"
210  << "\nn = " << n_
211  << "\nm = " << m_
212  << "\nm_bar = " << m_bar_
213  << "\ngamma_k = " << gamma_k_ << std::endl;
214  if( m_bar_ ) {
215  out << "\nS =\n" << *S()
216  << "\nY =\n" << *Y()
217  << "\nS'Y =\n" << STY_(1,m_bar_,1,m_bar_)
218  << "\nlower(S'S) \\ zero diagonal \\ upper(Y'Y) =\n"
219  << STSYTY_(1,m_bar_+1,1,m_bar_+1)
220  << "\nQ updated? = " << Q_updated_ << std::endl;
221  if(Q_updated_)
222  out << "\nCholesky of schur complement of Q, QJ =\n" << QJ_(1,m_bar_,1,m_bar_);
223  }
224  return out;
225 }
226 
227 MatrixOp& MatrixSymPosDefLBFGS::operator=(const MatrixOp& mwo)
228 {
229  const MatrixSymPosDefLBFGS *p_m = dynamic_cast<const MatrixSymPosDefLBFGS*>(&mwo);
230  if(p_m) {
231  if( p_m == this ) return *this; // assignment to self
232  // Important: Assign all members here.
233  auto_rescaling_ = p_m->auto_rescaling_;
234  maintain_original_ = p_m->maintain_original_;
235  original_is_updated_ = p_m->original_is_updated_;
236  maintain_inverse_ = p_m->maintain_inverse_;
237  inverse_is_updated_ = p_m->inverse_is_updated_;
238  vec_spc_ = p_m->vec_spc_.get() ? p_m->vec_spc_->clone() : Teuchos::null;
239  n_ = p_m->n_;
240  m_ = p_m->m_;
241  m_bar_ = p_m->m_bar_;
242  num_secant_updates_ = p_m->num_secant_updates_;
243  gamma_k_ = p_m->gamma_k_;
244  S_ = p_m->S_;
245  Y_ = p_m->Y_;
246  STY_ = p_m->STY_;
247  STSYTY_ = p_m->STSYTY_;
248  Q_updated_ = p_m->Q_updated_;
249  QJ_ = p_m->QJ_;
250  }
251  else {
253  true,std::invalid_argument
254  ,"MatrixSymPosDefLBFGS::operator=(const MatrixOp& mwo) : Error, "
255  "The concrete type of mwo \'" << typeName(mwo) << "\' is not "
256  "as subclass of MatrixSymPosDefLBFGS as required" );
257  }
258  return *this;
259 }
260 
261 // Level-2 BLAS
262 
264  VectorMutable* y, value_type alpha, BLAS_Cpp::Transp trans_rhs1
265  , const Vector& x, value_type beta
266  ) const
267 {
271  using LinAlgOpPack::V_StMtV;
272  using LinAlgOpPack::V_MtV;
273  typedef VectorDenseEncap vde;
274  typedef VectorDenseMutableEncap vdme;
275  using Teuchos::Workspace;
277 
278  assert_initialized();
279 
280  TEUCHOS_TEST_FOR_EXCEPT( !( original_is_updated_ ) ); // For now just always update
281 
282  // y = b*y + Bk * x
283  //
284  // y = b*y + (1/gk)*x - [ (1/gk)*S Y ] * inv(Q) * [ (1/gk)*S' ] * x
285  // [ Y' ]
286  // Perform the following operations (in order):
287  //
288  // y = b*y
289  //
290  // y += (1/gk)*x
291  //
292  // t1 = [ (1/gk)*S'*x ] <: R^(2*m)
293  // [ Y'*x ]
294  //
295  // t2 = inv(Q) * t1 <: R^(2*m)
296  //
297  // y += -(1/gk) * S * t2(1:m)
298  //
299  // y += -1.0 * Y * t2(m+1,2m)
300 
301  const value_type
302  invgk = 1.0 / gamma_k_;
303 
304  // y = b*y
305  Vt_S( y, beta );
306 
307  // y += (1/gk)*x
308  Vp_StV( y, invgk, x );
309 
310  if( !m_bar_ )
311  return; // No updates have been added yet.
312 
313  const multi_vec_ptr_t
314  S = this->S(),
315  Y = this->Y();
316 
317  // Get workspace
318 
319  const size_type
320  mb = m_bar_;
321 
322  Workspace<value_type> t1_ws(wss,2*mb);
323  DVectorSlice t1(&t1_ws[0],t1_ws.size());
324  Workspace<value_type> t2_ws(wss,2*mb);
325  DVectorSlice t2(&t2_ws[0],t2_ws.size());
326 
327  VectorSpace::vec_mut_ptr_t
328  t = S->space_rows().create_member();
329 
330  // t1 = [ (1/gk)*S'*x ]
331  // [ Y'*x ]
332 
333  V_StMtV( t.get(), invgk, *S, BLAS_Cpp::trans, x );
334  t1(1,mb) = vde(*t)();
335  V_MtV( t.get(), *Y, BLAS_Cpp::trans, x );
336  t1(mb+1,2*mb) = vde(*t)();
337 
338  // t2 = inv(Q) * t1
339  V_invQtV( &t2, t1 );
340 
341  // y += -(1/gk) * S * t2(1:m)
342  (vdme(*t)() = t2(1,mb));
343  Vp_StMtV( y, -invgk, *S, BLAS_Cpp::no_trans, *t );
344 
345  // y += -1.0 * Y * t2(m+1,2m
346  (vdme(*t)() = t2(mb+1,2*mb));
347  Vp_StMtV( y, -1.0, *Y, BLAS_Cpp::no_trans, *t );
348 
349 }
350 
351 // Overridden from MatrixOpNonsing
352 
354  VectorMutable* y, BLAS_Cpp::Transp trans_rhs1
355  , const Vector& x
356  ) const
357 {
360  using LinAlgOpPack::V_mV;
361  using LinAlgOpPack::V_StV;
362  using LinAlgOpPack::Vp_V;
363  using LinAlgOpPack::V_MtV;
364  using LinAlgOpPack::V_StMtV;
365  using LinAlgOpPack::Vp_MtV;
367  typedef VectorDenseEncap vde;
368  typedef VectorDenseMutableEncap vdme;
369  using Teuchos::Workspace;
371 
372  assert_initialized();
373 
374  TEUCHOS_TEST_FOR_EXCEPT( !( inverse_is_updated_ ) ); // For now just always update
375 
376  // y = inv(Bk) * x = Hk * x
377  //
378  // = gk*x + [S gk*Y] * [ inv(R')*(D+gk*Y'Y)*inv(R) -inv(R') ] * [ S' ] * x
379  // [ -inv(R) 0 ] [ gk*Y' ]
380  //
381  // Perform in the following (in order):
382  //
383  // y = gk*x
384  //
385  // t1 = [ S'*x ] <: R^(2*m)
386  // [ gk*Y'*x ]
387  //
388  // t2 = inv(R) * t1(1:m) <: R^(m)
389  //
390  // t3 = - inv(R') * t1(m+1,2*m) <: R^(m)
391  //
392  // t4 = gk * Y'Y * t2 <: R^(m)
393  //
394  // t4 += D*t2
395  //
396  // t5 = inv(R') * t4 <: R^(m)
397  //
398  // t5 += t3
399  //
400  // y += S*t5
401  //
402  // y += -gk*Y*t2
403 
404  // y = gk*x
405  V_StV( y, gamma_k_, x );
406 
407  const size_type
408  mb = m_bar_;
409 
410  if( !mb )
411  return; // No updates have been performed.
412 
413  const multi_vec_ptr_t
414  S = this->S(),
415  Y = this->Y();
416 
417  // Get workspace
418 
419  Workspace<value_type> t1_ws(wss,2*mb);
420  DVectorSlice t1(&t1_ws[0],t1_ws.size());
421  Workspace<value_type> t2_ws(wss,mb);
422  DVectorSlice t2(&t2_ws[0],t2_ws.size());
423  Workspace<value_type> t3_ws(wss,mb);
424  DVectorSlice t3(&t3_ws[0],t3_ws.size());
425  Workspace<value_type> t4_ws(wss,mb);
426  DVectorSlice t4(&t4_ws[0],t4_ws.size());
427  Workspace<value_type> t5_ws(wss,mb);
428  DVectorSlice t5(&t5_ws[0],t5_ws.size());
429 
430  VectorSpace::vec_mut_ptr_t
431  t = S->space_rows().create_member();
432 
433  const DMatrixSliceTri
434  &R = this->R();
435 
436  const DMatrixSliceSym
437  &YTY = this->YTY();
438 
439  // t1 = [ S'*x ]
440  // [ gk*Y'*x ]
441  V_MtV( t.get(), *S, BLAS_Cpp::trans, x );
442  t1(1,mb) = vde(*t)();
443  V_StMtV( t.get(), gamma_k_, *Y, BLAS_Cpp::trans, x );
444  t1(mb+1,2*mb) = vde(*t)();
445 
446  // t2 = inv(R) * t1(1:m)
447  V_InvMtV( &t2, R, BLAS_Cpp::no_trans, t1(1,mb) );
448 
449  // t3 = - inv(R') * t1(m+1,2*m)
450  V_mV( &t3, t1(mb+1,2*mb) );
451  V_InvMtV( &t3, R, BLAS_Cpp::trans, t3 );
452 
453  // t4 = gk * Y'Y * t2
454  V_StMtV( &t4, gamma_k_, YTY, BLAS_Cpp::no_trans, t2 );
455 
456  // t4 += D*t2
457  Vp_DtV( &t4, t2 );
458 
459  // t5 = inv(R') * t4
460  V_InvMtV( &t5, R, BLAS_Cpp::trans, t4 );
461 
462  // t5 += t3
463  Vp_V( &t5, t3 );
464 
465  // y += S*t5
466  (vdme(*t)() = t5);
467  Vp_MtV( y, *S, BLAS_Cpp::no_trans, *t );
468 
469  // y += -gk*Y*t2
470  (vdme(*t)() = t2);
471  Vp_StMtV( y, -gamma_k_, *Y, BLAS_Cpp::no_trans, *t );
472 
473 }
474 
475 // Overridden from MatrixSymSecant
476 
477 void MatrixSymPosDefLBFGS::init_identity( const VectorSpace& space_diag, value_type alpha )
478 {
479  // Validate input
481  alpha <= 0.0, std::invalid_argument
482  ,"MatrixSymPosDefLBFGS::init_identity(n,alpha) : Error, "
483  "alpha = " << alpha << " <= 0 is not allowed!" );
484 
485  // Set the vector space
486  vec_spc_ = space_diag.clone();
487  vec_spc_.get();
488 
489  // Set storage
490  S_ = vec_spc_->create_members(m_);
491  Y_ = vec_spc_->create_members(m_);
492  TEUCHOS_TEST_FOR_EXCEPT( !( S_.get() ) );
493  TEUCHOS_TEST_FOR_EXCEPT( !( Y_.get() ) );
494  STY_.resize( m_, m_ );
495  STSYTY_.resize( m_+1, m_+1 );
496  STSYTY_.diag(0) = 0.0;
497 
498  gamma_k_ = 1.0/alpha;
499 
500  // Initialize counters
501  m_bar_ = 0;
502 
503  n_ = vec_spc_->dim(); // initialized;
504  original_is_updated_ = true; // This will never change for now
505  inverse_is_updated_ = true; // This will never change for now
506  num_secant_updates_ = 0; // reset this to zero
507 }
508 
509 void MatrixSymPosDefLBFGS::init_diagonal(const Vector& diag)
510 {
511  init_identity( diag.space(), diag.norm_inf() );
512 }
513 
515  VectorMutable* s, VectorMutable* y, VectorMutable* Bs
516  )
517 {
518  using AbstractLinAlgPack::BFGS_sTy_suff_p_d;
520  using LinAlgOpPack::V_MtV;
521  using Teuchos::Workspace;
523 
524  assert_initialized();
525 
526  // Check skipping the BFGS update
527  const value_type
528  sTy = dot(*s,*y);
529  std::ostringstream omsg;
530  if( !BFGS_sTy_suff_p_d(*s,*y,&sTy,&omsg,"MatrixSymPosDefLBFGS::secant_update(...)" ) ) {
531  throw UpdateSkippedException( omsg.str() );
532  }
533 
534  try {
535 
536  // Update counters
537  if( m_bar_ == m_ ) {
538  // We are at the end of the storage so remove the oldest stored update
539  // and move updates to make room for the new update. This has to be done for the
540  // the matrix to behave properly
541  {for( size_type k = 1; k <= m_-1; ++k ) {
542  S_->col(k) = S_->col(k+1); // Shift S.col() to the left
543  Y_->col(k) = Y_->col(k+1); // Shift Y.col() to the left
544  STY_.col(k)(1,m_-1) = STY_.col(k+1)(2,m_); // Move submatrix STY(2,m-1,2,m-1) up and left
545  STSYTY_.col(k)(k+1,m_) = STSYTY_.col(k+1)(k+2,m_+1); // Move triangular submatrix STS(2,m-1,2,m-1) up and left
546  STSYTY_.col(k+1)(1,k) = STSYTY_.col(k+2)(2,k+1); // Move triangular submatrix YTY(2,m-1,2,m-1) up and left
547  }}
548  // ToDo: Create an abstract interface, call it MultiVectorShiftVecs, to rearrange S and Y all at once.
549  // This will be important for parallel performance.
550  }
551  else {
552  m_bar_++;
553  }
554  // Set the update vectors
555  *S_->col(m_bar_) = *s;
556  *Y_->col(m_bar_) = *y;
557 
558  // /////////////////////////////////////////////////////////////////////////////////////
559  // Update S'Y
560  //
561  // Update the row and column m_bar
562  //
563  // S'Y =
564  //
565  // [ s(1)'*y(1) ... s(1)'*y(m_bar) ... s(1)'*y(m_bar) ]
566  // [ . . . ] row
567  // [ s(m_bar)'*y(1) ... s(m_bar)'*y(m_bar) ... s(m_bar)'*y(m_bar) ] m_bar
568  // [ . . . ]
569  // [ s(m_bar)'*y(1) ... s(m_bar)'*y(m_bar) ... s(m_bar)'*y(m_bar) ]
570  //
571  // col m_bar
572  //
573  // Therefore we set:
574  // (S'Y)(:,m_bar) = S'*y(m_bar)
575  // (S'Y)(m_bar,:) = s(m_bar)'*Y
576 
577  const multi_vec_ptr_t
578  S = this->S(),
579  Y = this->Y();
580 
581  VectorSpace::vec_mut_ptr_t
582  t = S->space_rows().create_member(); // temporary workspace
583 
584  // (S'Y)(:,m_bar) = S'*y(m_bar)
585  V_MtV( t.get(), *S, BLAS_Cpp::trans, *y );
586  STY_.col(m_bar_)(1,m_bar_) = VectorDenseEncap(*t)();
587 
588  // (S'Y)(m_bar,:)' = Y'*s(m_bar)
589  V_MtV( t.get(), *Y, BLAS_Cpp::trans, *s );
590  STY_.row(m_bar_)(1,m_bar_) = VectorDenseEncap(*t)();
591 
592  // /////////////////////////////////////////////////////////////////
593  // Update S'S
594  //
595  // S'S =
596  //
597  // [ s(1)'*s(1) ... symmetric symmetric ]
598  // [ . . . ] row
599  // [ s(m_bar)'*s(1) ... s(m_bar)'*s(m_bar) ... symmetric ] m_bar
600  // [ . . . ]
601  // [ s(m_bar)'*s(1) ... s(m_bar)'*s(m_bar) ... s(m_bar)'*s(m_bar) ]
602  //
603  // col m_bar
604  //
605  // Here we will update the lower triangular part of S'S. To do this we
606  // only need to compute:
607  // t = S'*s(m_bar) = { s(m_bar)' * [ s(1),..,s(m_bar),..,s(m_bar) ] }'
608  // then set the appropriate rows and columns of S'S.
609 
610  Workspace<value_type> work_ws(wss,m_bar_);
611  DVectorSlice work(&work_ws[0],work_ws.size());
612 
613  // work = S'*s(m_bar)
614  V_MtV( t.get(), *S, BLAS_Cpp::trans, *s );
615  work = VectorDenseEncap(*t)();
616 
617  // Set row elements
618  STSYTY_.row(m_bar_+1)(1,m_bar_) = work;
619  // Set column elements
620  STSYTY_.col(m_bar_)(m_bar_+1,m_bar_+1) = work(m_bar_,m_bar_);
621 
622  // /////////////////////////////////////////////////////////////////////////////////////
623  // Update Y'Y
624  //
625  // Update the row and column m_bar
626  //
627  // Y'Y =
628  //
629  // [ y(1)'*y(1) ... y(1)'*y(m_bar) ... y(1)'*y(m_bar) ]
630  // [ . . . ] row
631  // [ symmetric ... y(m_bar)'*y(m_bar) ... y(m_bar)'*y(m_bar) ] m_bar
632  // [ . . . ]
633  // [ symmetric ... symmetric ... y(m_bar)'*y(m_bar) ]
634  //
635  // col m_bar
636  //
637  // Here we will update the upper triangular part of Y'Y. To do this we
638  // only need to compute:
639  // t = Y'*y(m_bar) = { y(m_bar)' * [ y(1),..,y(m_bar),..,y(m_bar) ] }'
640  // then set the appropriate rows and columns of Y'Y.
641 
642  // work = Y'*y(m_bar)
643  V_MtV( t.get(), *Y, BLAS_Cpp::trans, *y );
644  work = VectorDenseEncap(*t)();
645 
646  // Set row elements
647  STSYTY_.col(m_bar_+1)(1,m_bar_) = work;
648  // Set column elements
649  STSYTY_.row(m_bar_)(m_bar_+1,m_bar_+1) = work(m_bar_,m_bar_);
650 
651  // /////////////////////////////
652  // Update gamma_k
653 
654  // gamma_k = s'*y / y'*y
655  if(auto_rescaling_)
656  gamma_k_ = STY_(m_bar_,m_bar_) / STSYTY_(m_bar_,m_bar_+1);
657 
658  // We do not initially update Q unless we have to form a matrix-vector
659  // product later.
660 
661  Q_updated_ = false;
662  num_secant_updates_++;
663 
664  } // end try
665  catch(...) {
666  // If we throw any exception the we should make the matrix uninitialized
667  // so that we do not leave this object in an inconsistant state.
668  n_ = 0;
669  throw;
670  }
671 
672 }
673 
674 // Private member functions
675 
676 void MatrixSymPosDefLBFGS::Vp_DtV( DVectorSlice* y, const DVectorSlice& x ) const
677 {
678  DenseLinAlgPack::Vp_MtV_assert_sizes(
679  y->dim(), m_bar_, m_bar_, BLAS_Cpp::no_trans, x.dim() );
680 
681  DVectorSlice::const_iterator
682  d_itr = STY_.diag(0).begin(),
683  x_itr = x.begin();
684  DVectorSlice::iterator
685  y_itr = y->begin();
686 
687  while( y_itr != y->end() )
688  *y_itr++ += (*d_itr++) * (*x_itr++);
689 }
690 
691 //
692 // We need to update the factorizations to solve for:
693 //
694 // x = inv(Q) * y => Q * x = y
695 //
696 // [ (1/gk)*S'S L ] * [ x1 ] = [ y1 ]
697 // [ L' -D ] [ x2 ] [ y2 ]
698 //
699 // We will solve the above system using a schur complement:
700 //
701 // C = (1/gk)*S'S + L*inv(D)*L'
702 //
703 // According to the referenced paper, C is p.d. so:
704 //
705 // C = J*J'
706 //
707 // We then compute the solution as:
708 //
709 // x1 = inv(C) * ( y1 + L*inv(D)*y2 )
710 // x2 = - inv(D) * ( y2 - L'*x1 )
711 //
712 // Therefore we will just update the factorization C = J*J'
713 // where the factor J is stored in QJ_.
714 //
715 
716 void MatrixSymPosDefLBFGS::update_Q() const
717 {
718  using DenseLinAlgPack::tri;
719  using DenseLinAlgPack::tri_ele;
721 
722  //
723  // We need update the factorizations to solve for:
724  //
725  // x = inv(Q) * y
726  //
727  // [ y1 ] = [ (1/gk)*S'S L ] * [ x1 ]
728  // [ y2 ] [ L' -D ] [ x2 ]
729  //
730  // We will solve the above system using the schur complement:
731  //
732  // C = (1/gk)*S'S + L*inv(D)*L'
733  //
734  // According to the referenced paper, C is p.d. so:
735  //
736  // C = J*J'
737  //
738  // We then compute the solution as:
739  //
740  // x1 = inv(C) * ( y1 + L*inv(D)*y2 )
741  // x2 = - inv(D) * ( y2 - L'*x1 )
742  //
743  // Therefore we will just update the factorization C = J*J'
744  //
745 
746  // Form the upper triangular part of C which will become J
747  // which we are using storage of QJ
748 
749  if( QJ_.rows() < m_ )
750  QJ_.resize( m_, m_ );
751 
752  const size_type
753  mb = m_bar_;
754 
755  DMatrixSlice
756  C = QJ_(1,mb,1,mb);
757 
758  // C = L * inv(D) * L'
759  //
760  // Here L is a strictly lower triangular (zero diagonal) matrix where:
761  //
762  // L = [ 0 0 ]
763  // [ Lb 0 ]
764  //
765  // Lb is lower triangular (nonzero diagonal)
766  //
767  // Therefore we compute L*inv(D)*L' as:
768  //
769  // C = [ 0 0 ] * [ Db 0 ] * [ 0 Lb' ]
770  // [ Lb 0 ] [ 0 db ] [ 0 0 ]
771  //
772  // = [ 0 0 ] = [ 0 0 ]
773  // [ 0 Cb ] [ 0 Lb*Db*Lb' ]
774  //
775  // We need to compute the upper triangular part of Cb.
776 
777  C.row(1) = 0.0;
778  if( mb > 1 )
779  comp_Cb( STY_(2,mb,1,mb-1), STY_.diag(0)(1,mb-1), &C(2,mb,2,mb) );
780 
781  // C += (1/gk)*S'S
782 
783  const DMatrixSliceSym &STS = this->STS();
784  Mp_StM( &C, (1/gamma_k_), tri( STS.gms(), STS.uplo(), BLAS_Cpp::nonunit )
785  , BLAS_Cpp::trans );
786 
787  // Now perform a cholesky factorization of C
788  // After this factorization the upper triangular part of QJ
789  // (through C) will contain the cholesky factor.
790 
791  DMatrixSliceTriEle C_upper = tri_ele( C, BLAS_Cpp::upper );
792  try {
793  DenseLinAlgLAPack::potrf( &C_upper );
794  }
795  catch( const DenseLinAlgLAPack::FactorizationException &fe ) {
797  true, UpdateFailedException
798  ,"Error, the factorization of Q which should be s.p.d. failed with"
799  " the error message: {" << fe.what() << "}";
800  );
801  }
802 
803  Q_updated_ = true;
804 }
805 
806 void MatrixSymPosDefLBFGS::V_invQtV( DVectorSlice* x, const DVectorSlice& y ) const
807 {
808  using DenseLinAlgPack::sym;
809  using DenseLinAlgPack::tri;
812 
813  using LinAlgOpPack::Vp_V;
814  using LinAlgOpPack::V_MtV;
815 
816 
817  // Solve the system
818  //
819  // Q * x = y
820  //
821  // Using the schur complement factorization as described above.
822 
823  const size_type
824  mb = m_bar_;
825 
826  if(!Q_updated_) {
827  update_Q();
828  }
829 
830  DVectorSlice
831  x1 = (*x)(1,mb),
832  x2 = (*x)(mb+1,2*mb);
833 
834  const DVectorSlice
835  y1 = y(1,mb),
836  y2 = y(mb+1,2*mb);
837 
838  // //////////////////////////////////////
839  // x1 = inv(C) * ( y1 + L*inv(D)*y2 )
840  // = inv(J'*J) * r
841  // = inv(J) * inv(J') * r
842 
843  { // x1 = inv(D) * y2
844  DVectorSlice::const_iterator
845  d_itr = STY_.diag(0).begin(),
846  y2_itr = y2.begin();
847  DVectorSlice::iterator
848  x1_itr = x1.begin();
849  while( x1_itr != x1.end() )
850  *x1_itr++ = *y2_itr++ / *d_itr++;
851  }
852 
853  // x1 = L * x1
854  //
855  // = [ 0 0 ] * [ x1(1:mb-1) ]
856  // [ Lb 0 ] [ x1(mb) ]
857  //
858  // = [ 0 ]
859  // [ Lb*x1(1:mb-1) ]
860  //
861  if( mb > 1 ) {
862  // x1(2,mb) = x1(1,mb-1) ( copy from mb-1 to mb then mb-2 to mb-1
863  // etc. so that we don't overwrite the elements we need to copy ).
864  DVectorSlice
865  x1a = x1(1,mb-1),
866  x1b = x1(2,mb);
867  std::copy( x1a.rbegin(), x1a.rend(), x1b.rbegin() );
868  V_MtV( &x1b, Lb(), BLAS_Cpp::no_trans, x1b );
869  }
870  x1(1) = 0.0;
871 
872  // r = x1 += y1
873  Vp_V( &x1, y1 );
874 
875  // x1 = inv(J') * r
876  const DMatrixSliceTri J = tri( QJ_(1,mb,1,mb), BLAS_Cpp::upper, BLAS_Cpp::nonunit );
877  V_InvMtV( &x1, J, BLAS_Cpp::trans, x1 );
878 
879  // x1 = inv(J) * x1
880  V_InvMtV( &x1, J, BLAS_Cpp::no_trans, x1 );
881 
882  // /////////////////////////////////////
883  // x2 = inv(D) * ( - y2 + L'*x1 )
884 
885  // x2 = L'*x1
886  //
887  // = [ 0 Lb' ] * [ x1(1) ]
888  // [ 0 0 ] [ x1(2,mb) ]
889  //
890  // = [ Lb'*x1(2,mb) ]
891  // [ 0 ]
892  //
893  if( mb > 1 ) {
894  V_MtV( &x2(1,mb-1), Lb(), BLAS_Cpp::trans, x1(2,mb) );
895  }
896  x2(mb) = 0.0;
897 
898  // x2 += -y2
899  Vp_StV( &x2, -1.0, y2 );
900 
901  // x2 = inv(D) * x2
902  {
903  DVectorSlice::const_iterator
904  d_itr = STY_.diag(0).begin();
905  DVectorSlice::iterator
906  x2_itr = x2.begin();
907  while( x2_itr != x2.end() )
908  *x2_itr++ /= *d_itr++;
909  }
910 }
911 
912 void MatrixSymPosDefLBFGS::assert_initialized() const
913 {
914  if(!n_)
915  throw std::logic_error( "MatrixSymPosDefLBFGS::assert_initialized() : "
916  "Error, matrix not initialized" );
917 }
918 
919 } // end namespace ConstrainedOptPack
920 
921 namespace {
922 
923 void comp_Cb(
924  const DMatrixSlice& Lb, const DVectorSlice& Db_diag
925  ,DMatrixSlice* Cb
926  )
927 {
928  // Lb * inv(Db) * Lb =
929  //
930  // [ l(1,1) ] [ dd(1) ] [ l(1,1) l(2,1) ... l(p,1) ]
931  // [ l(2,1) l(2,2) ] [ dd(2) ] [ l(2,2) ... l(p,2) ]
932  // [ . . . ] * [ . ] * [ . . ]
933  // [ l(p,1) l(p,2) ... l(p,p) ] [ dd(p) ] [ l(p,p) ]
934  //
935  //
936  // [ l(1,1)*dd(1)*l(1,1) l(1,1)*dd(1)*l(2,1) ... l(1,1)*dd(1)*l(p,1) ]
937  // [ symmetric l(2,1)*dd(1)*l(2,1) + l(2,2)*dd(2)*l(2,2) ... l(2,1)*dd(1)*l(p,1) + l(2,2)*dd(2)*l(p,2) ]
938  // [ . . ...
939  // [ symmetric symmetric ... sum( l(p,i)*dd(i)*l(p,i), i=1,..,p ) ]
940  //
941  // Therefore we can express the upper triangular elemetns of Cb as:
942  //
943  // Cb(i,j) = sum( l(i,k)*dd(k)*l(j,k), k = 1,..,i )
944 
946  typedef DenseLinAlgPack::value_type value_type;
947 
948  TEUCHOS_TEST_FOR_EXCEPT( !( Lb.rows() == Cb->rows() && Cb->rows() == Db_diag.dim() ) ); // only a local error!
949 
950  const size_type p = Db_diag.dim();
951 
952  for( size_type i = 1; i <= p; ++i ) {
953  for( size_type j = i; j <= p; ++j ) {
954  value_type &c = (*Cb)(i,j) = 0.0;
955  for( size_type k = 1; k <= i; ++k ) {
956  c += Lb(i,k) * Lb(j,k) / Db_diag(k);
957  }
958  }
959  }
960 
961  // ToDo: Make the above operation more efficent if needed! (i.e. write
962  // it in fortran or something?).
963 }
964 
965 } // end namespace
void init_identity(const VectorSpace &space_diag, value_type alpha)
void secant_update(VectorMutable *s, VectorMutable *y, VectorMutable *Bs)
void initial_setup(size_type m=10, bool maintain_original=true, bool maintain_inverse=true, bool auto_rescaling=false)
Initial setup for the matrix.
void Vt_S(VectorMutable *v_lhs, const value_type &alpha)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
void Vp_StV(VectorMutable *v_lhs, const value_type &alpha, const Vector &v_rhs)
void V_InvMtV(VectorMutable *v_lhs, BLAS_Cpp::Transp trans_rhs1, const Vector &v_rhs2) const
void V_StV(VectorMutable *v_lhs, value_type alpha, const V &V_rhs)
MatrixSymPosDefLBFGS(size_type m=10, bool maintain_original=true, bool maintain_inverse=true, bool auto_rescaling=false)
Calls this->initial_setup()
void init_diagonal(const Vector &diag)
Actually this calls init_identity( diag.space(), diag.norm_inf() ).
void V_mV(VectorMutable *v_lhs, const V &V_rhs)
size_t size_type
void V_InvMtV(VectorMutable *v_lhs, const MatrixNonsing &M_rhs1, BLAS_Cpp::Transp trans_rhs1, const Vector &v_rhs2)
void Mp_StM(MatrixOp *mwo_lhs, value_type alpha, const MatrixOp &M_rhs, BLAS_Cpp::Transp trans_rhs)
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)
friend friend void Mp_StM(MatrixOp *mwo_lhs, value_type alpha, const MatrixOp &M_rhs, BLAS_Cpp::Transp trans_rhs)
value_type dot(const Vector &v_rhs1, const Vector &v_rhs2)
void Vp_StMtV(VectorMutable *v_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs1, const Vector &v_rhs2, value_type beta) const
Implementation of limited Memory BFGS matrix for arbitrary vector spaces.
Transp
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
TEUCHOSCORE_LIB_DLL_EXPORT Teuchos::RCP< WorkspaceStore > get_default_workspace_store()
std::string typeName(const T &t)
void Vp_V(VectorMutable *v_lhs, const V &V_rhs)