[Lme4-commits] r1474 - pkg/lme4Eigen/src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Dec 12 23:02:35 CET 2011
Author: dmbates
Date: 2011-12-12 23:02:35 +0100 (Mon, 12 Dec 2011)
New Revision: 1474
Added:
pkg/lme4Eigen/src/optimizer.cpp
pkg/lme4Eigen/src/optimizer.h
Modified:
pkg/lme4Eigen/src/glmFamily.h
Log:
Added simple optimizers (Golden search for bounded 1-dimensional and Nelder-Mead for multidimensional).
Modified: pkg/lme4Eigen/src/glmFamily.h
===================================================================
--- pkg/lme4Eigen/src/glmFamily.h 2011-12-07 23:25:59 UTC (rev 1473)
+++ pkg/lme4Eigen/src/glmFamily.h 2011-12-12 22:02:35 UTC (rev 1474)
@@ -88,6 +88,7 @@
qq = -std::exp(qq);
return lower_tail ? -expm1(qq) : std::exp(qq);
}
+
static inline double cloglogLinkInv(const double& x) {
return pgumbel2(x, 0., 1., 1);
}
Added: pkg/lme4Eigen/src/optimizer.cpp
===================================================================
--- pkg/lme4Eigen/src/optimizer.cpp (rev 0)
+++ pkg/lme4Eigen/src/optimizer.cpp 2011-12-12 22:02:35 UTC (rev 1474)
@@ -0,0 +1,312 @@
+//
+// Nelder_Mead.cpp: implementation of Nelder-Mead optimization algorithm
+//
+// Based on the files nldrmd.h, nldrmd.c, stop.c and nlopt-util.h from NLopt 2.2.4
+// Steven G. Johnson, The NLopt nonlinear-optimization package,
+// http://ab-initio.mit.edu/nlopt
+//
+// Original implementation Copyright (C) 2007-2011 Massachusetts Institute of Technology
+// Modifications Copyright (C) 2011 Douglas Bates, Martin Maechler and Ben Bolker
+//
+// This file is part of lme4Eigen.
+
+#include "optimizer.h"
+
+namespace optimizer {
+ using std::invalid_argument;
+ using std::runtime_error;
+
+ typedef VectorXd::Scalar Scalar;
+ typedef VectorXd::Index Index;
+
+/**
+ * Determine if two values are approximately equal relative to floating-point precision
+ *
+ * @param a first value to compare
+ * @param b second value to compare
+ *
+ * @return true if a and b are approximately equal else false
+ */
+ static bool close(const Scalar& a, const Scalar& b) {
+ return (std::abs(a - b) <= 1e-13 * (std::abs(a) + std::abs(b)));
+ }
+
+/**
+ *
+ *
+ * @param lb lower bounds
+ * @param ub upper bounds
+ * @param xstep initial step sizes
+ * @param x initial parameter vector
+ * @param f value of function at initial parameter vector
+ */
+ Nelder_Mead::Nelder_Mead(const VectorXd& lb, const VectorXd& ub,
+ const VectorXd& xstep, const VectorXd& x,
+ const nl_stop& stp)
+ : d_lb( lb),
+ d_ub( ub),
+ d_xstep( xstep),
+ d_x( x),
+ d_n( x.size()),
+ d_pts( d_n, d_n + 1),
+ d_vals( d_n + 1),
+ d_c( d_n),
+ d_xcur( d_n),
+ d_xeval( x),
+ d_minf( std::numeric_limits<Scalar>::infinity()),
+ d_stage( nm_restart),
+ d_stop( stp)
+ {
+ if (!d_n || d_lb.size() != d_n ||
+ d_ub.size() != d_n || d_xstep.size() != d_n)
+ throw invalid_argument("dimension mismatch");
+ if (((d_x - d_lb).array() < 0).any() || ((d_ub - d_x).array() < 0).any())
+ throw std::invalid_argument("initial x is not a feasible point");
+ d_stop.resetEvals();
+ init_pos = 0;
+ d_pts = d_x.replicate(1, d_n + 1);
+ for (Index i = 0; i < d_n; ++i) { // generate and check the initial positions
+ Index j(i + 1);
+ d_pts(i, j) += d_xstep[i];
+ if (d_pts(i, j) > d_ub[i]) {
+ d_pts(i, j) = (d_ub[i] - d_x[i] > std::abs(d_xstep[i]) * 0.1) ? d_ub[i]
+ // ub is too close to pt, go in other direction
+ : d_x[i] - std::abs(d_xstep[i]);
+ }
+ if (d_pts(i,j) < d_lb[i]) {
+ if (d_x[i] - d_lb[i] > std::abs(d_xstep[i]) * 0.1) d_pts(i,j) = d_lb[i];
+ else { // lb is too close to pt, go in other direction
+ d_pts(i,j) = d_x[i] + std::abs(d_xstep[i]);
+ if (d_pts(i, j) > d_ub[i]) // go towards farther of lb, ub */
+ d_pts(i, j) = 0.5 * ((d_ub[i] - d_x[i] > d_x[i] - d_lb[i] ?
+ d_ub[i] : d_lb[i]) + d_x[i]);
+ }
+ }
+ if (close(d_pts(i,j), d_x[i]))
+ throw std::invalid_argument("cannot generate feasible simplex");
+ }
+ }
+
+/**
+ * Install the function value at d_xeval
+ *
+ * @param f value of function at d_xeval
+ *
+ * @return status
+ */
+ nm_status Nelder_Mead::newf(const Scalar& f) {
+ Rcpp::Rcout << "f = " << f
+ << " at " << d_xeval.adjoint() << std::endl;
+ d_stop.incrEvals();
+ if (d_stop.forced()) return nm_forced;
+ if (f < d_minf) {
+ d_minf = f;
+ d_x = d_xeval; // save the value generating current minimum
+ if (d_minf < d_stop.minfMax()) return nm_minf_max;
+ }
+ if (d_stop.evals()) return nm_evals;
+ if (init_pos <= d_n) return init(f);
+ switch (d_stage) {
+ case nm_restart: return restart(f);
+ case nm_postreflect: return postreflect(f);
+ case nm_postexpand: return postexpand(f);
+ case nm_postcontract: return postcontract(f);
+ }
+ return nm_active; // -Wall
+ }
+
+/**
+ * Initialization of d_vals from the positions in d_pts;
+ *
+ * @param f function value
+ *
+ * @return status
+ */
+ nm_status Nelder_Mead::init(const Scalar& f) {
+ if (init_pos > d_n) throw std::runtime_error("init called after n evaluations");
+ d_vals[init_pos++] = f;
+ if (init_pos > d_n) return restart(f);
+ d_xeval = d_pts.col(init_pos);
+ return nm_active;
+ }
+
+/**
+ * Recompute the high/low function values (d_fh and d_fl) and indices
+ * (d_ih and d_il) plus the centroid of the n-1 simplex opposite the high
+ * vertex. Check if the simplex has collapsed. If not, attempt a reflection.
+ *
+ * @param f function value
+ *
+ * @return status
+ */
+ nm_status Nelder_Mead::restart(const Scalar& f) {
+ d_fl = d_vals.minCoeff(&d_il);
+ d_fh = d_vals.maxCoeff(&d_ih);
+ Rcpp::Rcout << "restart: neval = " << d_stop.ev()
+ << ", d_fl = " << d_fl
+ << ", d_fh = " << d_fh
+ << ", d_il = " << d_il
+ << ", d_ih = " << d_ih << std::endl;
+ d_c = (d_pts.rowwise().sum() - d_pts.col(d_ih)) / d_n; // compute centroid
+ Rcpp::Rcout << "d_c: " << d_c.adjoint() << std::endl;
+
+ // Check if the simplex has gotten to be too small. First
+ // calculate the coefficient-wise maximum abs deviation for
+ // the centroid.
+ d_xcur = (d_pts.colwise() - d_c).array().abs().rowwise().maxCoeff();
+ if (d_stop.x(VectorXd::Zero(d_n), d_xcur)) return nm_xcvg;
+
+ if (!reflectpt(d_xcur, d_c, alpha, d_pts.col(d_ih))) return nm_xcvg;
+ Rcpp::Rcout << "d_xcur (after reflectpt): " << d_xcur.adjoint() << std::endl;
+ d_xeval = d_xcur;
+ d_stage = nm_postreflect;
+ return nm_active;
+ }
+
+ nm_status Nelder_Mead::postreflect(const Scalar& f) {
+ Rcpp::Rcout << "postreflect: ";
+ if (f < d_fl) { // new best point, try to expand
+ if (!reflectpt(d_xeval, d_c, gamm, d_pts.col(d_ih))) return nm_xcvg;
+ Rcpp::Rcout << "New best point" << std::endl;
+ d_stage = nm_postexpand;
+ f_old = f;
+ return nm_active;
+ }
+ if (f < d_fh) { // accept new point
+ Rcpp::Rcout << "Accept new point" << std::endl;
+ d_vals[d_ih] = f;
+ d_pts.col(d_ih) = d_xeval;
+ return restart(f);
+ }
+ // new worst point, contract
+ Rcpp::Rcout << "New worst point" << std::endl;
+ if (!reflectpt(d_xcur, d_c, d_fh <= f ? -beta : beta, d_pts.col(d_ih))) return nm_xcvg;
+ f_old = f;
+ d_xeval = d_xcur;
+ d_stage = nm_postcontract;
+ return nm_active;
+ }
+
+ nm_status Nelder_Mead::postexpand(const Scalar& f) {
+ if (f < d_vals[d_ih]) { // expanding improved
+ Rcpp::Rcout << "successful expand" << std::endl;
+ d_pts.col(d_ih) = d_xeval;
+ d_vals[d_ih] = f;
+ } else {
+ Rcpp::Rcout << "unsuccessful expand" << std::endl;
+ d_pts.col(d_ih) = d_xcur;
+ d_vals[d_ih] = f_old;
+ }
+ return restart(f);
+ }
+
+ nm_status Nelder_Mead::postcontract(const Scalar& f) {
+ if (f < f_old && f < d_fh) {
+ Rcpp::Rcout << "successful contraction:" << std::endl;
+ d_pts.col(d_ih) = d_xeval;
+ d_vals[d_ih] = f;
+ return restart(f);
+ }
+ Rcpp::Rcout << "unsuccessful contraction, shrink simplex" << std::endl;
+ for (Index i = 0; i <= d_n; ++i) {
+ if (i != d_il) {
+ if (!reflectpt(d_xeval, d_pts.col(d_il), -delta, d_pts.col(i))) return nm_xcvg;
+ d_pts.col(i) = d_xeval;
+ }
+ }
+ init_pos = 0;
+ d_xeval = d_pts.col(0);
+ return nm_active;
+ }
+
+/* Perform the reflection xnew = c + scale * (c - xold),
+ returning 0 if xnew == c or xnew == xold (coincident points), 1 otherwise.
+
+ The reflected point xnew is "pinned" to the lower and upper bounds
+ (lb and ub), as suggested by J. A. Richardson and J. L. Kuester,
+ "The complex method for constrained optimization," Commun. ACM
+ 16(8), 487-489 (1973). This is probably a suboptimal way to handle
+ bound constraints, but I don't know a better way. The main danger
+ with this is that the simplex might collapse into a
+ lower-dimensional hyperplane; this danger can be ameliorated by
+ restarting (as in subplex), however. */
+ bool Nelder_Mead::reflectpt(VectorXd& xnew, const VectorXd& c,
+ const Scalar& scale, const VectorXd& xold) {
+ xnew = c + scale * (c - xold);
+ bool equalc = true, equalold = true;
+ for (Index i = 0; i < d_n; ++i) {
+ Scalar newx = std::min(std::max(xnew[i], d_lb[i]), d_ub[i]);
+ equalc = equalc && close(newx, d_c[i]);
+ equalold = equalold && close(newx, xold[i]);
+ xnew[i] = newx;
+ }
+ return !(equalc || equalold);
+ }
+
+ nl_stop::nl_stop(const VectorXd& xtol)
+ : xtol_abs( xtol),
+ maxeval( 300),
+ minf_max( std::numeric_limits<Scalar>::min()),
+ ftol_rel( 1e-15),
+ xtol_rel( 1e-7) {
+ }
+
+ bool nl_stop::x(const VectorXd& x, const VectorXd& oldx) const {
+ for (Index i = 0; i < x.size(); ++i)
+ if (!relstop(oldx[i], x[i], xtol_rel, xtol_abs[i])) return false;
+ return true;
+ }
+
+ bool nl_stop::dx(const VectorXd& x, const VectorXd& dx) const {
+ for (Index i = 0; i < x.size(); ++i)
+ if (!relstop(x[i] - dx[i], x[i], xtol_rel, xtol_abs[i])) return false;
+ return true;
+ }
+
+ bool nl_stop::xs(const VectorXd& xs, const VectorXd& oldxs,
+ const VectorXd& scale_min, const VectorXd& scale_max) const {
+ for (Index i = 0; i < xs.size(); ++i)
+ if (relstop(sc(oldxs[i], scale_min[i], scale_max[i]),
+ sc(xs[i], scale_min[i], scale_max[i]),
+ xtol_rel, xtol_abs[i]))
+ return true;
+ return false;
+ }
+
+ Golden::Golden(const Scalar& lower, const Scalar& upper)
+ : d_lower(lower), d_upper(upper) {
+ if (lower >= upper)
+ throw invalid_argument("lower >= upper");
+ d_invratio = 2./(1. + std::sqrt(5.));
+ double range = upper - lower;
+ d_x[0] = lower + range * (1. - d_invratio);
+ d_x[1] = lower + range * d_invratio;
+ d_init = true;
+ d_ll = true;
+ }
+
+ void Golden::newf(const Scalar& fv) {
+ Rcpp::Rcout << "f = " << fv << " at x = " << xeval() << std::endl;
+ d_f[d_ll ? 0 : 1] = fv;
+ if (d_init) {
+ d_init = false;
+ d_ll = false;
+ return;
+ }
+ if (d_f[0] > d_f[1]) { // discard left portion of interval
+ d_lower = d_x[0];
+ d_x[0] = d_x[1];
+ d_f[0] = d_f[1];
+ d_x[1] = d_lower + (d_upper - d_lower) * d_invratio;
+ d_ll = false;
+ } else {
+ d_upper = d_x[1];
+ d_x[1] = d_x[0];
+ d_f[1] = d_f[0];
+ d_x[0] = d_lower + (d_upper - d_lower) * (1 - d_invratio);
+ d_ll = true;
+ }
+ }
+
+}
+
Added: pkg/lme4Eigen/src/optimizer.h
===================================================================
--- pkg/lme4Eigen/src/optimizer.h (rev 0)
+++ pkg/lme4Eigen/src/optimizer.h 2011-12-12 22:02:35 UTC (rev 1474)
@@ -0,0 +1,139 @@
+// -*- mode: C++; c-indent-level: 4; c-basic-offset: 4; tab-width: 8 -*-
+//
+// Nelder_Mead.h: NLopt's Nelder-Mead optimizer, modified to use Eigen
+//
+// Copyright (C) 2011 Douglas Bates, Martin Maechler and Ben Bolker
+//
+// This file is part of lme4Eigen.
+
+#ifndef LME4_NELDER_MEAD_H
+#define LME4_NELDER_MEAD_H
+
+#include <RcppEigen.h>
+
+namespace optimizer {
+ using Eigen::MatrixXd;
+ using Eigen::VectorXd;
+ using Eigen::VectorXi;
+
+ typedef VectorXd::Scalar Scalar;
+ typedef VectorXd::Index Index;
+
+ class nl_stop {
+ private: // utilities
+ bool relstop(const Scalar& vold, const Scalar& vnew,
+ const Scalar& reltol, const Scalar& abstol) const;
+ Scalar sc(const Scalar& x, const Scalar& smin, const Scalar& smax) const {
+ return smin + x * (smax - smin);
+ }
+ protected:
+ const VectorXd xtol_abs;
+ unsigned n, nevals, maxeval;
+ Scalar minf_max, ftol_rel, ftol_abs, xtol_rel;
+ bool force_stop;
+ public:
+ nl_stop(const VectorXd&); // constructor
+
+ void resetEvals() {nevals = 0;} // setters
+ void incrEvals() {nevals++;}
+ void setMinf_max(const Scalar& mm) {minf_max = mm;}
+ void setFtol_rel(const Scalar& ftr) {ftol_rel = ftr;}
+ void setFtol_abs(const Scalar& fta) {ftol_abs = fta;}
+ void setForce_stop(const bool& stp) {force_stop = stp;}
+
+ bool f(const Scalar& f, const Scalar& oldf) const { // convergence checking
+ return (f <= minf_max || ftol(f, oldf));
+ }
+ bool ftol(const Scalar& f, const Scalar& oldf) const {
+ return relstop(oldf, f, ftol_rel, ftol_abs);
+ }
+ bool x(const VectorXd& x, const VectorXd& oldx) const;
+ bool dx(const VectorXd& x, const VectorXd& dx) const;
+ bool xs(const VectorXd& xs, const VectorXd& oldxs,
+ const VectorXd& scale_min, const VectorXd& scale_max) const;
+ bool evals() const {return maxeval > 0 && nevals > maxeval;}
+ bool forced() const {return force_stop;}
+
+ int ev() const {return nevals;}
+
+ Scalar minfMax() const {return minf_max;}
+ };
+
+ inline bool nl_stop::relstop(const Scalar& vold, const Scalar& vnew,
+ const Scalar& reltol, const Scalar& abstol) const {
+ if (std::abs(vold) == std::numeric_limits<Scalar>::infinity()) return false;
+ return std::abs(vnew - vold) < abstol
+ || std::abs(vnew - vold) < reltol * (std::abs(vnew) + std::abs(vold)) * 0.5
+ || (reltol > 0 && vnew == vold);
+ }
+
+ enum nm_status {nm_active, nm_x0notfeasible, nm_nofeasible, nm_forced, nm_minf_max,
+ nm_evals, nm_fcvg, nm_xcvg};
+
+ enum nm_stage {nm_restart, nm_postreflect, nm_postexpand, nm_postcontract};
+
+ /* heuristic "strategy" constants: */
+ static const double alpha = 1, beta = 0.5, gamm = 2, delta = 0.5;
+
+ class Nelder_Mead {
+ private:
+ Scalar f_old;
+ Index init_pos;
+ nm_status init(const Scalar&);
+ nm_status restart(const Scalar&);
+ bool reflectpt(VectorXd&, const VectorXd&, const Scalar&, const VectorXd&);
+ nm_status postreflect(const Scalar&);
+ nm_status postexpand(const Scalar&);
+ nm_status postcontract(const Scalar&);
+ protected:
+ const VectorXd d_lb; /*<< lower bounds */
+ const VectorXd d_ub; /*<< upper bounds */
+ const VectorXd d_xstep; /*<< initial step sizes */
+ VectorXd d_x; /*<< initial value and optimum */
+ Index d_ih; /**< index in d_vals of largest value */
+ Index d_il; /**< index in d_vals of smallest value */
+ Index d_n; /**< size of parameter vector */
+ MatrixXd d_pts; /*<< points */
+ VectorXd d_vals; /*<< function values */
+ VectorXd d_c; /*<< centroid */
+ VectorXd d_xcur; /*<< current x */
+ VectorXd d_xeval; /*<< x at which next evaluation is requested */
+ Scalar d_fl, d_fh, d_minf;
+ nm_status d_stat;
+ nm_stage d_stage;
+ nl_stop d_stop;
+ public:
+ Nelder_Mead(const VectorXd&, const VectorXd&, const VectorXd&, const VectorXd&,
+ const nl_stop&);
+
+ const MatrixXd& pts() const {return d_pts;}
+
+ const VectorXd& lb() const {return d_lb;}
+ const VectorXd& ub() const {return d_ub;}
+ const VectorXd& vals() const {return d_vals;}
+ const VectorXd& x() const {return d_x;}
+ const VectorXd& xstep() const {return d_xstep;}
+ const VectorXd& xeval() const {return d_xeval;}
+ Index ih() const {return d_ih;}
+ Index il() const {return d_il;}
+ Index evals() const {return d_stop.ev();}
+
+ Scalar minf() const {return d_minf;}
+ nm_status newf(const Scalar&);
+ };
+
+ class Golden {
+ protected:
+ Scalar d_invratio, d_lower, d_upper;
+ Eigen::Vector2d d_x, d_f;
+ bool d_init, d_ll;
+ public:
+ Golden(const Scalar&, const Scalar&);
+ void newf(const Scalar&);
+ Scalar xeval() const {return d_x[d_ll ? 0 : 1];}
+ Scalar value() const {return d_f[0];}
+ Scalar xpos() const {return d_x[0];}
+ };
+}
+
+#endif // LME4_NELDER_MEAD_H
More information about the Lme4-commits
mailing list