[Solution] Optimize function in Rcpp

[Xiurui_Zhu]

Recently I read about a topic regarding implementing stats::optimize() in C++: Optimize function in Rcpp. By coincidence, I just wrote a C++ version of this function, since calling R functions with Rcpp in C++ code is not a good practice and gains litter performace improbement. I woule like to share this in case anyone is in need.

Generally my script consists of three files:

• optimize.h: a header file defining the optimize() function, in a form close to stats::optimize()and an Optim class taking care of dot arguments in stats::optimize(). You do not need to apdapt this code.

• optimize.cpp: source file with a function adapted from stats C code: Brend_fmin(), the optimize() function and other helper (internal) functions. You do not need to apdapt this code.

• optim_test.cpp: source file with unit tests to demonstrate the use of defined functions. You may create your own subclass, which takes care of the function to optimize and the parameters for the function (passed on from dot arguments in R).

This is optimized.h:

#ifndef OPTIMIZE_R
#define OPTIMIZE_R

#include <iostream>
#include <cmath>
#include <cfloat>

// Optim class: virtual class to find max/min of univariate function in R manner
// Usually a subclass is used to define a substantiated function
// Member function (public): value, evaluate function with doulbe x
// Other parameters may be added into the subclass as members
class Optim
{
public:
	virtual double value(double x) = 0;
	virtual ~Optim() {}
};

// Define optimize function
double optimize(Optim* optim, double lower, double upper, bool maximum, double tol);

#endif

This is optimize.cpp:

#include "optimize.h"
using namespace std;

// This function is copied from R: stats/src/optimize.c
// Add argument maximum and its use in function pointer *f
// Define Brent optimization
double Brent_fmin(double ax, double bx, double (*f)(double, void *, bool),
		  void *info, bool maximum, double tol)
{
    /*  c is the squared inverse of the golden ratio */
    const double c = (3. - sqrt(5.)) * .5;

    /* Local variables */
    double a, b, d, e, p, q, r, u, v, w, x;
    double t2, fu, fv, fw, fx, xm, eps, tol1, tol3;

/*  eps is approximately the square root of the relative machine precision. */
    eps = DBL_EPSILON;
    tol1 = eps + 1.;/* the smallest 1.000... > 1 */
    eps = sqrt(eps);

    a = ax;
    b = bx;
    v = a + c * (b - a);
    w = v;
    x = v;

    d = 0.;/* -Wall */
    e = 0.;
    fx = (*f)(x, info, maximum);
    fv = fx;
    fw = fx;
    tol3 = tol / 3.;

/*  main loop starts here ----------------------------------- */

    for(;;) {
	xm = (a + b) * .5;
	tol1 = eps * fabs(x) + tol3;
	t2 = tol1 * 2.;

	/* check stopping criterion */

	if (fabs(x - xm) <= t2 - (b - a) * .5) break;
	p = 0.;
	q = 0.;
	r = 0.;
	if (fabs(e) > tol1) { /* fit parabola */

	    r = (x - w) * (fx - fv);
	    q = (x - v) * (fx - fw);
	    p = (x - v) * q - (x - w) * r;
	    q = (q - r) * 2.;
	    if (q > 0.) p = -p; else q = -q;
	    r = e;
	    e = d;
	}

	if (fabs(p) >= fabs(q * .5 * r) ||
	    p <= q * (a - x) || p >= q * (b - x)) { /* a golden-section step */

	    if (x < xm) e = b - x; else e = a - x;
	    d = c * e;
	}
	else { /* a parabolic-interpolation step */

	    d = p / q;
	    u = x + d;

	    /* f must not be evaluated too close to ax or bx */

	    if (u - a < t2 || b - u < t2) {
		d = tol1;
		if (x >= xm) d = -d;
	    }
	}

	/* f must not be evaluated too close to x */

	if (fabs(d) >= tol1)
	    u = x + d;
	else if (d > 0.)
	    u = x + tol1;
	else
	    u = x - tol1;

	fu = (*f)(u, info, maximum);

	/*  update  a, b, v, w, and x */

	if (fu <= fx) {
	    if (u < x) b = x; else a = x;
	    v = w;    w = x;   x = u;
	    fv = fw; fw = fx; fx = fu;
	} else {
	    if (u < x) a = u; else b = u;
	    if (fu <= fw || w == x) {
		v = w; fv = fw;
		w = u; fw = fu;
	    } else if (fu <= fv || v == x || v == w) {
		v = u; fv = fu;
	    }
	}
    }
    /* end of main loop */

    return x;
} // Brent_fmin()

// Optim's value function: wrapper around Optim::value
// This function is used in Brent_fmin as function pointer
double optim_value(double x, Optim* optim, bool maximum) {
	double out = (*optim).value(x);
	if (maximum == true) {
		out = -out;
	}
	return out;
}

// Optimize function: finding minimun of class-based univariate function
// optim: an object inheriting from Optim, with double value(double) member function
// Other parameters can be stored in members of Optim class
double optimize(Optim* optim, double lower, double upper, bool maximum = false,
				double tol = pow(DBL_EPSILON, 0.25)) {
	return Brent_fmin(lower, upper, (double (*)(double, void*, bool)) optim_value, optim,
		maximum, tol);
}

This is optim_test.cpp:

// Unit test
#include "optimize.h"
using namespace std;

// Define a class of hyperbolic function: f(x) = a * x^2 + b * x + c
class Parabol: public Optim
{
private:
	const double a;
	const double b;
	const double c;
public:
	Parabol(double a_, double b_, double c_) : a(a_), b(b_), c(c_) {}
	double value(double x) {
		return a * pow(x, 2) + b * x + c;
	}
};

// Main function to minimize a hyperbolic function
int main() {
	Parabol parabol1(1, -5, 3);
	double x_min = optimize(&parabol1, 0, 5, false, 1e-3);
	cout << x_min << endl;
	Parabol parabol2(-1, -5, 3);
	double x_max = optimize(&parabol2, -5, 0, true, 1e-3);
	cout << x_max << endl;
	return 0;
}

These files have been compiled with Visual C++ 6.0 and the output from main() is:

2.5
-2.5
Press any key to continue