# ConstrainedOptPack_MatrixSymPosDefLBFGSSerial.cpp

```00001 // @HEADER
00002 // ***********************************************************************
00003 //
00004 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
00005 //                  Copyright (2003) Sandia Corporation
00006 //
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 //
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov)
00025 //
00026 // ***********************************************************************
00027 // @HEADER
00028 //
00029 // Let's define a compact representation for the matrix B^{k} and
00030 // its inverse H^{k} = inv(B^{k}).
00031 //
00032 // Bk = (1/gk)*I - [ (1/gk)*S  Y ] * inv( [ (1/gk)*S'S   L ]     [ (1/gk)*S' ]
00033 //                                        [    L'       -D ] ) * [    Y'     ]
00034 //                                        \________________/
00035 //                                                Q
00036 //
00037 // Hk = gk*I + [ S  gk*Y ] * [ inv(R')*(D+gk*Y'Y)*inv(R)     -inv(R') ] * [   S'  ]
00038 //                           [            -inv(R)                0    ]   [ gk*Y' ]
00039 //
00040 // where:
00041 //
00042 // gk = gamma_k <: R
00043 //
00044 //      [ s^{1}'*y^{1}  s^{1}'*y^{2}  ...   s^{1}'*y^{m}  ]
00045 //  S'Y = [ s^{2}'*y^{1}  s^{2}'*y^{2}  ...   s^{2}'*y^{m}  ] <: R^(m x m)
00046 //      [ .       .           .       ]
00047 //      [ s^{m}'*y^{1}  s^{m}'*y^{2}  ...   s^{m}'*y^{m}  ]
00048 //
00049 //      [ s^{1}'*y^{1}  0       ...   0       ]
00050 //  D =   [ 0       s^{2}'*y^{2}  ...   0       ] <: R^(m x m)
00051 //      [ .       .           .       ]
00052 //      [ 0       0       ...   s^{m}'*y^{m}  ]
00053 //
00054 //  R = upper triangular part of S'Y
00055 //
00056 //  L = lower tirangular part of S'Y with zeros on the diagonal
00057 //
00058
00059 #include <assert.h>
00060
00061 #include "ConstrainedOptPack_MatrixSymPosDefLBFGS.hpp"
00062 #include "ConstrainedOptPack/src/AbstractLinAlgPack_BFGS_helpers.hpp"
00063 #include "DenseLinAlgPack_LinAlgOpPack.hpp"
00064 #include "DenseLinAlgPack_DMatrixOut.hpp"
00065 #include "DenseLinAlgLAPack.hpp"
00066
00067 namespace {
00068
00069   using DenseLinAlgPack::DVectorSlice;
00070   using DenseLinAlgPack::DMatrixSlice;
00071
00079   void comp_Cb( const DMatrixSlice& Lb, const DVectorSlice& Db_diag
00080     , DMatrixSlice* Cb );
00081
00082 } // end namespace
00083
00084 namespace ConstrainedOptPack {
00085
00086 // /////////////////////////////////
00087 // Inline private member functions
00088
00089 inline
00090 const DMatrixSliceTri MatrixSymPosDefLBFGS::R() const
00091 {
00092   return DenseLinAlgPack::tri( STY_(1,m_bar_,1,m_bar_), BLAS_Cpp::upper, BLAS_Cpp::nonunit );
00093 }
00094
00095 inline
00096 const DMatrixSliceTri MatrixSymPosDefLBFGS::Lb() const
00097 {
00098   return DenseLinAlgPack::tri( STY_(2,m_bar_,1,m_bar_-1), BLAS_Cpp::lower, BLAS_Cpp::nonunit );
00099 }
00100
00101 inline
00102 DMatrixSlice MatrixSymPosDefLBFGS::STY()
00103 {
00104   return STY_(1,m_bar_,1,m_bar_);
00105 }
00106
00107 inline
00108 const DMatrixSlice MatrixSymPosDefLBFGS::STY() const
00109 {
00110   return STY_(1,m_bar_,1,m_bar_);
00111 }
00112
00113 inline
00114 DMatrixSliceSym MatrixSymPosDefLBFGS::STS()
00115 {
00116   return DenseLinAlgPack::nonconst_sym( STSYTY_(2,m_bar_+1,1,m_bar_),BLAS_Cpp::lower );
00117 }
00118
00119 inline
00120 const DMatrixSliceSym MatrixSymPosDefLBFGS::STS() const
00121 {
00122   return DenseLinAlgPack::sym( STSYTY_(2,m_bar_+1,1,m_bar_),BLAS_Cpp::lower );
00123 }
00124
00125 inline
00126 DMatrixSliceSym MatrixSymPosDefLBFGS::YTY()
00127 {
00128   return DenseLinAlgPack::nonconst_sym( STSYTY_(1,m_bar_,2,m_bar_+1),BLAS_Cpp::upper );
00129 }
00130
00131 inline
00132 const DMatrixSliceSym MatrixSymPosDefLBFGS::YTY() const
00133 {
00134   return DenseLinAlgPack::sym( STSYTY_(1,m_bar_,2,m_bar_+1),BLAS_Cpp::upper );
00135 }
00136
00137 // ///////////////////////
00138 // Nonlinined functions
00139
00140 MatrixSymPosDefLBFGS::MatrixSymPosDefLBFGS(
00141   size_type   max_size
00142     ,size_type  m
00143   ,bool       maintain_original
00144   ,bool       maintain_inverse
00145   ,bool       auto_rescaling
00146   )
00147 {
00148   initial_setup(max_size,m,maintain_original,maintain_inverse,auto_rescaling);
00149 }
00150
00151 void MatrixSymPosDefLBFGS::initial_setup(
00152   size_type   max_size
00153     ,size_type  m
00154   ,bool       maintain_original
00155   ,bool       maintain_inverse
00156   ,bool       auto_rescaling
00157   )
00158 {
00159   // Validate input
00160   if( !maintain_original && !maintain_inverse )
00161     throw std::invalid_argument(
00162       "MatrixSymPosDefLBFGS::initial_setup(...) : "
00163       "Error, both maintain_original and maintain_inverse can not both be false!" );
00164   if( m < 1 )
00165     throw std::invalid_argument(
00166       "MatrixSymPosDefLBFGS::set_num_updates_stored(m) : "
00167       "Error, the number of storage locations must be > 0" );
00168   maintain_original_ = maintain_original;
00169   maintain_inverse_  = maintain_inverse;
00170   m_                 = m;
00171   n_                 = 0; // make uninitialized
00172   n_max_             = max_size;
00173   num_secant_updates_= 0;
00174 }
00175
00176 // Overridden from Matrix
00177
00178 size_type MatrixSymPosDefLBFGS::rows() const
00179 {
00180   return n_;
00181 }
00182
00183 // Overridden from MatrixOp
00184
00185 std::ostream& MatrixSymPosDefLBFGS::output(std::ostream& out) const
00186 {
00187   assert_initialized();
00188   out << "*** Limited Memory BFGS matrix.\n"
00189     << "Conversion to dense =\n";
00190   MatrixOp::output(out);
00191   out << "\n*** Stored quantities\n"
00192     << "\ngamma_k = " << gamma_k_ << std::endl;
00193   if( m_bar_ ) {
00194     out << "\nS =\n" << S()
00195       << "\nY =\n" << Y()
00196       << "\nS'Y =\n" << STY_(1,m_bar_,1,m_bar_)
00197       << "\nlower(S'S) \\ zero diagonal \\ upper(Y'Y) =\n"
00198         << STSYTY_(1,m_bar_+1,1,m_bar_+1)
00199       << "\nQ updated? = " << Q_updated_ << std::endl
00200       << "\nCholesky of schur complement of Q, QJ =\n" << QJ_(1,m_bar_,1,m_bar_);
00201   }
00202   return out;
00203 }
00204
00205 MatrixOp& MatrixSymPosDefLBFGS::operator=(const MatrixOp& m)
00206 {
00207   const MatrixSymPosDefLBFGS *p_m = dynamic_cast<const MatrixSymPosDefLBFGS*>(&m);
00208   if(p_m) {
00209     if( p_m == this ) return *this; // assignment to self
00210     // Important: Assign all members here.
00211     auto_rescaling_      = p_m->auto_rescaling_;
00212     maintain_original_   = p_m->maintain_original_;
00213     original_is_updated_ = p_m->original_is_updated_;
00214     maintain_inverse_    = p_m->maintain_inverse_;
00215     inverse_is_updated_  = p_m->inverse_is_updated_;
00216     n_max_             = p_m->n_max_;
00217     n_               = p_m->n_;
00218     m_               = p_m->m_;
00219     m_bar_             = p_m->m_bar_;
00220     k_bar_             = p_m->k_bar_;
00221     num_secant_updates_  = p_m->num_secant_updates_;
00222     gamma_k_           = p_m->gamma_k_;
00223     S_               = p_m->S_;
00224     Y_               = p_m->Y_;
00225     STY_             = p_m->STY_;
00226     STSYTY_            = p_m->STSYTY_;
00227     Q_updated_           = p_m->Q_updated_;
00228     QJ_              = p_m->QJ_;
00229   }
00230   else {
00231     throw std::invalid_argument("MatrixSymPosDefLBFGS::operator=(const MatrixOp& m)"
00232       " : The concrete type of m is not a subclass of MatrixSymPosDefLBFGS as expected" );
00233   }
00234   return *this;
00235 }
00236
00237 // Level-2 BLAS
00238
00239 void MatrixSymPosDefLBFGS::Vp_StMtV(
00240     DVectorSlice* y, value_type alpha, BLAS_Cpp::Transp trans_rhs1
00241   , const DVectorSlice& x, value_type beta) const
00242 {
00243   using LinAlgOpPack::V_StMtV;
00244   using LinAlgOpPack::V_MtV;
00245
00246   using DenseLinAlgPack::Vt_S;
00247   using DenseLinAlgPack::Vp_StV;
00248   using DenseLinAlgPack::Vp_StMtV;
00249
00250   assert_initialized();
00251
00252   assert( original_is_updated_ ); // For now just always update
00253
00254   // y = b*y + Bk * x
00255   //
00256   // y = b*y + (1/gk)*x - [ (1/gk)*S  Y ] * inv(Q) * [ (1/gk)*S' ] * x
00257   //                                                 [     Y'    ]
00258   // Perform the following operations (in order):
00259   //
00260   // y = b*y
00261   //
00262   // y += (1/gk)*x
00263   //
00264   // t1 = [ (1/gk)*S'*x ]   <: R^(2*m)
00265   //    [      Y'*x   ]
00266   //
00267   // t2 = inv(Q) * t1     <: R^(2*m)
00268   //
00269   // y += -(1/gk) * S * t2(1:m)
00270   //
00271   // y += -1.0 * Y * t2(m+1,2m)
00272
00273   const value_type
00274     invgk = 1.0 / gamma_k_;
00275
00276   // y = b*y
00277
00278   if( beta == 0.0 )
00279     *y = beta;
00280   else
00281     Vt_S( y, beta );
00282
00283   // y += (1/gk)*x
00284
00285   Vp_StV( y, invgk, x );
00286
00287   if( !m_bar_ )
00288     return; // No updates have been added yet.
00289
00290   // Get workspace
00291
00292   if( work_.size() < 4 * m_ )
00293     work_.resize( 4 * m_ );
00294
00295   const size_type
00296     mb = m_bar_;
00297
00298   const size_type
00299     t1s = 1,
00300     t1n = 2*mb,
00301     t2s = t1s+t1n,
00302     t2n = 2*mb;
00303
00304   DVectorSlice
00305     t1 = work_( t1s,  t1s + t1n - 1 ),
00306     t2 = work_( t2s,  t2s + t2n - 1   );
00307
00308   const DMatrixSlice
00309     &S = this->S(),
00310     &Y = this->Y();
00311
00312   // t1 = [ (1/gk)*S'*x ]
00313   //    [      Y'*x   ]
00314
00315   V_StMtV( &t1(1,mb), invgk, S, BLAS_Cpp::trans, x );
00316   V_MtV( &t1(mb+1,2*mb), Y, BLAS_Cpp::trans, x );
00317
00318   // t2 = inv(Q) * t1
00319
00320   V_invQtV( &t2, t1 );
00321
00322   // y += -(1/gk) * S * t2(1:m)
00323
00324   Vp_StMtV( y, -invgk, S, BLAS_Cpp::no_trans, t2(1,mb) );
00325
00326   // y += -1.0 * Y * t2(m+1,2m)
00327
00328   Vp_StMtV( y, -1.0, Y, BLAS_Cpp::no_trans, t2(mb+1,2*mb) );
00329
00330 }
00331
00332 // Overridden from MatrixWithOpFactorized
00333
00334 // Level-2 BLAS
00335
00336 void MatrixSymPosDefLBFGS::V_InvMtV( DVectorSlice* y, BLAS_Cpp::Transp trans_rhs1
00337   , const DVectorSlice& x ) const
00338 {
00339   using DenseLinAlgPack::V_mV;
00340   using DenseLinAlgPack::V_StV;
00341   using DenseLinAlgPack::V_InvMtV;
00342
00343   using LinAlgOpPack::Vp_V;
00344   using LinAlgOpPack::V_MtV;
00345   using LinAlgOpPack::V_StMtV;
00346   using LinAlgOpPack::Vp_MtV;
00347   using LinAlgOpPack::Vp_StMtV;
00348
00349   assert_initialized();
00350
00351   assert( inverse_is_updated_ ); // For now just always update
00352
00353   // y = inv(Bk) * x = Hk * x
00354   //
00355   // = gk*x + [S gk*Y] * [ inv(R')*(D+gk*Y'Y)*inv(R)     -inv(R') ] * [   S'  ] * x
00356   //                     [            -inv(R)                0    ]   [ gk*Y' ]
00357   //
00358   // Perform in the following (in order):
00359   //
00360   // y = gk*x
00361   //
00362   // t1 = [   S'*x  ]         <: R^(2*m)
00363   //      [ gk*Y'*x ]
00364   //
00365   // t2 = inv(R) * t1(1:m)      <: R^(m)
00366   //
00367   // t3 = - inv(R') * t1(m+1,2*m)   <: R^(m)
00368   //
00369   // t4 = gk * Y'Y * t2       <: R^(m)
00370   //
00371   // t4 += D*t2
00372   //
00373   // t5 = inv(R') * t4        <: R^(m)
00374   //
00375   // t5 += t3
00376   //
00377   // y += S*t5
00378   //
00379   // y += -gk*Y*t2
00380
00381   // y = gk*x
00382   V_StV( y, gamma_k_, x );
00383
00384   const size_type
00385     mb = m_bar_;
00386
00387   if( !mb )
00388     return; // No updates have been performed.
00389
00390   // Get workspace
00391
00392   if( work_.size() < 6*m_ )
00393     work_.resize( 6*m_ );
00394
00395   const size_type
00396     t1s   = 1,
00397     t1n   = 2*mb,
00398     t2s   = t1s + t1n,
00399     t2n   = mb,
00400     t3s   = t2s + t2n,
00401     t3n   = mb,
00402     t4s   = t3s + t3n,
00403     t4n   = mb,
00404     t5s   = t4s + t4n,
00405     t5n   = mb;
00406
00407   DVectorSlice
00408     t1  = work_( t1s, t1s + t1n - 1 ),
00409     t2  = work_( t2s, t2s + t2n - 1 ),
00410     t3  = work_( t3s, t3s + t3n - 1 ),
00411     t4  = work_( t4s, t4s + t4n - 1 ),
00412     t5  = work_( t5s, t5s + t5n - 1 );
00413
00414   const DMatrixSlice
00415     &S = this->S(),
00416     &Y = this->Y();
00417
00418   const DMatrixSliceTri
00419     &R = this->R();
00420
00421   const DMatrixSliceSym
00422     &YTY = this->YTY();
00423
00424   // t1 = [   S'*x  ]
00425   //      [ gk*Y'*x ]
00426   V_MtV( &t1(1,mb), S, BLAS_Cpp::trans, x );
00427   V_StMtV( &t1(mb+1,2*mb), gamma_k_, Y, BLAS_Cpp::trans, x );
00428
00429   // t2 = inv(R) * t1(1:m)
00430   V_InvMtV( &t2, R, BLAS_Cpp::no_trans, t1(1,mb) );
00431
00432   // t3 = - inv(R') * t1(m+1,2*m)
00433   V_mV( &t3, t1(mb+1,2*mb) );
00434   V_InvMtV( &t3, R, BLAS_Cpp::trans, t3 );
00435
00436   // t4 = gk * Y'Y * t2
00437   V_StMtV( &t4, gamma_k_, YTY, BLAS_Cpp::no_trans, t2 );
00438
00439   // t4 += D*t2
00440   Vp_DtV( &t4, t2 );
00441
00442   // t5 = inv(R') * t4
00443   V_InvMtV( &t5, R, BLAS_Cpp::trans, t4 );
00444
00445   // t5 += t3
00446   Vp_V( &t5, t3 );
00447
00448   // y += S*t5
00449   Vp_MtV( y, S, BLAS_Cpp::no_trans, t5 );
00450
00451   // y += -gk*Y*t2
00452   Vp_StMtV( y, -gamma_k_, Y, BLAS_Cpp::no_trans, t2 );
00453
00454 }
00455
00456 // Overridden from MatrixSymSecant
00457
00458 void MatrixSymPosDefLBFGS::init_identity(size_type n, value_type alpha)
00459 {
00460   // Validate input
00461   if( alpha <= 0.0 ) {
00462     std::ostringstream omsg;
00463     omsg
00464       << "MatrixSymPosDefLBFGS::init_identity(n,alpha) : Error, "
00465       << "alpha = " << alpha << " <= 0 is not allowed!";
00466     throw std::invalid_argument( omsg.str() );
00467   }
00468   if( n_max_ == 0 ) {
00469     n_max_ = n;
00470   }
00471   else if( n > n_max_ ) {
00472     std::ostringstream omsg;
00473     omsg
00474       << "MatrixSymPosDefLBFGS::init_identity(n,alpha) : Error, "
00475       << "n = " << n << " > max_size = " << n_max_;
00476     throw std::invalid_argument( omsg.str() );
00477   }
00478
00479   // Resize storage
00480   S_.resize( n_max_, m_ );
00481   Y_.resize( n_max_, m_ );
00482   STY_.resize( m_, m_ );
00483   STSYTY_.resize( m_+1, m_+1 );
00484   STSYTY_.diag(0) = 0.0;
00485
00486   gamma_k_ = 1.0/alpha;
00487
00488   // Initialize counters
00489   k_bar_  = 0;
00490   m_bar_  = 0;
00491
00492   n_ = n;  // initialized;
00493   original_is_updated_ = true; // This will never change for now
00494   inverse_is_updated_  = true; // This will never change for now
00495   num_secant_updates_  = 0;
00496 }
00497
00498 void MatrixSymPosDefLBFGS::init_diagonal(const DVectorSlice& diag)
00499 {
00500   using DenseLinAlgPack::norm_inf;
00501   init_identity( diag.size(), norm_inf(diag) );
00502 }
00503
00504 void MatrixSymPosDefLBFGS::secant_update(
00505   DVectorSlice* s, DVectorSlice* y, DVectorSlice* Bs)
00506 {
00507   using DenseLinAlgPack::dot;
00508   using DenseLinAlgPack::norm_2;
00509
00510   using LinAlgOpPack::V_MtV;
00511
00512   assert_initialized();
00513
00514   // Check skipping the BFGS update
00515   const value_type
00516     sTy       = dot(*s,*y);
00517   std::ostringstream omsg;
00518   if( !BFGS_sTy_suff_p_d(*s,*y,&sTy,&omsg,"MatrixSymPosDefLBFGS::secant_update(...)" ) ) {
00519     throw UpdateSkippedException( omsg.str() );
00520   }
00521
00522   try {
00523
00524   // Update counters
00525   if( k_bar_ == m_ ) {
00526 //    // We are at the end storage so loop back around again
00527 //    k_bar_ = 1;
00528     // We are at the end of the storage so remove the oldest stored update
00529     // and move updates to make room for the new update.  This has to be done for the
00530     // the matrix to behave properly
00531     {for( size_type k = 1; k <= m_-1; ++k ) {
00532       S_.col(k) = S_.col(k+1);                              // Shift S.col() to the left
00533       Y_.col(k) = Y_.col(k+1);                              // Shift Y.col() to the left
00534       STY_.col(k)(1,m_-1) = STY_.col(k+1)(2,m_);            // Move submatrix STY(2,m-1,2,m-1) up and left
00535       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
00536       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
00537     }}
00538   }
00539   else {
00540     k_bar_++;
00541   }
00542   if( m_bar_ < m_ ) {
00543     // This is the first few updates where we have not maxed out the storage.
00544     m_bar_++;
00545   }
00546
00547   // Set the update vectors
00548   S_.col(k_bar_)(1,n_) = *s;
00549   Y_.col(k_bar_)(1,n_) = *y;
00550
00551   // /////////////////////////////////////////////////////////////////////////////////////
00552   // Update S'Y
00553   //
00554   // Update the row and column k_bar
00555   //
00556   //  S'Y =
00557   //
00558   //  [ s(1)'*y(1)    ...   s(1)'*y(k_bar)    ...   s(1)'*y(m_bar)    ]
00559   //  [ .           .             .         ] row
00560   //  [ s(k_bar)'*y(1)  ...   s(k_bar)'*y(k_bar)  ...   s(k_bar)'*y(m_bar)  ] k_bar
00561   //  [ .           .             .         ]
00562   //  [ s(m_bar)'*y(1)  ...   s(m_bar)'*y(k_bar)  ...   s(m_bar)'*y(m_bar)  ]
00563   //
00564   //                col k_bar
00565   //
00566   // Therefore we set:
00567   //  (S'Y)(:,k_bar) =  S'*y(k_bar)
00568   //  (S'Y)(k_bar,:) =  s(k_bar)'*Y
00569
00570   const DMatrixSlice
00571     &S = this->S(),
00572     &Y = this->Y();
00573
00574   //  (S'Y)(:,k_bar) =  S'*y(k_bar)
00575   V_MtV( &STY_.col(k_bar_)(1,m_bar_), S, BLAS_Cpp::trans, Y.col(k_bar_) );
00576
00577   //  (S'Y)(k_bar,:)' =  Y'*s(k_bar)
00578   V_MtV( &STY_.row(k_bar_)(1,m_bar_), Y, BLAS_Cpp::trans, S.col(k_bar_) );
00579
00580   // /////////////////////////////////////////////////////////////////
00581   // Update S'S
00582   //
00583   //  S'S =
00584   //
00585   //  [ s(1)'*s(1)    ...   symmetric         symmetric     ]
00586   //  [ .           .             .         ] row
00587   //  [ s(k_bar)'*s(1)  ...   s(k_bar)'*s(k_bar)  ...   symmetric     ] k_bar
00588   //  [ .           .             .         ]
00589   //  [ s(m_bar)'*s(1)  ...   s(m_bar)'*s(k_bar)  ...   s(m_bar)'*s(m_bar)  ]
00590   //
00591   //                col k_bar
00592   //
00593   // Here we will update the lower triangular part of S'S.  To do this we
00594   // only need to compute:
00595   //    t = S'*s(k_bar) = { s(k_bar)' * [ s(1),..,s(k_bar),..,s(m_bar) ]  }'
00596   // then set the appropriate rows and columns of S'S.
00597
00598   if( work_.size() < m_ )
00599     work_.resize(m_);
00600
00601   // work = S'*s(k_bar)
00602   V_MtV( &work_(1,m_bar_), S, BLAS_Cpp::trans, S.col(k_bar_) );
00603
00604   // Set row elements
00605   STSYTY_.row(k_bar_+1)(1,k_bar_) = work_(1,k_bar_);
00606   // Set column elements
00607   STSYTY_.col(k_bar_)(k_bar_+1,m_bar_+1) = work_(k_bar_,m_bar_);
00608
00609   // /////////////////////////////////////////////////////////////////////////////////////
00610   // Update Y'Y
00611   //
00612   // Update the row and column k_bar
00613   //
00614   //  Y'Y =
00615   //
00616   //  [ y(1)'*y(1)    ...   y(1)'*y(k_bar)    ...   y(1)'*y(m_bar)    ]
00617   //  [ .           .             .         ] row
00618   //  [ symmetric   ...   y(k_bar)'*y(k_bar)  ...   y(k_bar)'*y(m_bar)  ] k_bar
00619   //  [ .           .             .         ]
00620   //  [ symmetric   ...   symmetric     ...   y(m_bar)'*y(m_bar)  ]
00621   //
00622   //                col k_bar
00623   //
00624   // Here we will update the upper triangular part of Y'Y.  To do this we
00625   // only need to compute:
00626   //    t = Y'*y(k_bar) = { y(k_bar)' * [ y(1),..,y(k_bar),..,y(m_bar) ]  }'
00627   // then set the appropriate rows and columns of Y'Y.
00628
00629   // work = Y'*y(k_bar)
00630   V_MtV( &work_(1,m_bar_), Y, BLAS_Cpp::trans, Y.col(k_bar_) );
00631
00632   // Set row elements
00633   STSYTY_.col(k_bar_+1)(1,k_bar_) = work_(1,k_bar_);
00634   // Set column elements
00635   STSYTY_.row(k_bar_)(k_bar_+1,m_bar_+1) = work_(k_bar_,m_bar_);
00636
00637   // /////////////////////////////
00638   // Update gamma_k
00639
00640   // gamma_k = s'*y / y'*y
00641   if(auto_rescaling_)
00642     gamma_k_ = STY_(k_bar_,k_bar_) / STSYTY_(k_bar_,k_bar_+1);
00643
00644   // We do not initially update Q unless we have to form a matrix-vector
00645   // product later.
00646
00647   Q_updated_ = false;
00648   num_secant_updates_++;
00649
00650   } //  end try
00651   catch(...) {
00652     // If we throw any exception the we should make the matrix uninitialized
00653     // so that we do not leave this object in an inconsistant state.
00654     n_ = 0;
00655     throw;
00656   }
00657 }
00658
00659 // Overridden from MatrixSymAddDelUpdateble
00660
00661 void MatrixSymPosDefLBFGS::initialize(
00662   value_type         alpha
00663   ,size_type         max_size
00664   )
00665 {
00666   // Validate input
00667   if( alpha <= 0.0 ) {
00668     std::ostringstream omsg;
00669     omsg
00670       << "MatrixSymPosDefLBFGS::initialize(alpha,max_size) : Error, "
00671       << "alpha = " << alpha << " <= 0 is not allowed!";
00672     throw std::invalid_argument( omsg.str() );
00673   }
00674   n_max_ = max_size;
00675   this->init_identity(1,alpha);
00676 }
00677
00678 void MatrixSymPosDefLBFGS::initialize(
00679   const DMatrixSliceSym      &A
00680   ,size_type         max_size
00681   ,bool              force_factorization
00682   ,Inertia           inertia
00683   ,PivotTolerances   pivot_tols
00684   )
00685 {
00686   throw std::runtime_error(
00687     "MatrixSymPosDefLBFGS::initialize(A,max_size,force_refactorization,inertia) : Error, "
00688     "This function is undefined for this subclass.  I am so sorry for this terrible hack!" );
00689 }
00690
00691 size_type MatrixSymPosDefLBFGS::max_size() const
00692 {
00693   return n_max_;
00694 }
00695
00696 MatrixSymAddDelUpdateable::Inertia MatrixSymPosDefLBFGS::inertia() const
00697 {
00698   return Inertia(0,0,n_);
00699 }
00700
00701 void MatrixSymPosDefLBFGS::set_uninitialized()
00702 {
00703   n_ = 0;
00704 }
00705
00706 void MatrixSymPosDefLBFGS::augment_update(
00707   const DVectorSlice  *t
00708   ,value_type        alpha
00709   ,bool              force_refactorization
00710   ,EEigenValType     add_eigen_val
00711   ,PivotTolerances   pivot_tols
00712   )
00713 {
00714   assert_initialized();
00715   if( n_ == n_max_ ) {
00716     std::ostringstream omsg;
00717     omsg
00718       << "MatrixSymPosDefLBFGS::augment_update(...) : Error, "
00719       << "this->rows() = " << n_ << " == this->max_size() = " << n_max_
00720       << " so we can't allow the matrix to grow!";
00721     throw std::invalid_argument( omsg.str() );
00722   }
00723   if( t ) {
00724     throw std::invalid_argument(
00725       "MatrixSymPosDefLBFGS::augment_update(...) : Error, "
00726       "t must be NULL in this implemention.  Sorry for this hack" );
00727   }
00728   if( alpha <= 0.0 ) {
00729     std::ostringstream omsg;
00730     omsg
00731       << "MatrixSymPosDefLBFGS::augment_update(...) : Error, "
00732       << "alpha = " << alpha << " <= 0 is not allowed!";
00733     throw std::invalid_argument( omsg.str() );
00734   }
00735   if( add_eigen_val == MatrixSymAddDelUpdateable::EIGEN_VAL_NEG ) {
00736     std::ostringstream omsg;
00737     omsg
00738       << "MatrixSymPosDefLBFGS::augment_update(...) : Error, "
00739       << "add_eigen_val == EIGEN_VAL_NEG is not allowed!";
00740     throw std::invalid_argument( omsg.str() );
00741   }
00742   //
00743   // Here we will do the simplest thing possible.  We will just  set:
00744   //
00745   // [ S ] -> S       [ Y ] -> Y
00746   // [ 0 ]            [ 0 ]
00747   //
00748   // and let the new matrix be:
00749   //
00750   // [ B      0     ] -> B
00751   // [ 0  1/gamma_k ]
00752   //
00753   // Nothing else, not even Q, needs to be updated!
00754   //
00755   S_.row(n_+1)(1,m_bar_) = 0.0;
00756   Y_.row(n_+1)(1,m_bar_) = 0.0;
00757   ++n_;
00758 }
00759
00760 void MatrixSymPosDefLBFGS::delete_update(
00761   size_type          jd
00762   ,bool              force_refactorization
00763   ,EEigenValType     drop_eigen_val
00764   ,PivotTolerances   pivot_tols
00765   )
00766 {
00767   assert_initialized();
00768   //
00769   // Removing a symmetric row and column jd is the same a removing row
00770   // S(jd,:) from S and row Y(jd,:) from Y.  At the same time we must
00771   // update S'*Y, S'*S and Y'*Y.  To see how to update these matrices
00772   // not that we can represent each column of S and Y as:
00773   //
00774   //           [ S(1:jd-1,k) ]                 [ Y(1:jd-1,k) ]
00775   // S(:,k) =  [ S(jd,k)     ]     , Y(:,k) =  [ Y(jd,k)     ]  , k = 1...m_bar
00776   //           [ S(jd+1:n,k) ]                 [ Y(jd+1:n,k) ]
00777   //
00778   // Using the above, we can write:
00779   //
00780   // (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)
00781   //     , for p = 1...m_bar, q = 1...m_bar
00782   //
00783   // We see that the new S'*Y and the old differ by only the term S(jd,p)*Y(jd,q).  Therefore, we
00784   // only need to subtract off this term in each of the elements in order to update S'*Y for the
00785   // deletion of this element jd.  To see how to do this with BLAS, first consider subtracting
00786   // of the terms by column as:
00787   //
00788   // (S'*Y)(:,q) <- (S'*Y)(:,q) - S(jd,:)'*Y(jd,q)
00789   //     , for q = 1...m_bar
00790   //
00791   // Then, if we put all of the columns together we get:
00792   //
00793   // (S'*Y)(:,:) <- (S'*Y)(:,:) - S(jd,:)'*Y(jd,:)
00794   // =>
00795   // (S'*Y) <- (S'*Y) - S.row(jd)*Y.row(jd)'
00796   //
00797   // In otherwords the above update operation is just an unsymmetric rank-1 update
00798   //
00799   // Similar updates for S'*S and Y'*Y are derived by just substituting matrices
00800   // in to the above update for S'*Y:
00801   //
00802   // (S'*S) <- (S'*S) - S.row(jd)*S.row(jd)'
00803   //
00804   // (Y'*Y) <- (Y'*Y) - Y.row(jd)*Y.row(jd)'
00805   //
00806   // These updates are symmetric rank-1 updates.
00807   //
00808   DMatrixSlice S = this->S();
00809   DMatrixSlice Y = this->Y();
00810   DMatrixSlice STY = this->STY();
00811   DMatrixSliceSym        STS = this->STS();
00812   DMatrixSliceSym        YTY = this->YTY();
00813   // (S'*Y) <- (S'*Y) - S.row(jd)*Y.row(jd)'
00814   DenseLinAlgPack::ger( -1.0, S.row(jd), Y.row(jd), &STY );
00815   // (S'*S) <- (S'*S) - S.row(jd)*S.row(jd)'
00816   DenseLinAlgPack::syr( -1.0, S.row(jd), &STS );
00817   // (Y'*Y) <- (Y'*Y) - Y.row(jd)*Y.row(jd)'
00818   DenseLinAlgPack::syr( -1.0, Y.row(jd), &YTY );
00819   // Remove row jd from S and Y one column at a time
00820   // (one element at a time!)
00821   if( jd < n_ ) {
00822     {for( size_type k = 1; k <= m_bar_; ++k ) {
00823       value_type *ptr = S.col_ptr(k);
00824       std::copy( ptr + jd, ptr + n_, ptr + jd - 1 );
00825     }}
00826     {for( size_type k = 1; k <= m_bar_; ++k ) {
00827       value_type *ptr = Y.col_ptr(k);
00828       std::copy( ptr + jd, ptr + n_, ptr + jd - 1 );
00829     }}
00830   }
00831   // Update the size
00832   --n_;
00833   Q_updated_ = false;
00834 }
00835
00836 // Private member functions
00837
00838 void MatrixSymPosDefLBFGS::Vp_DtV( DVectorSlice* y, const DVectorSlice& x ) const
00839 {
00840   DenseLinAlgPack::Vp_MtV_assert_sizes( y->size(), m_bar_, m_bar_
00841     , BLAS_Cpp::no_trans, x.size() );
00842
00843   DVectorSlice::const_iterator
00844     d_itr = STY_.diag(0).begin(),
00845     x_itr = x.begin();
00846   DVectorSlice::iterator
00847     y_itr = y->begin();
00848
00849   while( y_itr != y->end() )
00850     *y_itr++ += (*d_itr++) * (*x_itr++);
00851 }
00852
00853 //
00854 // We need to update the factorizations to solve for:
00855 //
00856 // x = inv(Q) * y   =>   Q * x = y
00857 //
00858 //  [ (1/gk)*S'S   L  ] * [ x1 ] = [ y1 ]
00859 //  [      L'   -D  ]   [ x2 ]   [ y2 ]
00860 //
00861 // We will solve the above system using a schur complement:
00862 //
00863 // C = (1/gk)*S'S + L*inv(D)*L'
00864 //
00865 // According to the referenced paper, C is p.d. so:
00866 //
00867 // C = J*J'
00868 //
00869 // We then compute the solution as:
00870 //
00871 // x1 = inv(C) * ( y1 + L*inv(D)*y2 )
00872 // x2 = - inv(D) * ( y2 - L'*x1 )
00873 //
00874 // Therefore we will just update the factorization C = J*J'
00875 // where the factor J is stored in QJ_.
00876 //
00877
00878 void MatrixSymPosDefLBFGS::update_Q() const
00879 {
00880   using DenseLinAlgPack::tri;
00881   using DenseLinAlgPack::tri_ele;
00882   using DenseLinAlgPack::Mp_StM;
00883
00884   //
00885   // We need update the factorizations to solve for:
00886   //
00887   // x = inv(Q) * y
00888   //
00889   //  [ y1 ]  = [ (1/gk)*S'S   L  ] * [ x1 ]
00890   //  [ y2 ]    [      L'   -D  ]   [ x2 ]
00891   //
00892   // We will solve the above system using the schur complement:
00893   //
00894   // C = (1/gk)*S'S + L*inv(D)*L'
00895   //
00896   // According to the referenced paper, C is p.d. so:
00897   //
00898   // C = J*J'
00899   //
00900   // We then compute the solution as:
00901   //
00902   // x1 = inv(C) * ( y1 + L*inv(D)*y2 )
00903   // x2 = - inv(D) * ( y2 - L'*x1 )
00904   //
00905   // Therefore we will just update the factorization C = J*J'
00906   //
00907
00908   // Form the upper triangular part of C which will become J
00909   // which we are using storage of QJ
00910
00911   if( QJ_.rows() < m_ )
00912     QJ_.resize( m_, m_ );
00913
00914   const size_type
00915     mb = m_bar_;
00916
00917   DMatrixSlice
00918     C = QJ_(1,mb,1,mb);
00919
00920   // C = L * inv(D) * L'
00921   //
00922   // Here L is a strictly lower triangular (zero diagonal) matrix where:
00923   //
00924   // L = [ 0  0 ]
00925   //     [ Lb 0 ]
00926   //
00927   // Lb is lower triangular (nonzero diagonal)
00928   //
00929   // Therefore we compute L*inv(D)*L' as:
00930   //
00931   // C = [ 0  0 ] * [ Db  0  ] * [ 0  Lb' ]
00932   //     [ Lb 0 ]   [ 0   db ]   [ 0   0  ]
00933   //
00934   //   = [ 0  0  ] = [ 0      0     ]
00935   //     [ 0  Cb ]   [ 0  Lb*Db*Lb' ]
00936   //
00937   // We need to compute the upper triangular part of Cb.
00938
00939   C.row(1) = 0.0;
00940   if( mb > 1 )
00941     comp_Cb( STY_(2,mb,1,mb-1), STY_.diag(0)(1,mb-1), &C(2,mb,2,mb) );
00942
00943   // C += (1/gk)*S'S
00944
00945   const DMatrixSliceSym &STS = this->STS();
00946   Mp_StM( &C, (1/gamma_k_), tri( STS.gms(), STS.uplo(), BLAS_Cpp::nonunit )
00947     , BLAS_Cpp::trans );
00948
00949   // Now perform a cholesky factorization of C
00950   // After this factorization the upper triangular part of QJ
00951   // (through C) will contain the cholesky factor.
00952
00953   DMatrixSliceTriEle C_upper = tri_ele( C, BLAS_Cpp::upper );
00954   DenseLinAlgLAPack::potrf( &C_upper );
00955
00956   Q_updated_ = true;
00957 }
00958
00959 void MatrixSymPosDefLBFGS::V_invQtV( DVectorSlice* x, const DVectorSlice& y ) const
00960 {
00961   using DenseLinAlgPack::sym;
00962   using DenseLinAlgPack::tri;
00963   using DenseLinAlgPack::Vp_StV;
00964   using DenseLinAlgPack::V_InvMtV;
00965
00966   using LinAlgOpPack::Vp_V;
00967   using LinAlgOpPack::V_MtV;
00968
00969
00970   // Solve the system
00971   //
00972   // Q * x = y
00973   //
00974   // Using the schur complement factorization as described above.
00975
00976   const size_type
00977     mb = m_bar_;
00978
00979   if(!Q_updated_) {
00980     update_Q();
00981   }
00982
00983   DVectorSlice
00984     x1 = (*x)(1,mb),
00985     x2 = (*x)(mb+1,2*mb);
00986
00987   const DVectorSlice
00988     y1 = y(1,mb),
00989     y2 = y(mb+1,2*mb);
00990
00991   // //////////////////////////////////////
00992   // x1 = inv(C) * ( y1 + L*inv(D)*y2 )
00993   //     = inv(J'*J) * r
00994   //     = inv(J) * inv(J') * r
00995
00996   { // x1 = inv(D) * y2
00997     DVectorSlice::const_iterator
00998       d_itr = STY_.diag(0).begin(),
00999       y2_itr = y2.begin();
01000     DVectorSlice::iterator
01001       x1_itr = x1.begin();
01002     while( x1_itr != x1.end() )
01003       *x1_itr++ = *y2_itr++ / *d_itr++;
01004   }
01005
01006   // x1 = L * x1
01007   //
01008   //    = [ 0  0 ] * [ x1(1:mb-1) ]
01009   //      [ Lb 0 ]   [ x1(mb)     ]
01010   //
01011   //    = [ 0             ]
01012   //      [ Lb*x1(1:mb-1) ]
01013   //
01014   if( mb > 1 ) {
01015     // x1(2,mb) = x1(1,mb-1) ( copy from mb-1 to mb then mb-2 to mb-1
01016     // etc. so that we don't overwrite the elements we need to copy ).
01017     DVectorSlice
01018       x1a = x1(1,mb-1),
01019       x1b = x1(2,mb);
01020     std::copy( x1a.rbegin(), x1a.rend(), x1b.rbegin() );
01021     V_MtV( &x1b, Lb(), BLAS_Cpp::no_trans, x1b );
01022   }
01023   x1(1) = 0.0;
01024
01025   // r = x1 += y1
01026   Vp_V( &x1, y1 );
01027
01028   // x1 = inv(J') * r
01029   const DMatrixSliceTri J = tri( QJ_(1,mb,1,mb), BLAS_Cpp::upper, BLAS_Cpp::nonunit );
01030   V_InvMtV( &x1, J, BLAS_Cpp::trans, x1 );
01031
01032   // x1 = inv(J) * x1
01033   V_InvMtV( &x1, J, BLAS_Cpp::no_trans, x1 );
01034
01035   // /////////////////////////////////////
01036   // x2 = inv(D) * ( - y2 + L'*x1 )
01037
01038   // x2 = L'*x1
01039   //
01040   //    = [ 0  Lb' ] * [ x1(1)    ]
01041   //      [ 0  0   ]   [ x1(2,mb) ]
01042   //
01043   //    = [ Lb'*x1(2,mb) ]
01044   //      [      0       ]
01045   //
01046   if( mb > 1 ) {
01047     V_MtV( &x2(1,mb-1), Lb(), BLAS_Cpp::trans, x1(2,mb) );
01048   }
01049   x2(mb) = 0.0;
01050
01051   // x2 += -y2
01052   Vp_StV( &x2, -1.0, y2 );
01053
01054   // x2 = inv(D) * x2
01055   {
01056     DVectorSlice::const_iterator
01057       d_itr = STY_.diag(0).begin();
01058     DVectorSlice::iterator
01059       x2_itr = x2.begin();
01060     while( x2_itr != x2.end() )
01061       *x2_itr++ /= *d_itr++;
01062   }
01063 }
01064
01065 void MatrixSymPosDefLBFGS::assert_initialized() const
01066 {
01067   if(!n_)
01068     throw std::logic_error( "MatrixSymPosDefLBFGS::assert_initialized() : "
01069       "Error, matrix not initialized" );
01070 }
01071
01072 } // end namespace ConstrainedOptPack
01073
01074 namespace {
01075
01076 void comp_Cb( const DMatrixSlice& Lb, const DVectorSlice& Db_diag
01077   , DMatrixSlice* Cb )
01078 {
01079   // Lb * inv(Db) * Lb =
01080   //
01081   // [ l(1,1)           ]   [ dd(1)         ]   [ l(1,1)  l(2,1)  ... l(p,1)  ]
01082   // [ l(2,1) l(2,2)        ]   [   dd(2)     ]   [     l(2,2)  ... l(p,2)  ]
01083   // [ .    .        .      ] * [     .     ] * [         . .   ]
01084   // [ l(p,1) l(p,2)  ... l(p,p)  ]   [       dd(p) ]   [           l(p,p)  ]
01085   //
01086   //
01087   // [ l(1,1)*dd(1)*l(1,1)  l(1,1)*dd(1)*l(2,1)             ...   l(1,1)*dd(1)*l(p,1)             ]
01088   // [ 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) ]
01089   // [ .            .                     ...
01090   // [ symmetric        symmetric                 ...   sum( l(p,i)*dd(i)*l(p,i), i=1,..,p )    ]
01091   //
01092   // Therefore we can express the upper triangular elemetns of Cb as:
01093   //
01094   // Cb(i,j) = sum( l(i,k)*dd(k)*l(j,k), k = 1,..,i )
01095
01096   typedef DenseLinAlgPack::size_type size_type;
01097   typedef DenseLinAlgPack::value_type value_type;
01098
01099   assert( Lb.rows() == Cb->rows() && Cb->rows() == Db_diag.size() );
01100
01101   const size_type p = Db_diag.size();
01102
01103   for( size_type i = 1; i <= p; ++i ) {
01104     for( size_type j = i; j <= p; ++j ) {
01105       value_type &c = (*Cb)(i,j) = 0.0;
01106       for( size_type k = 1; k <= i; ++k ) {
01107         c += Lb(i,k) * Lb(j,k) / Db_diag(k);
01108       }
01109     }
01110   }
01111
01112   // ToDo: Make the above operation more efficent if needed!
01113 }
01114
01115 } // end namespace
```

Generated on Thu Sep 18 12:34:15 2008 for ConstrainedOptPack: C++ Tools for Constrained (and Unconstrained) Optimization by  1.3.9.1