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