ConstrainedOptPack_MatrixSymPosDefLBFGS.cpp

Go to the documentation of this file.
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 <typeinfo>
00062 
00063 #include "ConstrainedOptPack_MatrixSymPosDefLBFGS.hpp"
00064 #include "AbstractLinAlgPack_BFGS_helpers.hpp"
00065 #include "AbstractLinAlgPack_VectorDenseEncap.hpp"
00066 #include "AbstractLinAlgPack_MatrixOpOut.hpp"
00067 #include "AbstractLinAlgPack_VectorStdOps.hpp"
00068 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
00069 #include "DenseLinAlgPack_LinAlgOpPack.hpp"
00070 #include "DenseLinAlgPack_DMatrixOut.hpp"
00071 #include "DenseLinAlgLAPack.hpp"
00072 #include "Teuchos_Workspace.hpp"
00073 #include "Teuchos_TestForException.hpp"
00074 
00075 namespace {
00076 
00077   using DenseLinAlgPack::DVectorSlice;
00078   using DenseLinAlgPack::DMatrixSlice;
00079 
00087   void comp_Cb( const DMatrixSlice& Lb, const DVectorSlice& Db_diag
00088     , DMatrixSlice* Cb );
00089 
00090 } // end namespace
00091 
00092 namespace ConstrainedOptPack {
00093 
00094 // /////////////////////////////////
00095 // Inline private member functions
00096 
00097 inline
00098 const DMatrixSliceTri MatrixSymPosDefLBFGS::R() const
00099 {
00100   return DenseLinAlgPack::tri( STY_(1,m_bar_,1,m_bar_), BLAS_Cpp::upper, BLAS_Cpp::nonunit );
00101 }
00102 
00103 inline
00104 const DMatrixSliceTri MatrixSymPosDefLBFGS::Lb() const
00105 {
00106   return DenseLinAlgPack::tri( STY_(2,m_bar_,1,m_bar_-1), BLAS_Cpp::lower, BLAS_Cpp::nonunit );
00107 }
00108 
00109 inline
00110 DMatrixSlice MatrixSymPosDefLBFGS::STY()
00111 {
00112   return STY_(1,m_bar_,1,m_bar_);
00113 }
00114 
00115 inline
00116 const DMatrixSlice MatrixSymPosDefLBFGS::STY() const
00117 {
00118   return STY_(1,m_bar_,1,m_bar_);
00119 }
00120 
00121 inline
00122 DMatrixSliceSym MatrixSymPosDefLBFGS::STS()
00123 {
00124   return DenseLinAlgPack::nonconst_sym( STSYTY_(2,m_bar_+1,1,m_bar_),BLAS_Cpp::lower );
00125 }
00126 
00127 inline
00128 const DMatrixSliceSym MatrixSymPosDefLBFGS::STS() const
00129 {
00130   return DenseLinAlgPack::sym( STSYTY_(2,m_bar_+1,1,m_bar_),BLAS_Cpp::lower );
00131 }
00132 
00133 inline
00134 DMatrixSliceSym MatrixSymPosDefLBFGS::YTY()
00135 {
00136   return DenseLinAlgPack::nonconst_sym( STSYTY_(1,m_bar_,2,m_bar_+1),BLAS_Cpp::upper );
00137 }
00138 
00139 inline
00140 const DMatrixSliceSym MatrixSymPosDefLBFGS::YTY() const
00141 {
00142   return DenseLinAlgPack::sym( STSYTY_(1,m_bar_,2,m_bar_+1),BLAS_Cpp::upper );
00143 }
00144 
00145 // ///////////////////////
00146 // Nonlinined functions
00147 
00148 MatrixSymPosDefLBFGS::MatrixSymPosDefLBFGS(
00149     size_type   m
00150   ,bool       maintain_original
00151   ,bool       maintain_inverse
00152   ,bool       auto_rescaling
00153   )
00154 {
00155   initial_setup(m,maintain_original,maintain_inverse,auto_rescaling);
00156 }
00157 
00158 void MatrixSymPosDefLBFGS::initial_setup(
00159   size_type m,
00160   bool maintain_original,
00161   bool maintain_inverse,
00162   bool auto_rescaling
00163   )
00164 {
00165   // Validate input
00166   TEST_FOR_EXCEPTION(
00167     !maintain_original && !maintain_inverse, std::invalid_argument
00168     ,"MatrixSymPosDefLBFGS::initial_setup(...) : "
00169     "Error, both maintain_original and maintain_inverse can not both be false!" );
00170   TEST_FOR_EXCEPTION(
00171     m < 1, std::invalid_argument
00172     ,"MatrixSymPosDefLBFGS::set_num_updates_stored(m) : "
00173     "Error, the number of storage locations must be > 0" );
00174   vec_spc_           = Teuchos::null;
00175   maintain_original_ = maintain_original;
00176   maintain_inverse_  = maintain_inverse;
00177   m_                 = m;
00178   n_                 = 0; // make uninitialized
00179   num_secant_updates_= 0;
00180   auto_rescaling_    = auto_rescaling;
00181 }
00182 
00183 // Overridden from MatrixOp
00184 
00185 const VectorSpace& MatrixSymPosDefLBFGS::space_cols() const
00186 {
00187   return *vec_spc_;
00188 }
00189 
00190 std::ostream& MatrixSymPosDefLBFGS::output(std::ostream& out) const
00191 {
00192   assert_initialized();
00193   out << "\n*** Limited Memory BFGS matrix\n";
00194   out << "\nConversion to dense =\n";
00195   MatrixOp::output(out);
00196   out << "\nStored quantities\n"
00197     << "\nn       = " << n_
00198     << "\nm       = " << m_
00199     << "\nm_bar   = " << m_bar_
00200     << "\ngamma_k = " << gamma_k_ << std::endl;
00201   if( m_bar_ ) {
00202     out << "\nS =\n" << *S()
00203       << "\nY =\n" << *Y()
00204       << "\nS'Y =\n" << STY_(1,m_bar_,1,m_bar_)
00205       << "\nlower(S'S) \\ zero diagonal \\ upper(Y'Y) =\n"
00206         << STSYTY_(1,m_bar_+1,1,m_bar_+1)
00207       << "\nQ updated? = " << Q_updated_ << std::endl;
00208     if(Q_updated_)
00209       out << "\nCholesky of schur complement of Q, QJ =\n" << QJ_(1,m_bar_,1,m_bar_);
00210   }
00211   return out;
00212 }
00213 
00214 MatrixOp& MatrixSymPosDefLBFGS::operator=(const MatrixOp& mwo)
00215 { 
00216   const MatrixSymPosDefLBFGS *p_m = dynamic_cast<const MatrixSymPosDefLBFGS*>(&mwo);
00217   if(p_m) {
00218     if( p_m == this ) return *this; // assignment to self
00219     // Important: Assign all members here.
00220     auto_rescaling_      = p_m->auto_rescaling_;
00221     maintain_original_   = p_m->maintain_original_;
00222     original_is_updated_ = p_m->original_is_updated_;
00223     maintain_inverse_    = p_m->maintain_inverse_;
00224     inverse_is_updated_  = p_m->inverse_is_updated_;
00225     vec_spc_             = p_m->vec_spc_.get() ? p_m->vec_spc_->clone() : Teuchos::null;
00226     n_               = p_m->n_;
00227     m_               = p_m->m_;
00228     m_bar_             = p_m->m_bar_;
00229     num_secant_updates_  = p_m->num_secant_updates_;
00230     gamma_k_           = p_m->gamma_k_;
00231     S_               = p_m->S_;
00232     Y_               = p_m->Y_;
00233     STY_             = p_m->STY_;
00234     STSYTY_            = p_m->STSYTY_;
00235     Q_updated_           = p_m->Q_updated_;
00236     QJ_              = p_m->QJ_;
00237   }
00238   else {
00239     TEST_FOR_EXCEPTION(
00240       true,std::invalid_argument
00241       ,"MatrixSymPosDefLBFGS::operator=(const MatrixOp& mwo) : Error, "
00242       "The concrete type of mwo \'" << typeName(mwo) << "\' is not "
00243       "as subclass of MatrixSymPosDefLBFGS as required" );
00244   }
00245   return *this;
00246 }
00247 
00248 // Level-2 BLAS
00249 
00250 void MatrixSymPosDefLBFGS::Vp_StMtV(
00251     VectorMutable* y, value_type alpha, BLAS_Cpp::Transp trans_rhs1
00252   , const Vector& x, value_type beta
00253   ) const
00254 {
00255   using AbstractLinAlgPack::Vt_S;
00256   using AbstractLinAlgPack::Vp_StV;
00257   using AbstractLinAlgPack::Vp_StMtV;
00258   using LinAlgOpPack::V_StMtV;
00259   using LinAlgOpPack::V_MtV;
00260   typedef VectorDenseEncap         vde;
00261   typedef VectorDenseMutableEncap  vdme;
00262   using Teuchos::Workspace;
00263   Teuchos::WorkspaceStore* wss = Teuchos::get_default_workspace_store().get();
00264 
00265   assert_initialized();
00266 
00267   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   Vt_S( y, beta );
00293 
00294   // y += (1/gk)*x
00295   Vp_StV( y, invgk, x );
00296 
00297   if( !m_bar_ )
00298     return; // No updates have been added yet.
00299 
00300   const multi_vec_ptr_t
00301     S = this->S(),
00302     Y = this->Y();
00303 
00304   // Get workspace
00305 
00306   const size_type
00307     mb = m_bar_;
00308 
00309   Workspace<value_type>  t1_ws(wss,2*mb);
00310   DVectorSlice                 t1(&t1_ws[0],t1_ws.size());
00311   Workspace<value_type>  t2_ws(wss,2*mb);
00312   DVectorSlice                 t2(&t2_ws[0],t2_ws.size());
00313 
00314   VectorSpace::vec_mut_ptr_t
00315     t = S->space_rows().create_member();
00316 
00317   // t1 = [ (1/gk)*S'*x ]
00318   //    [      Y'*x   ]
00319 
00320   V_StMtV( t.get(), invgk, *S, BLAS_Cpp::trans, x );
00321   t1(1,mb) = vde(*t)();
00322   V_MtV( t.get(), *Y, BLAS_Cpp::trans, x );
00323   t1(mb+1,2*mb) = vde(*t)();
00324 
00325   // t2 = inv(Q) * t1
00326   V_invQtV( &t2, t1 );
00327 
00328   // y += -(1/gk) * S * t2(1:m)
00329   (vdme(*t)() = t2(1,mb));
00330   Vp_StMtV( y, -invgk, *S, BLAS_Cpp::no_trans, *t );
00331 
00332   // y += -1.0 * Y * t2(m+1,2m
00333   (vdme(*t)() = t2(mb+1,2*mb));
00334   Vp_StMtV( y, -1.0, *Y, BLAS_Cpp::no_trans, *t );
00335 
00336 }
00337 
00338 // Overridden from MatrixOpNonsing
00339 
00340 void MatrixSymPosDefLBFGS::V_InvMtV(
00341   VectorMutable* y, BLAS_Cpp::Transp trans_rhs1
00342   , const Vector& x
00343   ) const
00344 {
00345   using AbstractLinAlgPack::Vp_StMtV;
00346   using DenseLinAlgPack::V_InvMtV;
00347   using LinAlgOpPack::V_mV;
00348   using LinAlgOpPack::V_StV;
00349   using LinAlgOpPack::Vp_V;
00350   using LinAlgOpPack::V_MtV;
00351   using LinAlgOpPack::V_StMtV;
00352   using LinAlgOpPack::Vp_MtV;
00353   using DenseLinAlgPack::Vp_StMtV;
00354   typedef VectorDenseEncap         vde;
00355   typedef VectorDenseMutableEncap  vdme;
00356   using Teuchos::Workspace;
00357   Teuchos::WorkspaceStore* wss = Teuchos::get_default_workspace_store().get();
00358 
00359   assert_initialized();
00360 
00361   TEST_FOR_EXCEPT( !(  inverse_is_updated_  ) ); // For now just always update
00362 
00363   // y = inv(Bk) * x = Hk * x
00364   //
00365   // = gk*x + [S gk*Y] * [ inv(R')*(D+gk*Y'Y)*inv(R)     -inv(R') ] * [   S'  ] * x
00366   //                     [            -inv(R)                0    ]   [ gk*Y' ]
00367   //
00368   // Perform in the following (in order):
00369   //
00370   // y = gk*x
00371   //
00372   // t1 = [   S'*x  ]         <: R^(2*m)
00373   //      [ gk*Y'*x ]
00374   //
00375   // t2 = inv(R) * t1(1:m)      <: R^(m)
00376   //
00377   // t3 = - inv(R') * t1(m+1,2*m)   <: R^(m)
00378   //
00379   // t4 = gk * Y'Y * t2       <: R^(m)
00380   //
00381   // t4 += D*t2
00382   //
00383   // t5 = inv(R') * t4        <: R^(m)
00384   //
00385   // t5 += t3
00386   //
00387   // y += S*t5
00388   //
00389   // y += -gk*Y*t2
00390 
00391   // y = gk*x
00392   V_StV( y, gamma_k_, x );
00393 
00394   const size_type
00395     mb = m_bar_;  
00396   
00397   if( !mb )
00398     return; // No updates have been performed.
00399 
00400   const multi_vec_ptr_t
00401     S = this->S(),
00402     Y = this->Y();
00403 
00404   // Get workspace
00405 
00406   Workspace<value_type>    t1_ws(wss,2*mb);
00407   DVectorSlice                   t1(&t1_ws[0],t1_ws.size());
00408   Workspace<value_type>    t2_ws(wss,mb);
00409   DVectorSlice                   t2(&t2_ws[0],t2_ws.size());
00410   Workspace<value_type>    t3_ws(wss,mb);
00411   DVectorSlice                   t3(&t3_ws[0],t3_ws.size());
00412   Workspace<value_type>    t4_ws(wss,mb);
00413   DVectorSlice                   t4(&t4_ws[0],t4_ws.size());
00414   Workspace<value_type>    t5_ws(wss,mb);
00415   DVectorSlice                   t5(&t5_ws[0],t5_ws.size());
00416 
00417   VectorSpace::vec_mut_ptr_t
00418     t = S->space_rows().create_member();
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( t.get(), *S, BLAS_Cpp::trans, x );
00429   t1(1,mb) = vde(*t)();
00430   V_StMtV( t.get(), gamma_k_, *Y, BLAS_Cpp::trans, x );
00431   t1(mb+1,2*mb) = vde(*t)();
00432 
00433   // t2 = inv(R) * t1(1:m)
00434   V_InvMtV( &t2, R, BLAS_Cpp::no_trans, t1(1,mb) );
00435 
00436   // t3 = - inv(R') * t1(m+1,2*m)
00437   V_mV( &t3, t1(mb+1,2*mb) );
00438   V_InvMtV( &t3, R, BLAS_Cpp::trans, t3 );
00439 
00440   // t4 = gk * Y'Y * t2
00441   V_StMtV( &t4, gamma_k_, YTY, BLAS_Cpp::no_trans, t2 );
00442 
00443   // t4 += D*t2
00444   Vp_DtV( &t4, t2 );
00445 
00446   // t5 = inv(R') * t4
00447   V_InvMtV( &t5, R, BLAS_Cpp::trans, t4 );
00448 
00449   // t5 += t3
00450   Vp_V( &t5, t3 );
00451 
00452   // y += S*t5
00453   (vdme(*t)() = t5);
00454   Vp_MtV( y, *S, BLAS_Cpp::no_trans, *t );
00455 
00456   // y += -gk*Y*t2
00457   (vdme(*t)() = t2);
00458   Vp_StMtV( y, -gamma_k_, *Y, BLAS_Cpp::no_trans, *t );
00459 
00460 }
00461 
00462 // Overridden from MatrixSymSecant
00463 
00464 void MatrixSymPosDefLBFGS::init_identity( const VectorSpace& space_diag, value_type alpha )
00465 {
00466   // Validate input
00467   TEST_FOR_EXCEPTION(
00468     alpha <= 0.0, std::invalid_argument
00469     ,"MatrixSymPosDefLBFGS::init_identity(n,alpha) : Error, "
00470     "alpha = " << alpha << " <= 0 is not allowed!" );
00471   
00472   // Set the vector space
00473   vec_spc_ = space_diag.clone();
00474   vec_spc_.get();
00475 
00476   // Set storage
00477   S_ = vec_spc_->create_members(m_);
00478   Y_ = vec_spc_->create_members(m_);
00479   TEST_FOR_EXCEPT( !( S_.get() ) );
00480   TEST_FOR_EXCEPT( !( Y_.get() ) );
00481   STY_.resize( m_, m_ );
00482   STSYTY_.resize( m_+1, m_+1 );
00483   STSYTY_.diag(0) = 0.0;
00484 
00485   gamma_k_ = 1.0/alpha;
00486 
00487   // Initialize counters
00488   m_bar_  = 0;
00489 
00490   n_ = vec_spc_->dim();        // initialized;
00491   original_is_updated_ = true; // This will never change for now
00492   inverse_is_updated_  = true; // This will never change for now
00493   num_secant_updates_  = 0;    // reset this to zero
00494 }
00495 
00496 void MatrixSymPosDefLBFGS::init_diagonal(const Vector& diag)
00497 {
00498   init_identity( diag.space(), diag.norm_inf() );
00499 }
00500 
00501 void MatrixSymPosDefLBFGS::secant_update(
00502   VectorMutable* s, VectorMutable* y, VectorMutable* Bs
00503   )
00504 {
00505   using AbstractLinAlgPack::BFGS_sTy_suff_p_d;
00506   using AbstractLinAlgPack::dot;
00507   using LinAlgOpPack::V_MtV;
00508   using Teuchos::Workspace;
00509   Teuchos::WorkspaceStore* wss = Teuchos::get_default_workspace_store().get();
00510 
00511   assert_initialized();
00512 
00513   // Check skipping the BFGS update
00514   const value_type
00515     sTy       = dot(*s,*y);
00516   std::ostringstream omsg;
00517   if( !BFGS_sTy_suff_p_d(*s,*y,&sTy,&omsg,"MatrixSymPosDefLBFGS::secant_update(...)" ) ) {
00518     throw UpdateSkippedException( omsg.str() ); 
00519   }
00520 
00521   try {
00522 
00523   // Update counters
00524   if( m_bar_ == m_ ) {
00525     // We are at the end of the storage so remove the oldest stored update
00526     // and move updates to make room for the new update.  This has to be done for the
00527     // the matrix to behave properly
00528     {for( size_type k = 1; k <= m_-1; ++k ) {
00529       S_->col(k) = S_->col(k+1);                            // Shift S.col() to the left
00530       Y_->col(k) = Y_->col(k+1);                            // Shift Y.col() to the left
00531       STY_.col(k)(1,m_-1) = STY_.col(k+1)(2,m_);            // Move submatrix STY(2,m-1,2,m-1) up and left
00532       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
00533       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
00534     }}
00535     // ToDo: Create an abstract interface, call it MultiVectorShiftVecs, to rearrange S and Y all at once.
00536     // This will be important for parallel performance.
00537   }
00538   else {
00539     m_bar_++;
00540   }
00541   // Set the update vectors
00542   *S_->col(m_bar_) = *s;
00543   *Y_->col(m_bar_) = *y;
00544 
00545   // /////////////////////////////////////////////////////////////////////////////////////
00546   // Update S'Y
00547   //
00548   // Update the row and column m_bar
00549   //
00550   //  S'Y = 
00551   //
00552   //  [ s(1)'*y(1)    ...   s(1)'*y(m_bar)    ...   s(1)'*y(m_bar)    ]
00553   //  [ .           .             .         ] row
00554   //  [ s(m_bar)'*y(1)  ...   s(m_bar)'*y(m_bar)  ...   s(m_bar)'*y(m_bar)  ] m_bar
00555   //  [ .           .             .         ]
00556   //  [ s(m_bar)'*y(1)  ...   s(m_bar)'*y(m_bar)  ...   s(m_bar)'*y(m_bar)  ]
00557   //
00558   //                col m_bar
00559   //
00560   // Therefore we set:
00561   //  (S'Y)(:,m_bar) =  S'*y(m_bar)
00562   //  (S'Y)(m_bar,:) =  s(m_bar)'*Y
00563 
00564   const multi_vec_ptr_t
00565     S = this->S(),
00566     Y = this->Y();
00567 
00568   VectorSpace::vec_mut_ptr_t
00569     t = S->space_rows().create_member();  // temporary workspace
00570 
00571   //  (S'Y)(:,m_bar) =  S'*y(m_bar)
00572   V_MtV( t.get(), *S, BLAS_Cpp::trans, *y );
00573   STY_.col(m_bar_)(1,m_bar_) = VectorDenseEncap(*t)();
00574 
00575   //  (S'Y)(m_bar,:)' =  Y'*s(m_bar)
00576   V_MtV( t.get(), *Y, BLAS_Cpp::trans, *s );
00577   STY_.row(m_bar_)(1,m_bar_) = VectorDenseEncap(*t)();
00578 
00579   // /////////////////////////////////////////////////////////////////
00580   // Update S'S
00581   //
00582   //  S'S = 
00583   //
00584   //  [ s(1)'*s(1)    ...   symmetric         symmetric     ]
00585   //  [ .           .             .         ] row
00586   //  [ s(m_bar)'*s(1)  ...   s(m_bar)'*s(m_bar)  ...   symmetric     ] m_bar
00587   //  [ .           .             .         ]
00588   //  [ s(m_bar)'*s(1)  ...   s(m_bar)'*s(m_bar)  ...   s(m_bar)'*s(m_bar)  ]
00589   //
00590   //                col m_bar
00591   //
00592   // Here we will update the lower triangular part of S'S.  To do this we
00593   // only need to compute:
00594   //    t = S'*s(m_bar) = { s(m_bar)' * [ s(1),..,s(m_bar),..,s(m_bar) ]  }'
00595   // then set the appropriate rows and columns of S'S.
00596 
00597   Workspace<value_type>   work_ws(wss,m_bar_);
00598   DVectorSlice                  work(&work_ws[0],work_ws.size());
00599 
00600   // work = S'*s(m_bar)
00601   V_MtV( t.get(), *S, BLAS_Cpp::trans, *s );
00602   work = VectorDenseEncap(*t)();
00603 
00604   // Set row elements
00605   STSYTY_.row(m_bar_+1)(1,m_bar_) = work;
00606   // Set column elements
00607   STSYTY_.col(m_bar_)(m_bar_+1,m_bar_+1) = work(m_bar_,m_bar_);
00608 
00609   // /////////////////////////////////////////////////////////////////////////////////////
00610   // Update Y'Y
00611   //
00612   // Update the row and column m_bar
00613   //
00614   //  Y'Y = 
00615   //
00616   //  [ y(1)'*y(1)    ...   y(1)'*y(m_bar)    ...   y(1)'*y(m_bar)    ]
00617   //  [ .           .             .         ] row
00618   //  [ symmetric   ...   y(m_bar)'*y(m_bar)  ...   y(m_bar)'*y(m_bar)  ] m_bar
00619   //  [ .           .             .         ]
00620   //  [ symmetric   ...   symmetric     ...   y(m_bar)'*y(m_bar)  ]
00621   //
00622   //                col m_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(m_bar) = { y(m_bar)' * [ y(1),..,y(m_bar),..,y(m_bar) ]  }'
00627   // then set the appropriate rows and columns of Y'Y.
00628 
00629   // work = Y'*y(m_bar)
00630   V_MtV( t.get(), *Y, BLAS_Cpp::trans, *y );
00631   work = VectorDenseEncap(*t)();
00632 
00633   // Set row elements
00634   STSYTY_.col(m_bar_+1)(1,m_bar_) = work;
00635   // Set column elements
00636   STSYTY_.row(m_bar_)(m_bar_+1,m_bar_+1) = work(m_bar_,m_bar_);
00637 
00638   // /////////////////////////////
00639   // Update gamma_k
00640 
00641   // gamma_k = s'*y / y'*y
00642   if(auto_rescaling_)
00643     gamma_k_ = STY_(m_bar_,m_bar_) / STSYTY_(m_bar_,m_bar_+1);
00644 
00645   // We do not initially update Q unless we have to form a matrix-vector
00646   // product later.
00647   
00648   Q_updated_ = false;
00649   num_secant_updates_++;
00650 
00651   } //  end try
00652   catch(...) {
00653     // If we throw any exception the we should make the matrix uninitialized
00654     // so that we do not leave this object in an inconsistant state.
00655     n_ = 0;
00656     throw;
00657   }
00658 
00659 }
00660 
00661 // Private member functions
00662 
00663 void MatrixSymPosDefLBFGS::Vp_DtV( DVectorSlice* y, const DVectorSlice& x ) const
00664 {
00665   DenseLinAlgPack::Vp_MtV_assert_sizes(
00666     y->dim(), m_bar_, m_bar_, BLAS_Cpp::no_trans, x.dim() );
00667 
00668   DVectorSlice::const_iterator
00669     d_itr = STY_.diag(0).begin(),
00670     x_itr = x.begin();
00671   DVectorSlice::iterator
00672     y_itr = y->begin();
00673 
00674   while( y_itr != y->end() )
00675     *y_itr++ += (*d_itr++) * (*x_itr++);    
00676 }
00677 
00678 //
00679 // We need to update the factorizations to solve for:
00680 //
00681 // x = inv(Q) * y   =>   Q * x = y
00682 //
00683 //  [ (1/gk)*S'S   L  ] * [ x1 ] = [ y1 ]
00684 //  [      L'   -D  ]   [ x2 ]   [ y2 ]
00685 //
00686 // We will solve the above system using a schur complement:
00687 //
00688 // C = (1/gk)*S'S + L*inv(D)*L'
00689 //
00690 // According to the referenced paper, C is p.d. so:
00691 //
00692 // C = J*J'
00693 //
00694 // We then compute the solution as:
00695 //
00696 // x1 = inv(C) * ( y1 + L*inv(D)*y2 )
00697 // x2 = - inv(D) * ( y2 - L'*x1 )
00698 //
00699 // Therefore we will just update the factorization C = J*J'
00700 // where the factor J is stored in QJ_.
00701 //
00702 
00703 void MatrixSymPosDefLBFGS::update_Q() const
00704 {
00705   using DenseLinAlgPack::tri;
00706   using DenseLinAlgPack::tri_ele;
00707   using DenseLinAlgPack::Mp_StM;
00708 
00709   //
00710   // We need update the factorizations to solve for:
00711   //
00712   // x = inv(Q) * y
00713   //
00714   //  [ y1 ]  = [ (1/gk)*S'S   L  ] * [ x1 ]
00715   //  [ y2 ]    [      L'   -D  ]   [ x2 ]
00716   //
00717   // We will solve the above system using the schur complement:
00718   //
00719   // C = (1/gk)*S'S + L*inv(D)*L'
00720   //
00721   // According to the referenced paper, C is p.d. so:
00722   //
00723   // C = J*J'
00724   //
00725   // We then compute the solution as:
00726   //
00727   // x1 = inv(C) * ( y1 + L*inv(D)*y2 )
00728   // x2 = - inv(D) * ( y2 - L'*x1 )
00729   //
00730   // Therefore we will just update the factorization C = J*J'
00731   //
00732 
00733   // Form the upper triangular part of C which will become J
00734   // which we are using storage of QJ
00735 
00736   if( QJ_.rows() < m_ )
00737     QJ_.resize( m_, m_ );
00738 
00739   const size_type
00740     mb = m_bar_;
00741 
00742   DMatrixSlice
00743     C = QJ_(1,mb,1,mb);
00744 
00745   // C = L * inv(D) * L'
00746   //
00747   // Here L is a strictly lower triangular (zero diagonal) matrix where:
00748   //
00749   // L = [ 0  0 ]
00750   //     [ Lb 0 ]
00751   //
00752   // Lb is lower triangular (nonzero diagonal)
00753   //
00754   // Therefore we compute L*inv(D)*L' as:
00755   //
00756   // C = [ 0  0 ] * [ Db  0  ] * [ 0  Lb' ]
00757   //     [ Lb 0 ]   [ 0   db ]   [ 0   0  ]
00758   //
00759   //   = [ 0  0  ] = [ 0      0     ]
00760   //     [ 0  Cb ]   [ 0  Lb*Db*Lb' ]
00761   //
00762   // We need to compute the upper triangular part of Cb.
00763 
00764   C.row(1) = 0.0;
00765   if( mb > 1 )
00766     comp_Cb( STY_(2,mb,1,mb-1), STY_.diag(0)(1,mb-1), &C(2,mb,2,mb) );
00767 
00768   // C += (1/gk)*S'S
00769 
00770   const DMatrixSliceSym &STS = this->STS();
00771   Mp_StM( &C, (1/gamma_k_), tri( STS.gms(), STS.uplo(), BLAS_Cpp::nonunit )
00772     , BLAS_Cpp::trans );
00773 
00774   // Now perform a cholesky factorization of C
00775   // After this factorization the upper triangular part of QJ
00776   // (through C) will contain the cholesky factor.
00777 
00778   DMatrixSliceTriEle C_upper = tri_ele( C, BLAS_Cpp::upper );
00779   try {
00780     DenseLinAlgLAPack::potrf( &C_upper );
00781   }
00782   catch( const DenseLinAlgLAPack::FactorizationException &fe ) {
00783     TEST_FOR_EXCEPTION(
00784       true, UpdateFailedException
00785       ,"Error, the factorization of Q which should be s.p.d. failed with"
00786       " the error message: {" << fe.what() << "}";
00787       );
00788   }
00789 
00790   Q_updated_ = true;
00791 }
00792 
00793 void MatrixSymPosDefLBFGS::V_invQtV( DVectorSlice* x, const DVectorSlice& y ) const
00794 {
00795   using DenseLinAlgPack::sym;
00796   using DenseLinAlgPack::tri;
00797   using DenseLinAlgPack::Vp_StV;
00798   using DenseLinAlgPack::V_InvMtV;
00799 
00800   using LinAlgOpPack::Vp_V;
00801   using LinAlgOpPack::V_MtV;
00802 
00803 
00804   // Solve the system 
00805   //
00806   // Q * x = y
00807   //
00808   // Using the schur complement factorization as described above.
00809 
00810   const size_type
00811     mb = m_bar_;
00812 
00813   if(!Q_updated_) {
00814     update_Q();
00815   }
00816 
00817   DVectorSlice
00818     x1 = (*x)(1,mb),
00819     x2 = (*x)(mb+1,2*mb);
00820 
00821   const DVectorSlice
00822     y1 = y(1,mb),
00823     y2 = y(mb+1,2*mb);
00824 
00825   // //////////////////////////////////////
00826   // x1 = inv(C) * ( y1 + L*inv(D)*y2 )
00827   //     = inv(J'*J) * r
00828   //     = inv(J) * inv(J') * r
00829 
00830   { // x1 = inv(D) * y2
00831     DVectorSlice::const_iterator
00832       d_itr = STY_.diag(0).begin(),
00833       y2_itr = y2.begin();
00834     DVectorSlice::iterator
00835       x1_itr = x1.begin();
00836     while( x1_itr != x1.end() )
00837       *x1_itr++ = *y2_itr++ / *d_itr++;
00838   }
00839 
00840   // x1 = L * x1
00841   //
00842   //    = [ 0  0 ] * [ x1(1:mb-1) ]
00843   //      [ Lb 0 ]   [ x1(mb)     ]
00844   //
00845   //    = [ 0             ]
00846   //      [ Lb*x1(1:mb-1) ]
00847   //
00848   if( mb > 1 ) {
00849     // x1(2,mb) = x1(1,mb-1) ( copy from mb-1 to mb then mb-2 to mb-1
00850     // etc. so that we don't overwrite the elements we need to copy ).
00851     DVectorSlice
00852       x1a = x1(1,mb-1),
00853       x1b = x1(2,mb);
00854     std::copy( x1a.rbegin(), x1a.rend(), x1b.rbegin() );
00855     V_MtV( &x1b, Lb(), BLAS_Cpp::no_trans, x1b );
00856   }
00857   x1(1) = 0.0;
00858 
00859   // r = x1 += y1
00860   Vp_V( &x1, y1 );
00861 
00862   // x1 = inv(J') * r
00863   const DMatrixSliceTri J = tri( QJ_(1,mb,1,mb), BLAS_Cpp::upper, BLAS_Cpp::nonunit );
00864   V_InvMtV( &x1, J, BLAS_Cpp::trans, x1 );
00865 
00866   // x1 = inv(J) * x1
00867   V_InvMtV( &x1, J, BLAS_Cpp::no_trans, x1 );
00868 
00869   // /////////////////////////////////////
00870   // x2 = inv(D) * ( - y2 + L'*x1 )
00871 
00872   // x2 = L'*x1
00873   //
00874   //    = [ 0  Lb' ] * [ x1(1)    ]
00875   //      [ 0  0   ]   [ x1(2,mb) ]
00876   //
00877   //    = [ Lb'*x1(2,mb) ]
00878   //      [      0       ]
00879   //
00880   if( mb > 1 ) {
00881     V_MtV( &x2(1,mb-1), Lb(), BLAS_Cpp::trans, x1(2,mb) );
00882   }
00883   x2(mb) = 0.0;
00884 
00885   // x2 += -y2
00886   Vp_StV( &x2, -1.0, y2 );
00887 
00888   // x2 = inv(D) * x2
00889   {
00890     DVectorSlice::const_iterator
00891       d_itr = STY_.diag(0).begin();
00892     DVectorSlice::iterator
00893       x2_itr = x2.begin();
00894     while( x2_itr != x2.end() )
00895       *x2_itr++ /= *d_itr++;
00896   }
00897 }
00898 
00899 void MatrixSymPosDefLBFGS::assert_initialized() const
00900 {
00901   if(!n_)
00902     throw std::logic_error( "MatrixSymPosDefLBFGS::assert_initialized() : "
00903       "Error, matrix not initialized" );
00904 }
00905 
00906 } // end namespace ConstrainedOptPack 
00907 
00908 namespace {
00909 
00910 void comp_Cb(
00911   const DMatrixSlice& Lb, const DVectorSlice& Db_diag
00912   ,DMatrixSlice* Cb
00913   )
00914 {
00915   // Lb * inv(Db) * Lb =
00916   //
00917   // [ l(1,1)           ]   [ dd(1)         ]   [ l(1,1)  l(2,1)  ... l(p,1)  ]
00918   // [ l(2,1) l(2,2)        ]   [   dd(2)     ]   [     l(2,2)  ... l(p,2)  ]
00919   // [ .    .        .      ] * [     .     ] * [         . .   ]
00920   // [ l(p,1) l(p,2)  ... l(p,p)  ]   [       dd(p) ]   [           l(p,p)  ]
00921   //
00922   //
00923   // [ l(1,1)*dd(1)*l(1,1)  l(1,1)*dd(1)*l(2,1)             ...   l(1,1)*dd(1)*l(p,1)             ]
00924   // [ 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) ]   
00925   // [ .            .                     ...
00926   // [ symmetric        symmetric                 ...   sum( l(p,i)*dd(i)*l(p,i), i=1,..,p )    ]
00927   //
00928   // Therefore we can express the upper triangular elemetns of Cb as:
00929   //
00930   // Cb(i,j) = sum( l(i,k)*dd(k)*l(j,k), k = 1,..,i )
00931 
00932   typedef DenseLinAlgPack::size_type size_type;
00933   typedef DenseLinAlgPack::value_type value_type;
00934 
00935   TEST_FOR_EXCEPT( !(  Lb.rows() == Cb->rows() && Cb->rows() == Db_diag.dim()  ) ); // only a local error!
00936 
00937   const size_type p = Db_diag.dim();
00938 
00939   for( size_type i = 1; i <= p; ++i ) {
00940     for( size_type j = i; j <= p; ++j ) {
00941       value_type &c = (*Cb)(i,j) = 0.0;
00942       for( size_type k = 1; k <= i; ++k ) {
00943         c += Lb(i,k) * Lb(j,k) / Db_diag(k);
00944       }
00945     }
00946   }
00947 
00948   // ToDo: Make the above operation more efficent if needed! (i.e. write
00949   // it in fortran or something?).
00950 }
00951 
00952 } // end namespace

Generated on Tue Jul 13 09:30:51 2010 for MOOCHO (Single Doxygen Collection) by  doxygen 1.4.7