ConstrainedOptPack_MatrixSymPosDefLBFGSSerial.cpp

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