optimize.h

使用R语言的马尔科夫链蒙特卡洛模拟（MCMC）源代码程序。

    T alpha_last = (T) 0.0;
    T alpha_cur = (T) 1.0;
    T alpha_max = (T) 10.0;
    T c1 = (T) 1e-4;
    T c2 = (T) 0.5;
    unsigned int max_iter = 50;
    T fgradalpha0 = gradfdifls(fun, (T) 0, theta, p);

    for (unsigned int i = 0; i < max_iter; ++i) {
      T phi_cur = fun(theta + alpha_cur * p);
      T phi_last = fun(theta + alpha_last * p);
     
      if ((phi_cur > (fun(theta) + c1 * alpha_cur * fgradalpha0))
          || ((phi_cur >= phi_last) && (i > 0))) {
        T alphastar = zoom(fun, alpha_last, alpha_cur, theta, p);
        return alphastar;
      }

      T fgradalpha_cur = gradfdifls(fun, alpha_cur, theta, p);
      if (std::fabs(fgradalpha_cur) <= -1 * c2 * fgradalpha0)
        return alpha_cur;

      if ( fgradalpha_cur >= (T) 0.0) {
        T alphastar = zoom(fun, alpha_cur, alpha_last, theta, p);
        return alphastar;
      }
      
      alpha_last = alpha_cur;
      // runif stuff below is probably not correc KQ 12/08/2006
      // I think it should work now DBP 01/02/2007
      alpha_cur = runif() * (alpha_max - alpha_cur) + alpha_cur;
    }

    return 0.001;
  }

   /*! \brief Find minimum of a function once bracketed.
    *
    * This function finds the minimum of a function, once bracketed.
    *
    * \param fun The function to minimize.  This function should take
    * one Matrix (vector) argument of type T and return a single value
    * of type T.
    * \param alpha_lo The lower bracket.
    * \param alpha_hi The upper bracket.
    * \param theta A column vector of parameter values that anchor the
    * 1-dimensional function.
    * \param p A direction vector that creates the 1-dimensional
    *
    * \see linesearch1(FUNCTOR fun, const Matrix<T,PO1,PS1>& theta, const Matrix<T,PO2,PS2>& p)
    * \see linesearch2(FUNCTOR fun, const Matrix<T,PO1,PS1>& theta, const Matrix<T,PO2,PS2>& p, rng<RNGTYPE>& runif)
    * \see BFGS(FUNCTOR fun, const Matrix<T,PO,PS>& theta, rng<RNGTYPE>& runif, unsigned int maxit, T tolerance, bool trace = false)
    *
    * \throw scythe_dimension_error (Level 1)
    *
    * \note
    * Users will typically wish to implement \a fun in terms of a
    * functor.  Using a functor provides a generic way in which to
    * evaluate functions with more than one parameter.  Furthermore,
    * although one can pass a function pointer to this routine,
    * the compiler cannot inline and fully optimize code
    * referenced by function pointers.
    * 
    */
  template <typename T, matrix_order PO1, matrix_style PS1,
            matrix_order PO2, matrix_style PS2, typename FUNCTOR>
  T zoom (FUNCTOR fun, T alpha_lo, T alpha_hi,
          const Matrix<T,PO1,PS1>& theta, const Matrix<T,PO2,PS2>& p)
  {
    SCYTHE_CHECK_10(! theta.isColVector(), scythe_dimension_error,
        "Theta not column vector");
    SCYTHE_CHECK_10(! p.isColVector(), scythe_dimension_error,
        "p not column vector");

    T alpha_j = (alpha_lo + alpha_hi) / 2.0;
    T phi_0 = fun(theta);
    T c1 = (T) 1e-4;
    T c2 = (T) 0.5;
    T fgrad0 = gradfdifls(fun, (T) 0, theta, p);

    unsigned int count = 0;
    unsigned int maxit = 20;
    while(count < maxit) {
      T phi_j = fun(theta + alpha_j * p);
      T phi_lo = fun(theta + alpha_lo * p);
     
      if ((phi_j > (phi_0 + c1 * alpha_j * fgrad0))
          || (phi_j >= phi_lo)){
        alpha_hi = alpha_j;
      } else {
        T fgradj = gradfdifls(fun, alpha_j, theta, p);
        if (std::fabs(fgradj) <= -1 * c2 * fgrad0){ 
          return alpha_j;
        }
        if ( fgradj * (alpha_hi - alpha_lo) >= 0){
          alpha_hi = alpha_lo;
        }
        alpha_lo = alpha_j;
      }
      ++count;
    }
   
    return alpha_j;
  }


   /*! \brief Find function minimum using the BFGS algorithm.
    *
    * Numerically find the minimum of a function using the BFGS
    * algorithm.
    *
    * \param fun The function to minimize.  This function should take
    * one Matrix (vector) argument of type T and return a single value
    * of type T.
    * \param theta A column vector of parameter values that anchor the
    * 1-dimensional function.
    * \param runif A random uniform number generator function object
    * (an object that returns a random uniform variate on (0,1) when
    * its () operator is invoked).
    * \param maxit The maximum number of iterations.
    * \param tolerance The convergence tolerance.
    * \param trace Boolean value determining whether BFGS should print 
    *              to stdout (defaults to false).
    *
    * \see linesearch1(FUNCTOR fun, const Matrix<T,PO1,PS1>& theta, const Matrix<T,PO2,PS2>& p)
    * \see linesearch2(FUNCTOR fun, const Matrix<T,PO1,PS1>& theta, const Matrix<T,PO2,PS2>& p, rng<RNGTYPE>& runif)
    * \see zoom(FUNCTOR fun, T alpha_lo, T alpha_hi, const Matrix<T,PO1,PS1>& theta, const Matrix<T,PO2,PS2>& p)
    *
    * \throw scythe_dimension_error (Level 1)
    * \throw scythe_convergence_error (Level 0)
    *
    * \note
    * Users will typically wish to implement \a fun in terms of a
    * functor.  Using a functor provides a generic way in which to
    * evaluate functions with more than one parameter.  Furthermore,
    * although one can pass a function pointer to this routine,
    * the compiler cannot inline and fully optimize code
    * referenced by function pointers.
    */
  // there were 2 versions of linesearch1-- the latter was what we
  // had been calling linesearch2
  template <matrix_order RO, matrix_style RS, typename T, 
            matrix_order PO, matrix_style PS,
            typename FUNCTOR, typename RNGTYPE>
  Matrix<T,RO,RS>
  BFGS (FUNCTOR fun, const Matrix<T,PO,PS>& theta, rng<RNGTYPE>& runif, 
        unsigned int maxit, T tolerance, bool trace = false)
  {
    SCYTHE_CHECK_10(! theta.isColVector(), scythe_dimension_error,
        "Theta not column vector");

    unsigned int n = theta.size();

    // H is initial inverse hessian
    Matrix<T,RO> H = inv(hesscdif(fun, theta));

    // gradient at starting values
    Matrix<T,RO> fgrad = gradfdif(fun, theta);
    Matrix<T,RO> thetamin = theta;
    Matrix<T,RO> fgrad_new = fgrad;
    Matrix<T,RO> I = eye<T,RO>(n); 
    Matrix<T,RO> s;
    Matrix<T,RO> y;

    unsigned int count = 0;
    while( (t(fgrad_new)*fgrad_new)[0] > tolerance) {
      Matrix<T> p = -1.0 * H * fgrad;
      //std::cout << "initial H * fgrad = " << H * fgrad << "\n";
      //std::cout << "initial p = " << p << "\n";

      T alpha = linesearch2(fun, thetamin, p, runif);
      //T alpha = linesearch1(fun, thetamin, p);

      //std::cout << "after linesearch p = " << p << "\n";


      Matrix<T> thetamin_new = thetamin + alpha * p;
      fgrad_new = gradfdif(fun, thetamin_new);
      s = thetamin_new - thetamin;
      y = fgrad_new - fgrad;
      T rho = 1.0 / (t(y) * s)[0];
      H = (I - rho * s * t(y)) * H *(I - rho * y * t(s))
        + rho * s * t(s);

      thetamin = thetamin_new;
      fgrad = fgrad_new;
      ++count;

      if (trace){
	std::cout << "BFGS iteration = " << count << std::endl;
	std::cout << "thetamin = " << (t(thetamin)) ;
	std::cout << "gradient = " << (t(fgrad)) ;
	std::cout << "t(gradient) * gradient = " << 
	  (t(fgrad) * fgrad) ;
	std::cout << "function value = " << fun(thetamin) << 
	  std::endl << std::endl;
      }
      //std::cout << "Hessian = " << hesscdif(fun, theta) << "\n";
      //std::cout << "H = " << H << "\n";
      //std::cout << "alpha = " << alpha << std::endl;
      //std::cout << "p = " << p << "\n";
      //std::cout << "-1 * H * fgrad = " << -1.0 * H * fgrad << "\n";

      SCYTHE_CHECK(count > maxit, scythe_convergence_error,
          "Failed to converge.  Try better starting values");
    }
   
    return thetamin;
  }

  // Default template
  template <typename T, matrix_order O, matrix_style S,
            typename FUNCTOR, typename RNGTYPE>
  Matrix<T,O,Concrete>
  BFGS (FUNCTOR fun, const Matrix<T,O,S>& theta, rng<RNGTYPE>& runif,
        unsigned int maxit, T tolerance, bool trace = false)
  {
    return BFGS<O,Concrete> (fun, theta, runif, maxit, tolerance, trace);
  }

  
  /* Solves a system of n nonlinear equations in n unknowns of the form
   * fun(thetastar) = 0 for thetastar given the function, starting 
   * value theta, max number of iterations, and tolerance.
   * Uses Broyden's method.
   */
   /*! \brief Solve a system of nonlinear equations.
    *
    * Solves a system of n nonlinear equations in n unknowns of the form
    * \f$fun(\theta^*) = 0\f$ for \f$\theta^*\f$.
    *
    * \param fun The function to solve.  The function should both take
    * and return a Matrix of type T.
    * \param theta A column vector of parameter values at which to
    * start the solve procedure.
    * \param maxit The maximum number of iterations.
    * \param tolerance The convergence tolerance.
    *
    * \throw scythe_dimension_error (Level 1)
    * \throw scythe_convergence_error (Level 1)
    *
    * \note
    * Users will typically wish to implement \a fun in terms of a
    * functor.  Using a functor provides a generic way in which to
    * evaluate functions with more than one parameter.  Furthermore,
    * although one can pass a function pointer to this routine,
    * the compiler cannot inline and fully optimize code
    * referenced by function pointers.
    */
  template <matrix_order RO, matrix_style RS, typename T,
            matrix_order PO, matrix_style PS, typename FUNCTOR>
  Matrix<T,RO,RS>
  nls_broyden(FUNCTOR fun, const Matrix<T,PO,PS>& theta,
              unsigned int maxit = 5000, T tolerance = 1e-6)
  {
    SCYTHE_CHECK_10(! theta.isColVector(), scythe_dimension_error,
        "Theta not column vector");


    Matrix<T,RO> thetastar = theta;
    Matrix<T,RO> B = jacfdif(fun, thetastar);

    Matrix<T,RO> fthetastar;
    Matrix<T,RO> p;
    Matrix<T,RO> thetastar_new;
    Matrix<T,RO> fthetastar_new;
    Matrix<T,RO> s;
    Matrix<T,RO> y;

    for (unsigned int i = 0; i < maxit; ++i) {
      fthetastar = fun(thetastar);
      p = lu_solve(B, -1 * fthetastar);
      T alpha = (T) 1.0;
      thetastar_new = thetastar + alpha*p;
      fthetastar_new = fun(thetastar_new);
      s = thetastar_new - thetastar;
      y = fthetastar_new - fthetastar;
      B = B + ((y - B * s) * t(s)) / (t(s) * s);
      thetastar = thetastar_new;
      if (max(fabs(fthetastar_new)) < tolerance)
        return thetastar;
    }
 
    SCYTHE_THROW_10(scythe_convergence_error,  "Failed to converge.  Try better starting values or increase maxit");

    return thetastar;
  }

  // default template
  template <typename T, matrix_order O, matrix_style S,
            typename FUNCTOR>
  Matrix<T,O,Concrete>
  nls_broyden (FUNCTOR fun, const Matrix<T,O,S>& theta,
               unsigned int maxit = 5000, T tolerance = 1e-6)
  {
    return nls_broyden<O,Concrete>(fun, theta, maxit, tolerance);
  }

} // namespace scythe

#endif /* SCYTHE_OPTIMIZE_H */