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_MatrixSymPosDefLBFGSSerial.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 // Let's define a compact representation for the matrix B^{k} and
45 // its inverse H^{k} = inv(B^{k}).
46 //
47 // Bk = (1/gk)*I - [ (1/gk)*S Y ] * inv( [ (1/gk)*S'S L ] [ (1/gk)*S' ]
48 // [ L' -D ] ) * [ Y' ]
49 // \________________/
50 // Q
51 //
52 // Hk = gk*I + [ S gk*Y ] * [ inv(R')*(D+gk*Y'Y)*inv(R) -inv(R') ] * [ S' ]
53 // [ -inv(R) 0 ] [ gk*Y' ]
54 //
55 // where:
56 //
57 // gk = gamma_k <: R
58 //
59 // [ s^{1}'*y^{1} s^{1}'*y^{2} ... s^{1}'*y^{m} ]
60 // S'Y = [ s^{2}'*y^{1} s^{2}'*y^{2} ... s^{2}'*y^{m} ] <: R^(m x m)
61 // [ . . . ]
62 // [ s^{m}'*y^{1} s^{m}'*y^{2} ... s^{m}'*y^{m} ]
63 //
64 // [ s^{1}'*y^{1} 0 ... 0 ]
65 // D = [ 0 s^{2}'*y^{2} ... 0 ] <: R^(m x m)
66 // [ . . . ]
67 // [ 0 0 ... s^{m}'*y^{m} ]
68 //
69 // R = upper triangular part of S'Y
70 //
71 // L = lower tirangular part of S'Y with zeros on the diagonal
72 //
73 
74 #include <assert.h>
75 
76 #include "ConstrainedOptPack_MatrixSymPosDefLBFGS.hpp"
77 #include "ConstrainedOptPack/src/AbstractLinAlgPack_BFGS_helpers.hpp"
78 #include "DenseLinAlgPack_LinAlgOpPack.hpp"
79 #include "DenseLinAlgPack_DMatrixOut.hpp"
80 #include "DenseLinAlgLAPack.hpp"
81 
82 namespace {
83 
84  using DenseLinAlgPack::DVectorSlice;
86 
94  void comp_Cb( const DMatrixSlice& Lb, const DVectorSlice& Db_diag
95  , DMatrixSlice* Cb );
96 
97 } // end namespace
98 
99 namespace ConstrainedOptPack {
100 
101 // /////////////////////////////////
102 // Inline private member functions
103 
104 inline
105 const DMatrixSliceTri MatrixSymPosDefLBFGS::R() const
106 {
107  return DenseLinAlgPack::tri( STY_(1,m_bar_,1,m_bar_), BLAS_Cpp::upper, BLAS_Cpp::nonunit );
108 }
109 
110 inline
111 const DMatrixSliceTri MatrixSymPosDefLBFGS::Lb() const
112 {
113  return DenseLinAlgPack::tri( STY_(2,m_bar_,1,m_bar_-1), BLAS_Cpp::lower, BLAS_Cpp::nonunit );
114 }
115 
116 inline
117 DMatrixSlice MatrixSymPosDefLBFGS::STY()
118 {
119  return STY_(1,m_bar_,1,m_bar_);
120 }
121 
122 inline
123 const DMatrixSlice MatrixSymPosDefLBFGS::STY() const
124 {
125  return STY_(1,m_bar_,1,m_bar_);
126 }
127 
128 inline
129 DMatrixSliceSym MatrixSymPosDefLBFGS::STS()
130 {
131  return DenseLinAlgPack::nonconst_sym( STSYTY_(2,m_bar_+1,1,m_bar_),BLAS_Cpp::lower );
132 }
133 
134 inline
135 const DMatrixSliceSym MatrixSymPosDefLBFGS::STS() const
136 {
137  return DenseLinAlgPack::sym( STSYTY_(2,m_bar_+1,1,m_bar_),BLAS_Cpp::lower );
138 }
139 
140 inline
141 DMatrixSliceSym MatrixSymPosDefLBFGS::YTY()
142 {
143  return DenseLinAlgPack::nonconst_sym( STSYTY_(1,m_bar_,2,m_bar_+1),BLAS_Cpp::upper );
144 }
145 
146 inline
147 const DMatrixSliceSym MatrixSymPosDefLBFGS::YTY() const
148 {
149  return DenseLinAlgPack::sym( STSYTY_(1,m_bar_,2,m_bar_+1),BLAS_Cpp::upper );
150 }
151 
152 // ///////////////////////
153 // Nonlinined functions
154 
156  size_type max_size
157  ,size_type m
158  ,bool maintain_original
159  ,bool maintain_inverse
160  ,bool auto_rescaling
161  )
162 {
163  initial_setup(max_size,m,maintain_original,maintain_inverse,auto_rescaling);
164 }
165 
167  size_type max_size
168  ,size_type m
169  ,bool maintain_original
170  ,bool maintain_inverse
171  ,bool auto_rescaling
172  )
173 {
174  // Validate input
175  if( !maintain_original && !maintain_inverse )
176  throw std::invalid_argument(
177  "MatrixSymPosDefLBFGS::initial_setup(...) : "
178  "Error, both maintain_original and maintain_inverse can not both be false!" );
179  if( m < 1 )
180  throw std::invalid_argument(
181  "MatrixSymPosDefLBFGS::set_num_updates_stored(m) : "
182  "Error, the number of storage locations must be > 0" );
183  maintain_original_ = maintain_original;
184  maintain_inverse_ = maintain_inverse;
185  m_ = m;
186  n_ = 0; // make uninitialized
187  n_max_ = max_size;
188  num_secant_updates_= 0;
189 }
190 
191 // Overridden from Matrix
192 
194 {
195  return n_;
196 }
197 
198 // Overridden from MatrixOp
199 
200 std::ostream& MatrixSymPosDefLBFGS::output(std::ostream& out) const
201 {
202  assert_initialized();
203  out << "*** Limited Memory BFGS matrix.\n"
204  << "Conversion to dense =\n";
205  MatrixOp::output(out);
206  out << "\n*** Stored quantities\n"
207  << "\ngamma_k = " << gamma_k_ << std::endl;
208  if( m_bar_ ) {
209  out << "\nS =\n" << S()
210  << "\nY =\n" << Y()
211  << "\nS'Y =\n" << STY_(1,m_bar_,1,m_bar_)
212  << "\nlower(S'S) \\ zero diagonal \\ upper(Y'Y) =\n"
213  << STSYTY_(1,m_bar_+1,1,m_bar_+1)
214  << "\nQ updated? = " << Q_updated_ << std::endl
215  << "\nCholesky of schur complement of Q, QJ =\n" << QJ_(1,m_bar_,1,m_bar_);
216  }
217  return out;
218 }
219 
220 MatrixOp& MatrixSymPosDefLBFGS::operator=(const MatrixOp& m)
221 {
222  const MatrixSymPosDefLBFGS *p_m = dynamic_cast<const MatrixSymPosDefLBFGS*>(&m);
223  if(p_m) {
224  if( p_m == this ) return *this; // assignment to self
225  // Important: Assign all members here.
226  auto_rescaling_ = p_m->auto_rescaling_;
227  maintain_original_ = p_m->maintain_original_;
228  original_is_updated_ = p_m->original_is_updated_;
229  maintain_inverse_ = p_m->maintain_inverse_;
230  inverse_is_updated_ = p_m->inverse_is_updated_;
231  n_max_ = p_m->n_max_;
232  n_ = p_m->n_;
233  m_ = p_m->m_;
234  m_bar_ = p_m->m_bar_;
235  k_bar_ = p_m->k_bar_;
236  num_secant_updates_ = p_m->num_secant_updates_;
237  gamma_k_ = p_m->gamma_k_;
238  S_ = p_m->S_;
239  Y_ = p_m->Y_;
240  STY_ = p_m->STY_;
241  STSYTY_ = p_m->STSYTY_;
242  Q_updated_ = p_m->Q_updated_;
243  QJ_ = p_m->QJ_;
244  }
245  else {
246  throw std::invalid_argument("MatrixSymPosDefLBFGS::operator=(const MatrixOp& m)"
247  " : The concrete type of m is not a subclass of MatrixSymPosDefLBFGS as expected" );
248  }
249  return *this;
250 }
251 
252 // Level-2 BLAS
253 
255  DVectorSlice* y, value_type alpha, BLAS_Cpp::Transp trans_rhs1
256  , const DVectorSlice& x, value_type beta) const
257 {
258  using LinAlgOpPack::V_StMtV;
259  using LinAlgOpPack::V_MtV;
260 
261  using DenseLinAlgPack::Vt_S;
264 
265  assert_initialized();
266 
267  TEUCHOS_TEST_FOR_EXCEPT( !( original_is_updated_ ) ); // For now just always update
268 
269  // y = b*y + Bk * x
270  //
271  // y = b*y + (1/gk)*x - [ (1/gk)*S Y ] * inv(Q) * [ (1/gk)*S' ] * x
272  // [ Y' ]
273  // Perform the following operations (in order):
274  //
275  // y = b*y
276  //
277  // y += (1/gk)*x
278  //
279  // t1 = [ (1/gk)*S'*x ] <: R^(2*m)
280  // [ Y'*x ]
281  //
282  // t2 = inv(Q) * t1 <: R^(2*m)
283  //
284  // y += -(1/gk) * S * t2(1:m)
285  //
286  // y += -1.0 * Y * t2(m+1,2m)
287 
288  const value_type
289  invgk = 1.0 / gamma_k_;
290 
291  // y = b*y
292 
293  if( beta == 0.0 )
294  *y = beta;
295  else
296  Vt_S( y, beta );
297 
298  // y += (1/gk)*x
299 
300  Vp_StV( y, invgk, x );
301 
302  if( !m_bar_ )
303  return; // No updates have been added yet.
304 
305  // Get workspace
306 
307  if( work_.size() < 4 * m_ )
308  work_.resize( 4 * m_ );
309 
310  const size_type
311  mb = m_bar_;
312 
313  const size_type
314  t1s = 1,
315  t1n = 2*mb,
316  t2s = t1s+t1n,
317  t2n = 2*mb;
318 
319  DVectorSlice
320  t1 = work_( t1s, t1s + t1n - 1 ),
321  t2 = work_( t2s, t2s + t2n - 1 );
322 
323  const DMatrixSlice
324  &S = this->S(),
325  &Y = this->Y();
326 
327  // t1 = [ (1/gk)*S'*x ]
328  // [ Y'*x ]
329 
330  V_StMtV( &t1(1,mb), invgk, S, BLAS_Cpp::trans, x );
331  V_MtV( &t1(mb+1,2*mb), Y, BLAS_Cpp::trans, x );
332 
333  // t2 = inv(Q) * t1
334 
335  V_invQtV( &t2, t1 );
336 
337  // y += -(1/gk) * S * t2(1:m)
338 
339  Vp_StMtV( y, -invgk, S, BLAS_Cpp::no_trans, t2(1,mb) );
340 
341  // y += -1.0 * Y * t2(m+1,2m)
342 
343  Vp_StMtV( y, -1.0, Y, BLAS_Cpp::no_trans, t2(mb+1,2*mb) );
344 
345 }
346 
347 // Overridden from MatrixWithOpFactorized
348 
349 // Level-2 BLAS
350 
351 void MatrixSymPosDefLBFGS::V_InvMtV( DVectorSlice* y, BLAS_Cpp::Transp trans_rhs1
352  , const DVectorSlice& x ) const
353 {
354  using DenseLinAlgPack::V_mV;
357 
358  using LinAlgOpPack::Vp_V;
359  using LinAlgOpPack::V_MtV;
360  using LinAlgOpPack::V_StMtV;
361  using LinAlgOpPack::Vp_MtV;
363 
364  assert_initialized();
365 
366  TEUCHOS_TEST_FOR_EXCEPT( !( inverse_is_updated_ ) ); // For now just always update
367 
368  // y = inv(Bk) * x = Hk * x
369  //
370  // = gk*x + [S gk*Y] * [ inv(R')*(D+gk*Y'Y)*inv(R) -inv(R') ] * [ S' ] * x
371  // [ -inv(R) 0 ] [ gk*Y' ]
372  //
373  // Perform in the following (in order):
374  //
375  // y = gk*x
376  //
377  // t1 = [ S'*x ] <: R^(2*m)
378  // [ gk*Y'*x ]
379  //
380  // t2 = inv(R) * t1(1:m) <: R^(m)
381  //
382  // t3 = - inv(R') * t1(m+1,2*m) <: R^(m)
383  //
384  // t4 = gk * Y'Y * t2 <: R^(m)
385  //
386  // t4 += D*t2
387  //
388  // t5 = inv(R') * t4 <: R^(m)
389  //
390  // t5 += t3
391  //
392  // y += S*t5
393  //
394  // y += -gk*Y*t2
395 
396  // y = gk*x
397  V_StV( y, gamma_k_, x );
398 
399  const size_type
400  mb = m_bar_;
401 
402  if( !mb )
403  return; // No updates have been performed.
404 
405  // Get workspace
406 
407  if( work_.size() < 6*m_ )
408  work_.resize( 6*m_ );
409 
410  const size_type
411  t1s = 1,
412  t1n = 2*mb,
413  t2s = t1s + t1n,
414  t2n = mb,
415  t3s = t2s + t2n,
416  t3n = mb,
417  t4s = t3s + t3n,
418  t4n = mb,
419  t5s = t4s + t4n,
420  t5n = mb;
421 
422  DVectorSlice
423  t1 = work_( t1s, t1s + t1n - 1 ),
424  t2 = work_( t2s, t2s + t2n - 1 ),
425  t3 = work_( t3s, t3s + t3n - 1 ),
426  t4 = work_( t4s, t4s + t4n - 1 ),
427  t5 = work_( t5s, t5s + t5n - 1 );
428 
429  const DMatrixSlice
430  &S = this->S(),
431  &Y = this->Y();
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( &t1(1,mb), S, BLAS_Cpp::trans, x );
442  V_StMtV( &t1(mb+1,2*mb), gamma_k_, Y, BLAS_Cpp::trans, x );
443 
444  // t2 = inv(R) * t1(1:m)
445  V_InvMtV( &t2, R, BLAS_Cpp::no_trans, t1(1,mb) );
446 
447  // t3 = - inv(R') * t1(m+1,2*m)
448  V_mV( &t3, t1(mb+1,2*mb) );
449  V_InvMtV( &t3, R, BLAS_Cpp::trans, t3 );
450 
451  // t4 = gk * Y'Y * t2
452  V_StMtV( &t4, gamma_k_, YTY, BLAS_Cpp::no_trans, t2 );
453 
454  // t4 += D*t2
455  Vp_DtV( &t4, t2 );
456 
457  // t5 = inv(R') * t4
458  V_InvMtV( &t5, R, BLAS_Cpp::trans, t4 );
459 
460  // t5 += t3
461  Vp_V( &t5, t3 );
462 
463  // y += S*t5
464  Vp_MtV( y, S, BLAS_Cpp::no_trans, t5 );
465 
466  // y += -gk*Y*t2
467  Vp_StMtV( y, -gamma_k_, Y, BLAS_Cpp::no_trans, t2 );
468 
469 }
470 
471 // Overridden from MatrixSymSecant
472 
473 void MatrixSymPosDefLBFGS::init_identity(size_type n, value_type alpha)
474 {
475  // Validate input
476  if( alpha <= 0.0 ) {
477  std::ostringstream omsg;
478  omsg
479  << "MatrixSymPosDefLBFGS::init_identity(n,alpha) : Error, "
480  << "alpha = " << alpha << " <= 0 is not allowed!";
481  throw std::invalid_argument( omsg.str() );
482  }
483  if( n_max_ == 0 ) {
484  n_max_ = n;
485  }
486  else if( n > n_max_ ) {
487  std::ostringstream omsg;
488  omsg
489  << "MatrixSymPosDefLBFGS::init_identity(n,alpha) : Error, "
490  << "n = " << n << " > max_size = " << n_max_;
491  throw std::invalid_argument( omsg.str() );
492  }
493 
494  // Resize storage
495  S_.resize( n_max_, m_ );
496  Y_.resize( n_max_, m_ );
497  STY_.resize( m_, m_ );
498  STSYTY_.resize( m_+1, m_+1 );
499  STSYTY_.diag(0) = 0.0;
500 
501  gamma_k_ = 1.0/alpha;
502 
503  // Initialize counters
504  k_bar_ = 0;
505  m_bar_ = 0;
506 
507  n_ = n; // initialized;
508  original_is_updated_ = true; // This will never change for now
509  inverse_is_updated_ = true; // This will never change for now
510  num_secant_updates_ = 0;
511 }
512 
513 void MatrixSymPosDefLBFGS::init_diagonal(const DVectorSlice& diag)
514 {
515  using DenseLinAlgPack::norm_inf;
516  init_identity( diag.size(), norm_inf(diag) );
517 }
518 
520  DVectorSlice* s, DVectorSlice* y, DVectorSlice* Bs)
521 {
522  using DenseLinAlgPack::dot;
523  using DenseLinAlgPack::norm_2;
524 
525  using LinAlgOpPack::V_MtV;
526 
527  assert_initialized();
528 
529  // Check skipping the BFGS update
530  const value_type
531  sTy = dot(*s,*y);
532  std::ostringstream omsg;
533  if( !BFGS_sTy_suff_p_d(*s,*y,&sTy,&omsg,"MatrixSymPosDefLBFGS::secant_update(...)" ) ) {
534  throw UpdateSkippedException( omsg.str() );
535  }
536 
537  try {
538 
539  // Update counters
540  if( k_bar_ == m_ ) {
541 // // We are at the end storage so loop back around again
542 // k_bar_ = 1;
543  // We are at the end of the storage so remove the oldest stored update
544  // and move updates to make room for the new update. This has to be done for the
545  // the matrix to behave properly
546  {for( size_type k = 1; k <= m_-1; ++k ) {
547  S_.col(k) = S_.col(k+1); // Shift S.col() to the left
548  Y_.col(k) = Y_.col(k+1); // Shift Y.col() to the left
549  STY_.col(k)(1,m_-1) = STY_.col(k+1)(2,m_); // Move submatrix STY(2,m-1,2,m-1) up and left
550  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
551  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
552  }}
553  }
554  else {
555  k_bar_++;
556  }
557  if( m_bar_ < m_ ) {
558  // This is the first few updates where we have not maxed out the storage.
559  m_bar_++;
560  }
561 
562  // Set the update vectors
563  S_.col(k_bar_)(1,n_) = *s;
564  Y_.col(k_bar_)(1,n_) = *y;
565 
566  // /////////////////////////////////////////////////////////////////////////////////////
567  // Update S'Y
568  //
569  // Update the row and column k_bar
570  //
571  // S'Y =
572  //
573  // [ s(1)'*y(1) ... s(1)'*y(k_bar) ... s(1)'*y(m_bar) ]
574  // [ . . . ] row
575  // [ s(k_bar)'*y(1) ... s(k_bar)'*y(k_bar) ... s(k_bar)'*y(m_bar) ] k_bar
576  // [ . . . ]
577  // [ s(m_bar)'*y(1) ... s(m_bar)'*y(k_bar) ... s(m_bar)'*y(m_bar) ]
578  //
579  // col k_bar
580  //
581  // Therefore we set:
582  // (S'Y)(:,k_bar) = S'*y(k_bar)
583  // (S'Y)(k_bar,:) = s(k_bar)'*Y
584 
585  const DMatrixSlice
586  &S = this->S(),
587  &Y = this->Y();
588 
589  // (S'Y)(:,k_bar) = S'*y(k_bar)
590  V_MtV( &STY_.col(k_bar_)(1,m_bar_), S, BLAS_Cpp::trans, Y.col(k_bar_) );
591 
592  // (S'Y)(k_bar,:)' = Y'*s(k_bar)
593  V_MtV( &STY_.row(k_bar_)(1,m_bar_), Y, BLAS_Cpp::trans, S.col(k_bar_) );
594 
595  // /////////////////////////////////////////////////////////////////
596  // Update S'S
597  //
598  // S'S =
599  //
600  // [ s(1)'*s(1) ... symmetric symmetric ]
601  // [ . . . ] row
602  // [ s(k_bar)'*s(1) ... s(k_bar)'*s(k_bar) ... symmetric ] k_bar
603  // [ . . . ]
604  // [ s(m_bar)'*s(1) ... s(m_bar)'*s(k_bar) ... s(m_bar)'*s(m_bar) ]
605  //
606  // col k_bar
607  //
608  // Here we will update the lower triangular part of S'S. To do this we
609  // only need to compute:
610  // t = S'*s(k_bar) = { s(k_bar)' * [ s(1),..,s(k_bar),..,s(m_bar) ] }'
611  // then set the appropriate rows and columns of S'S.
612 
613  if( work_.size() < m_ )
614  work_.resize(m_);
615 
616  // work = S'*s(k_bar)
617  V_MtV( &work_(1,m_bar_), S, BLAS_Cpp::trans, S.col(k_bar_) );
618 
619  // Set row elements
620  STSYTY_.row(k_bar_+1)(1,k_bar_) = work_(1,k_bar_);
621  // Set column elements
622  STSYTY_.col(k_bar_)(k_bar_+1,m_bar_+1) = work_(k_bar_,m_bar_);
623 
624  // /////////////////////////////////////////////////////////////////////////////////////
625  // Update Y'Y
626  //
627  // Update the row and column k_bar
628  //
629  // Y'Y =
630  //
631  // [ y(1)'*y(1) ... y(1)'*y(k_bar) ... y(1)'*y(m_bar) ]
632  // [ . . . ] row
633  // [ symmetric ... y(k_bar)'*y(k_bar) ... y(k_bar)'*y(m_bar) ] k_bar
634  // [ . . . ]
635  // [ symmetric ... symmetric ... y(m_bar)'*y(m_bar) ]
636  //
637  // col k_bar
638  //
639  // Here we will update the upper triangular part of Y'Y. To do this we
640  // only need to compute:
641  // t = Y'*y(k_bar) = { y(k_bar)' * [ y(1),..,y(k_bar),..,y(m_bar) ] }'
642  // then set the appropriate rows and columns of Y'Y.
643 
644  // work = Y'*y(k_bar)
645  V_MtV( &work_(1,m_bar_), Y, BLAS_Cpp::trans, Y.col(k_bar_) );
646 
647  // Set row elements
648  STSYTY_.col(k_bar_+1)(1,k_bar_) = work_(1,k_bar_);
649  // Set column elements
650  STSYTY_.row(k_bar_)(k_bar_+1,m_bar_+1) = work_(k_bar_,m_bar_);
651 
652  // /////////////////////////////
653  // Update gamma_k
654 
655  // gamma_k = s'*y / y'*y
656  if(auto_rescaling_)
657  gamma_k_ = STY_(k_bar_,k_bar_) / STSYTY_(k_bar_,k_bar_+1);
658 
659  // We do not initially update Q unless we have to form a matrix-vector
660  // product later.
661 
662  Q_updated_ = false;
663  num_secant_updates_++;
664 
665  } // end try
666  catch(...) {
667  // If we throw any exception the we should make the matrix uninitialized
668  // so that we do not leave this object in an inconsistant state.
669  n_ = 0;
670  throw;
671  }
672 }
673 
674 // Overridden from MatrixSymAddDelUpdateble
675 
677  value_type alpha
678  ,size_type max_size
679  )
680 {
681  // Validate input
682  if( alpha <= 0.0 ) {
683  std::ostringstream omsg;
684  omsg
685  << "MatrixSymPosDefLBFGS::initialize(alpha,max_size) : Error, "
686  << "alpha = " << alpha << " <= 0 is not allowed!";
687  throw std::invalid_argument( omsg.str() );
688  }
689  n_max_ = max_size;
690  this->init_identity(1,alpha);
691 }
692 
694  const DMatrixSliceSym &A
695  ,size_type max_size
696  ,bool force_factorization
697  ,Inertia inertia
698  ,PivotTolerances pivot_tols
699  )
700 {
701  throw std::runtime_error(
702  "MatrixSymPosDefLBFGS::initialize(A,max_size,force_refactorization,inertia) : Error, "
703  "This function is undefined for this subclass. I am so sorry for this terrible hack!" );
704 }
705 
707 {
708  return n_max_;
709 }
710 
711 MatrixSymAddDelUpdateable::Inertia MatrixSymPosDefLBFGS::inertia() const
712 {
713  return Inertia(0,0,n_);
714 }
715 
717 {
718  n_ = 0;
719 }
720 
722  const DVectorSlice *t
723  ,value_type alpha
724  ,bool force_refactorization
725  ,EEigenValType add_eigen_val
726  ,PivotTolerances pivot_tols
727  )
728 {
729  assert_initialized();
730  if( n_ == n_max_ ) {
731  std::ostringstream omsg;
732  omsg
733  << "MatrixSymPosDefLBFGS::augment_update(...) : Error, "
734  << "this->rows() = " << n_ << " == this->max_size() = " << n_max_
735  << " so we can't allow the matrix to grow!";
736  throw std::invalid_argument( omsg.str() );
737  }
738  if( t ) {
739  throw std::invalid_argument(
740  "MatrixSymPosDefLBFGS::augment_update(...) : Error, "
741  "t must be NULL in this implemention. Sorry for this hack" );
742  }
743  if( alpha <= 0.0 ) {
744  std::ostringstream omsg;
745  omsg
746  << "MatrixSymPosDefLBFGS::augment_update(...) : Error, "
747  << "alpha = " << alpha << " <= 0 is not allowed!";
748  throw std::invalid_argument( omsg.str() );
749  }
750  if( add_eigen_val == MatrixSymAddDelUpdateable::EIGEN_VAL_NEG ) {
751  std::ostringstream omsg;
752  omsg
753  << "MatrixSymPosDefLBFGS::augment_update(...) : Error, "
754  << "add_eigen_val == EIGEN_VAL_NEG is not allowed!";
755  throw std::invalid_argument( omsg.str() );
756  }
757  //
758  // Here we will do the simplest thing possible. We will just set:
759  //
760  // [ S ] -> S [ Y ] -> Y
761  // [ 0 ] [ 0 ]
762  //
763  // and let the new matrix be:
764  //
765  // [ B 0 ] -> B
766  // [ 0 1/gamma_k ]
767  //
768  // Nothing else, not even Q, needs to be updated!
769  //
770  S_.row(n_+1)(1,m_bar_) = 0.0;
771  Y_.row(n_+1)(1,m_bar_) = 0.0;
772  ++n_;
773 }
774 
776  size_type jd
777  ,bool force_refactorization
778  ,EEigenValType drop_eigen_val
779  ,PivotTolerances pivot_tols
780  )
781 {
782  assert_initialized();
783  //
784  // Removing a symmetric row and column jd is the same a removing row
785  // S(jd,:) from S and row Y(jd,:) from Y. At the same time we must
786  // update S'*Y, S'*S and Y'*Y. To see how to update these matrices
787  // not that we can represent each column of S and Y as:
788  //
789  // [ S(1:jd-1,k) ] [ Y(1:jd-1,k) ]
790  // S(:,k) = [ S(jd,k) ] , Y(:,k) = [ Y(jd,k) ] , k = 1...m_bar
791  // [ S(jd+1:n,k) ] [ Y(jd+1:n,k) ]
792  //
793  // Using the above, we can write:
794  //
795  // (S'*Y)(p,q) = S(1:jd-1,p)'*Y(1:jd-1,q) + S(jd,p)*Y(jd,q) + S(jd+1:n,p)'*Y(jd+1:n,q)
796  // , for p = 1...m_bar, q = 1...m_bar
797  //
798  // We see that the new S'*Y and the old differ by only the term S(jd,p)*Y(jd,q). Therefore, we
799  // only need to subtract off this term in each of the elements in order to update S'*Y for the
800  // deletion of this element jd. To see how to do this with BLAS, first consider subtracting
801  // of the terms by column as:
802  //
803  // (S'*Y)(:,q) <- (S'*Y)(:,q) - S(jd,:)'*Y(jd,q)
804  // , for q = 1...m_bar
805  //
806  // Then, if we put all of the columns together we get:
807  //
808  // (S'*Y)(:,:) <- (S'*Y)(:,:) - S(jd,:)'*Y(jd,:)
809  // =>
810  // (S'*Y) <- (S'*Y) - S.row(jd)*Y.row(jd)'
811  //
812  // In otherwords the above update operation is just an unsymmetric rank-1 update
813  //
814  // Similar updates for S'*S and Y'*Y are derived by just substituting matrices
815  // in to the above update for S'*Y:
816  //
817  // (S'*S) <- (S'*S) - S.row(jd)*S.row(jd)'
818  //
819  // (Y'*Y) <- (Y'*Y) - Y.row(jd)*Y.row(jd)'
820  //
821  // These updates are symmetric rank-1 updates.
822  //
823  DMatrixSlice S = this->S();
824  DMatrixSlice Y = this->Y();
825  DMatrixSlice STY = this->STY();
826  DMatrixSliceSym STS = this->STS();
827  DMatrixSliceSym YTY = this->YTY();
828  // (S'*Y) <- (S'*Y) - S.row(jd)*Y.row(jd)'
829  DenseLinAlgPack::ger( -1.0, S.row(jd), Y.row(jd), &STY );
830  // (S'*S) <- (S'*S) - S.row(jd)*S.row(jd)'
831  DenseLinAlgPack::syr( -1.0, S.row(jd), &STS );
832  // (Y'*Y) <- (Y'*Y) - Y.row(jd)*Y.row(jd)'
833  DenseLinAlgPack::syr( -1.0, Y.row(jd), &YTY );
834  // Remove row jd from S and Y one column at a time
835  // (one element at a time!)
836  if( jd < n_ ) {
837  {for( size_type k = 1; k <= m_bar_; ++k ) {
838  value_type *ptr = S.col_ptr(k);
839  std::copy( ptr + jd, ptr + n_, ptr + jd - 1 );
840  }}
841  {for( size_type k = 1; k <= m_bar_; ++k ) {
842  value_type *ptr = Y.col_ptr(k);
843  std::copy( ptr + jd, ptr + n_, ptr + jd - 1 );
844  }}
845  }
846  // Update the size
847  --n_;
848  Q_updated_ = false;
849 }
850 
851 // Private member functions
852 
853 void MatrixSymPosDefLBFGS::Vp_DtV( DVectorSlice* y, const DVectorSlice& x ) const
854 {
855  DenseLinAlgPack::Vp_MtV_assert_sizes( y->size(), m_bar_, m_bar_
856  , BLAS_Cpp::no_trans, x.size() );
857 
858  DVectorSlice::const_iterator
859  d_itr = STY_.diag(0).begin(),
860  x_itr = x.begin();
861  DVectorSlice::iterator
862  y_itr = y->begin();
863 
864  while( y_itr != y->end() )
865  *y_itr++ += (*d_itr++) * (*x_itr++);
866 }
867 
868 //
869 // We need to update the factorizations to solve for:
870 //
871 // x = inv(Q) * y => Q * x = y
872 //
873 // [ (1/gk)*S'S L ] * [ x1 ] = [ y1 ]
874 // [ L' -D ] [ x2 ] [ y2 ]
875 //
876 // We will solve the above system using a schur complement:
877 //
878 // C = (1/gk)*S'S + L*inv(D)*L'
879 //
880 // According to the referenced paper, C is p.d. so:
881 //
882 // C = J*J'
883 //
884 // We then compute the solution as:
885 //
886 // x1 = inv(C) * ( y1 + L*inv(D)*y2 )
887 // x2 = - inv(D) * ( y2 - L'*x1 )
888 //
889 // Therefore we will just update the factorization C = J*J'
890 // where the factor J is stored in QJ_.
891 //
892 
893 void MatrixSymPosDefLBFGS::update_Q() const
894 {
895  using DenseLinAlgPack::tri;
896  using DenseLinAlgPack::tri_ele;
898 
899  //
900  // We need update the factorizations to solve for:
901  //
902  // x = inv(Q) * y
903  //
904  // [ y1 ] = [ (1/gk)*S'S L ] * [ x1 ]
905  // [ y2 ] [ L' -D ] [ x2 ]
906  //
907  // We will solve the above system using the schur complement:
908  //
909  // C = (1/gk)*S'S + L*inv(D)*L'
910  //
911  // According to the referenced paper, C is p.d. so:
912  //
913  // C = J*J'
914  //
915  // We then compute the solution as:
916  //
917  // x1 = inv(C) * ( y1 + L*inv(D)*y2 )
918  // x2 = - inv(D) * ( y2 - L'*x1 )
919  //
920  // Therefore we will just update the factorization C = J*J'
921  //
922 
923  // Form the upper triangular part of C which will become J
924  // which we are using storage of QJ
925 
926  if( QJ_.rows() < m_ )
927  QJ_.resize( m_, m_ );
928 
929  const size_type
930  mb = m_bar_;
931 
932  DMatrixSlice
933  C = QJ_(1,mb,1,mb);
934 
935  // C = L * inv(D) * L'
936  //
937  // Here L is a strictly lower triangular (zero diagonal) matrix where:
938  //
939  // L = [ 0 0 ]
940  // [ Lb 0 ]
941  //
942  // Lb is lower triangular (nonzero diagonal)
943  //
944  // Therefore we compute L*inv(D)*L' as:
945  //
946  // C = [ 0 0 ] * [ Db 0 ] * [ 0 Lb' ]
947  // [ Lb 0 ] [ 0 db ] [ 0 0 ]
948  //
949  // = [ 0 0 ] = [ 0 0 ]
950  // [ 0 Cb ] [ 0 Lb*Db*Lb' ]
951  //
952  // We need to compute the upper triangular part of Cb.
953 
954  C.row(1) = 0.0;
955  if( mb > 1 )
956  comp_Cb( STY_(2,mb,1,mb-1), STY_.diag(0)(1,mb-1), &C(2,mb,2,mb) );
957 
958  // C += (1/gk)*S'S
959 
960  const DMatrixSliceSym &STS = this->STS();
961  Mp_StM( &C, (1/gamma_k_), tri( STS.gms(), STS.uplo(), BLAS_Cpp::nonunit )
962  , BLAS_Cpp::trans );
963 
964  // Now perform a cholesky factorization of C
965  // After this factorization the upper triangular part of QJ
966  // (through C) will contain the cholesky factor.
967 
968  DMatrixSliceTriEle C_upper = tri_ele( C, BLAS_Cpp::upper );
969  DenseLinAlgLAPack::potrf( &C_upper );
970 
971  Q_updated_ = true;
972 }
973 
974 void MatrixSymPosDefLBFGS::V_invQtV( DVectorSlice* x, const DVectorSlice& y ) const
975 {
976  using DenseLinAlgPack::sym;
977  using DenseLinAlgPack::tri;
980 
981  using LinAlgOpPack::Vp_V;
982  using LinAlgOpPack::V_MtV;
983 
984 
985  // Solve the system
986  //
987  // Q * x = y
988  //
989  // Using the schur complement factorization as described above.
990 
991  const size_type
992  mb = m_bar_;
993 
994  if(!Q_updated_) {
995  update_Q();
996  }
997 
998  DVectorSlice
999  x1 = (*x)(1,mb),
1000  x2 = (*x)(mb+1,2*mb);
1001 
1002  const DVectorSlice
1003  y1 = y(1,mb),
1004  y2 = y(mb+1,2*mb);
1005 
1006  // //////////////////////////////////////
1007  // x1 = inv(C) * ( y1 + L*inv(D)*y2 )
1008  // = inv(J'*J) * r
1009  // = inv(J) * inv(J') * r
1010 
1011  { // x1 = inv(D) * y2
1012  DVectorSlice::const_iterator
1013  d_itr = STY_.diag(0).begin(),
1014  y2_itr = y2.begin();
1015  DVectorSlice::iterator
1016  x1_itr = x1.begin();
1017  while( x1_itr != x1.end() )
1018  *x1_itr++ = *y2_itr++ / *d_itr++;
1019  }
1020 
1021  // x1 = L * x1
1022  //
1023  // = [ 0 0 ] * [ x1(1:mb-1) ]
1024  // [ Lb 0 ] [ x1(mb) ]
1025  //
1026  // = [ 0 ]
1027  // [ Lb*x1(1:mb-1) ]
1028  //
1029  if( mb > 1 ) {
1030  // x1(2,mb) = x1(1,mb-1) ( copy from mb-1 to mb then mb-2 to mb-1
1031  // etc. so that we don't overwrite the elements we need to copy ).
1032  DVectorSlice
1033  x1a = x1(1,mb-1),
1034  x1b = x1(2,mb);
1035  std::copy( x1a.rbegin(), x1a.rend(), x1b.rbegin() );
1036  V_MtV( &x1b, Lb(), BLAS_Cpp::no_trans, x1b );
1037  }
1038  x1(1) = 0.0;
1039 
1040  // r = x1 += y1
1041  Vp_V( &x1, y1 );
1042 
1043  // x1 = inv(J') * r
1044  const DMatrixSliceTri J = tri( QJ_(1,mb,1,mb), BLAS_Cpp::upper, BLAS_Cpp::nonunit );
1045  V_InvMtV( &x1, J, BLAS_Cpp::trans, x1 );
1046 
1047  // x1 = inv(J) * x1
1048  V_InvMtV( &x1, J, BLAS_Cpp::no_trans, x1 );
1049 
1050  // /////////////////////////////////////
1051  // x2 = inv(D) * ( - y2 + L'*x1 )
1052 
1053  // x2 = L'*x1
1054  //
1055  // = [ 0 Lb' ] * [ x1(1) ]
1056  // [ 0 0 ] [ x1(2,mb) ]
1057  //
1058  // = [ Lb'*x1(2,mb) ]
1059  // [ 0 ]
1060  //
1061  if( mb > 1 ) {
1062  V_MtV( &x2(1,mb-1), Lb(), BLAS_Cpp::trans, x1(2,mb) );
1063  }
1064  x2(mb) = 0.0;
1065 
1066  // x2 += -y2
1067  Vp_StV( &x2, -1.0, y2 );
1068 
1069  // x2 = inv(D) * x2
1070  {
1071  DVectorSlice::const_iterator
1072  d_itr = STY_.diag(0).begin();
1073  DVectorSlice::iterator
1074  x2_itr = x2.begin();
1075  while( x2_itr != x2.end() )
1076  *x2_itr++ /= *d_itr++;
1077  }
1078 }
1079 
1080 void MatrixSymPosDefLBFGS::assert_initialized() const
1081 {
1082  if(!n_)
1083  throw std::logic_error( "MatrixSymPosDefLBFGS::assert_initialized() : "
1084  "Error, matrix not initialized" );
1085 }
1086 
1087 } // end namespace ConstrainedOptPack
1088 
1089 namespace {
1090 
1091 void comp_Cb( const DMatrixSlice& Lb, const DVectorSlice& Db_diag
1092  , DMatrixSlice* Cb )
1093 {
1094  // Lb * inv(Db) * Lb =
1095  //
1096  // [ l(1,1) ] [ dd(1) ] [ l(1,1) l(2,1) ... l(p,1) ]
1097  // [ l(2,1) l(2,2) ] [ dd(2) ] [ l(2,2) ... l(p,2) ]
1098  // [ . . . ] * [ . ] * [ . . ]
1099  // [ l(p,1) l(p,2) ... l(p,p) ] [ dd(p) ] [ l(p,p) ]
1100  //
1101  //
1102  // [ l(1,1)*dd(1)*l(1,1) l(1,1)*dd(1)*l(2,1) ... l(1,1)*dd(1)*l(p,1) ]
1103  // [ 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) ]
1104  // [ . . ...
1105  // [ symmetric symmetric ... sum( l(p,i)*dd(i)*l(p,i), i=1,..,p ) ]
1106  //
1107  // Therefore we can express the upper triangular elemetns of Cb as:
1108  //
1109  // Cb(i,j) = sum( l(i,k)*dd(k)*l(j,k), k = 1,..,i )
1110 
1112  typedef DenseLinAlgPack::value_type value_type;
1113 
1114  TEUCHOS_TEST_FOR_EXCEPT( !( Lb.rows() == Cb->rows() && Cb->rows() == Db_diag.size() ) );
1115 
1116  const size_type p = Db_diag.size();
1117 
1118  for( size_type i = 1; i <= p; ++i ) {
1119  for( size_type j = i; j <= p; ++j ) {
1120  value_type &c = (*Cb)(i,j) = 0.0;
1121  for( size_type k = 1; k <= i; ++k ) {
1122  c += Lb(i,k) * Lb(j,k) / Db_diag(k);
1123  }
1124  }
1125  }
1126 
1127  // ToDo: Make the above operation more efficent if needed!
1128 }
1129 
1130 } // end namespace
1131 
1132 #endif // 0
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)
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 initialize(value_type alpha, size_type max_size)
This is fine as long as alpha > 0.0.
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 delete_update(size_type jd, bool force_refactorization, EEigenValType drop_eigen_val, PivotTolerances pivot_tols)
Should always succeed unless user gives a wrong value for drop_eigen_val.
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)
Inertia inertia() const
Returns (0,0,rows())
void Vp_StMtV(VectorMutable *v_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs1, const Vector &v_rhs2, value_type beta) const
Transp
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
void augment_update(const DVectorSlice *t, value_type alpha, bool force_refactorization, EEigenValType add_eigen_val, PivotTolerances pivot_tols)
Augment the matrix to add a row and column.
void set_uninitialized()
Will set rows() == 0.
void Vp_V(VectorMutable *v_lhs, const V &V_rhs)