apaam/mlpack
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

optim:GradientDescent

Browse Source
This commit is contained in:
tqlong
2010-07-14 21:37:40 +00:00
parent e4efc7eb33
commit 61d4fbfa53
6 changed files with 323 additions and 40 deletions
Download Patch File Download Diff File Expand all files Collapse all files
fastlib/trunk/contrib/tqlong/optimization/CMakeLists.txt
+1
View File
@@ -5,6 +5,7 @@ add_executable(optim_tqlong_test
test1.cpp
FunctionTemplate
NelderMead
GradientDescent
)
# link dependencies of executable
target_link_libraries(optim_tqlong_test
fastlib/trunk/contrib/tqlong/optimization/FunctionTemplate.cpp
+3 -1
View File
@@ -4,6 +4,7 @@
double optim::LengthEuclidianSquare::CalculateValue(const Vector &x)
{
DEBUG_ASSERT(x.length() == dim);
double s = 0;
for (int i = 0; i < x.length(); i++)
s += x[i]*x[i];
@@ -12,13 +13,14 @@ double optim::LengthEuclidianSquare::CalculateValue(const Vector &x)
void optim::LengthEuclidianSquare::CalculateGradient(const Vector &x, Vector &gradient)
{
DEBUG_ASSERT(x.length() == dim && gradient.length() == dim);
DEBUG_ASSERT(x.length() == gradient.length());
la::ScaleOverwrite(2.0, x, &gradient);
}
void optim::LengthEuclidianSquare::CalculateHessian(const Vector &x, Matrix &hessian)
{
DEBUG_ASSERT(x.length() == hessian.n_rows() && x.length() == hessian.n_cols());
DEBUG_ASSERT(x.length() == hessian.n_rows() && x.length() == hessian.n_cols() && x.length() == dim);
hessian.SetAll(0.0);
for (int i = 0; i < hessian.n_rows(); i++)
hessian.ref(i, i) = 2.0;
fastlib/trunk/contrib/tqlong/optimization/FunctionTemplate.h
+11 -1
View File
@@ -14,16 +14,26 @@ BEGIN_OPTIM_NAMESPACE;
/*************************************************************
* A function template, implement CalculateXXXXX methods
* 0-order smooth function: CalculateValue
* Function value : CalculateValue
* 1-order smooth function: CalculateValue, CalculateGradient
* 2-order smooth function: CalculateValue, CalculateGradient, CalculateHessian
* Required:
int dimension()
void Init(variable_type*)
variable_type (typedef)
*************************************************************/
class LengthEuclidianSquare {
int dim;
public:
typedef Vector variable_type; // required
int dimension() { return dim; } // required
void Init(Vector* x) { x->Init(dim); } // required
double CalculateValue(const Vector& x);
void CalculateGradient(const Vector& x, Vector& gradient);
void CalculateHessian(const Vector& x, Matrix& hessian);
public:
LengthEuclidianSquare(int dim_ = 2) : dim(dim_) {}
};
END_OPTIM_NAMESPACE;
fastlib/trunk/contrib/tqlong/optimization/GradientDescent.h
+234
View File
@@ -0,0 +1,234 @@
#ifndef GRADIENTDESCENT_H
#define GRADIENTDESCENT_H
#include <fastlib/fastlib.h>
#ifndef BEGIN_OPTIM_NAMESPACE
#define BEGIN_OPTIM_NAMESPACE namespace optim {
#endif
#ifndef END_OPTIM_NAMESPACE
#define END_OPTIM_NAMESPACE }
#endif
BEGIN_OPTIM_NAMESPACE;
/**********************************************************************
Implement 1-order Gradient Descent optimization method with Wolfe line search method
Function: CalculateValue(const variable_type&)
CalculateGradient(const variable_type&, variable_type& gradient)
Function::variable_type:
implement la::AddExpert(double, const variable_type&, variable_type*)
la::ScaleOverwrite(double, const variable_type&, variable_type*);
la::LengthEuclidean(const variable_type&)
la::Dot(const variable_type&, const variable_type&)
variable_type.CopyValues(const variable_type&)
Parameter:
General params
maxIter, rTol, aTol
Specific for Gradient Descent & Wolfe line search
c1, c2, beta
**********************************************************************/
template<typename Function> class GradientDescent
{
public:
typedef Function function_type;
typedef typename Function::variable_type variable_type;
typedef double* OptimizationParameters;
//maximum # of iterations : maxIter = (int) param[0];
//relative tolerance : rTol = param[1];
//absolute tolerance : aTol = param[2];
//c1 : Wolfe 1st condition : c1 = param[3];
//c2 : Wolfe 2nd condition : c2 = param[4];
//beta : scale parameter : beta = param[5];
struct HistoryRecord {
int iter;
int n_evals;
int n_grads;
double best_val;
double residual;
OT_DEF(HistoryRecord) {
OT_MY_OBJECT(iter);
OT_MY_OBJECT(n_evals);
OT_MY_OBJECT(n_grads);
OT_MY_OBJECT(best_val);
OT_MY_OBJECT(residual);
}
public:
HistoryRecord(int iter_, int n_evals_, int n_grads_, double best_val_, double residual_)
: iter(iter_), n_evals(n_evals_), n_grads(n_grads_), best_val(best_val_), residual(residual_) {}
};
protected:
static double default_gradient_descent_parameter[6]; // = {100, 0.001, 0.01, 1e-4, 0.9, 0.5};
public:
GradientDescent(function_type& f_, OptimizationParameters param_ = default_gradient_descent_parameter);
void setParam(OptimizationParameters param);
void setX0(const variable_type& x0_);
double optimize(variable_type& sol);
void printHistory();
ArrayList<HistoryRecord> history;
protected:
function_type& f;
variable_type x0;
// parameters
OptimizationParameters param;
int maxIter;
double aTol, rTol;
double c1, c2, beta;
// progress
int iter;
int n_evals;
int n_grads;
double best_val;
double residual;
double WolfeStep(const variable_type& x, double val_x, const variable_type& grad, const variable_type& p);
double CalculateValue(const variable_type& x);
void CalculateGradient(const variable_type& x, variable_type& grad);
void recordProgress();
void printProgress();
};
template<typename F>
double GradientDescent<F>::CalculateValue(const variable_type& x) {
double val = f.CalculateValue(x);
n_evals++;
return val;
}
template<typename F>
void GradientDescent<F>::CalculateGradient(const variable_type& x, variable_type& grad) {
f.CalculateGradient(x, grad);
n_grads++;
}
template<typename F>
void GradientDescent<F>::recordProgress() {
history.PushBackCopy(HistoryRecord(iter, n_evals, n_grads, best_val, residual));
printProgress();
}
template<typename F>
void GradientDescent<F>::printProgress() {
if (history.size() == 0) return;
printf("iter = %d n_evals = %d n_grads = %d best_val = %f residual = %f\n",
history.back().iter, history.back().n_evals, history.back().n_grads,
history.back().best_val, history.back().residual);
}
template<typename F>
void GradientDescent<F>::printHistory() {
printf("History:\n");
for (int i = 0; i < history.size(); i++) {
printf("iter = %d n_evals = %d n_grads = %d best_val = %f residual = %f\n",
history[i].iter, history[i].n_evals, history[i].n_grads, history[i].best_val, history[i].residual);
}
}
template<typename F>
double GradientDescent<F>::default_gradient_descent_parameter[6] = {100, 0.001, 0.01, 1e-4, 0.9, 0.5};
template<typename F>
GradientDescent<F>::GradientDescent(function_type &f_, OptimizationParameters param_)
: f(f_), param(param_) {
setParam(param);
f.Init(&x0);
history.Init();
}
template<typename F>
void GradientDescent<F>::setParam(OptimizationParameters param) {
maxIter = (int) param[0];
rTol = param[1];
aTol = param[2];
c1 = param[3];
c2 = param[4];
beta = param[5];
}
template<typename F>
void GradientDescent<F>::setX0(const variable_type& x0_) {
x0.CopyValues(x0_);
}
template<typename F>
double GradientDescent<F>::optimize(variable_type &sol) {
history.Clear();
iter = 0;
n_evals = 0;
n_grads = 0;
variable_type x; // the current search variable
variable_type grad, d; // the current search direction
f.Init(&x); // initialize search variable
f.Init(&d); // and direction
f.Init(&grad); // and direction
x.CopyValues(x0);
double val = CalculateValue(x);
CalculateGradient(x, grad); // Calculate gradient
double r0 = la::LengthEuclidean(grad);
sol.CopyValues(x0);
best_val = val;
residual = r0;
recordProgress();
for (iter = 1; iter < maxIter; iter++) {
// Calculate search direction: negative gradient
la::ScaleOverwrite(-1.0, grad, &d);
// Calculate step size by Wolfe's conditions
double lambda = WolfeStep(x, val, grad, d);
// Calculate new variable
la::AddExpert(lambda, d, &x);
// Update best value
val = CalculateValue(x);
if (val < best_val) {
best_val = val;
sol.CopyValues(x);
}
// Calculate gradient and check termination condition
CalculateGradient(x, grad);
residual = la::LengthEuclidean(grad);
recordProgress();
if (residual < rTol*r0+aTol) break;
}
return best_val;
}
template<typename F>
double GradientDescent<F>::WolfeStep(const variable_type& x, double val_x, const variable_type& grad, const variable_type& p) {
double lambda = 1.0/beta, val_xp;
double dot_grad_p;
variable_type x_p, grad_p;
f.Init(&x_p);
f.Init(&grad_p);
dot_grad_p = la::Dot(grad, p);
while (1) {
lambda *= beta;
if (lambda < 1e-10) {
printf("Line search results in a too small step size, try increasing c2 (param[4]).\n");
break;
}
x_p.CopyValues(x);
la::AddExpert(lambda, p, &x_p); // x_p = x + lambda*p
val_xp = CalculateValue(x_p); // f(x_p)
// first Wolfe condition
if (val_xp - val_x <= c1*lambda*dot_grad_p) { // f(x_p) - f(x) <= c1 * lambda * <grad, p>
// second Wolfe condition
CalculateGradient(x_p, grad_p); // grad f(x_p)
double dot_grad_x_p = la::Dot(grad_p, p);
if (dot_grad_x_p >= c2*dot_grad_p) break; // <p,grad f(x_p)> >= c2 * <p, grad>
}
}
return lambda;
}
END_OPTIM_NAMESPACE;
#endif // GRADIENTDESCENT_H
fastlib/trunk/contrib/tqlong/optimization/NelderMead.h
+38 -32
View File
@@ -24,20 +24,26 @@ struct OptimizationParameter {
BEGIN_OPTIM_NAMESPACE;
/**********************************************************************
Implement 0-order Nelder-Mead optimization method
Variable: implement la::AddExpert(alpha, const Variable& X, Variable* Y) and Copy() for this type
Function: implement CalculateValue(const Variable&)
Implement 0-order Nelder-Mead simplex optimization method
Function: implement CalculateValue(const variable_type&)
Function::variable_type:
implement la::AddExpert(double, const variable_type&, variable_type*)
Copy(const variable_type&), CopyValues(const variable_type&)
for this type
Parameter:
alpha : Reflection (1.0)
gamma : Expansion (2.0)
rho : Contraction (0.5)
sigma : Reduction (0.5)
General params
maxIter, rTol, aTol
Specific for Nelder-Mead
alpha : Reflection (1.0)
gamma : Expansion (2.0)
rho : Contraction (0.5)
sigma : Reduction (0.5)
**********************************************************************/
template<typename Variable, typename Function> class NelderMead
template<typename Function> class NelderMead
{
public:
typedef Function function_type;
typedef Variable variable_type;
typedef typename Function::variable_type variable_type;
public:
NelderMead(function_type& f_, OptimizationParameter param_ = OptimizationParameter(100, 0.01, 0.01));
void setParam(double alpha_, double gamma_, double rho_, double sigma_) {
@@ -64,33 +70,33 @@ protected:
void updateCenter(const variable_type& x1, const variable_type& x2);
};
template<typename V, typename F>
NelderMead<V, F>::NelderMead(function_type& f_, OptimizationParameter param_)
template<typename F>
NelderMead<F>::NelderMead(function_type& f_, OptimizationParameter param_)
: f(f_), param(param_)
{
memory.Init();
val.Init();
f.Init(&center);
alpha = 1.5;
gamma = 2.0;
rho = 0.5;
sigma = 0.5;
}
template<typename V, typename F>
void NelderMead<V, F>::add(const variable_type &x)
template<typename F>
void NelderMead<F>::add(const variable_type &x)
{
double v = f.CalculateValue(x);
add(x, v);
}
template<typename V, typename F>
void NelderMead<V, F>::addSeed(const ArrayList<variable_type>& vX) {
for (int i = 0; i < vX.size(); i++)
add(vX[i]);
template<typename F>
void NelderMead<F>::addSeed(const ArrayList<variable_type>& vX) {
for (int i = 0; i < vX.size(); i++) add(vX[i]);
}
template<typename V, typename F>
void NelderMead<V, F>::add(const variable_type &x, double v)
template<typename F>
void NelderMead<F>::add(const variable_type &x, double v)
{
double pos = findPos(v);
if (pos < memory.size()) {
@@ -104,19 +110,19 @@ void NelderMead<V, F>::add(const variable_type &x, double v)
updateCenter(x);
}
template<typename V, typename F>
int NelderMead<V, F>::findPos(double v)
template<typename F>
int NelderMead<F>::findPos(double v)
{
int i = 0;
while (i < val.size() && val[i] < v) i++;
return i;
}
template<typename V, typename F>
void NelderMead<V, F>::updateCenter(const variable_type &x)
template<typename F>
void NelderMead<F>::updateCenter(const variable_type &x)
{
if (memory.size() <= 1)
center.Copy(x);
center.CopyValues(x);
else {
variable_type tmp;
tmp.Copy(center);
@@ -126,8 +132,8 @@ void NelderMead<V, F>::updateCenter(const variable_type &x)
}
}
template<typename V, typename F>
void NelderMead<V, F>::updateCenter(const variable_type &x1, const variable_type &x2)
template<typename F>
void NelderMead<F>::updateCenter(const variable_type &x1, const variable_type &x2)
{
int n = memory.size();
DEBUG_ASSERT(n > 0);
@@ -135,8 +141,8 @@ void NelderMead<V, F>::updateCenter(const variable_type &x1, const variable_type
la::AddExpert(1.0/n, x2, &center);
}
template<typename V, typename F>
double NelderMead<V, F>::optimize(variable_type &sol)
template<typename F>
double NelderMead<F>::optimize(variable_type &sol)
{
int n = memory.size();
DEBUG_ASSERT(n >= 3);
@@ -195,8 +201,8 @@ double NelderMead<V, F>::optimize(variable_type &sol)
return best_val;
}
template<typename V, typename F>
void NelderMead<V, F>::replace(int pos, const variable_type &x, double v)
template<typename F>
void NelderMead<F>::replace(int pos, const variable_type &x, double v)
{
int n = memory.size();
DEBUG_ASSERT(pos < n);
@@ -218,8 +224,8 @@ void NelderMead<V, F>::replace(int pos, const variable_type &x, double v)
val[pos] = v;
}
template<typename V, typename F>
void NelderMead<V, F>::updateAll()
template<typename F>
void NelderMead<F>::updateAll()
{
int n = memory.size();
DEBUG_ASSERT(n > 0);
fastlib/trunk/contrib/tqlong/optimization/test1.cpp
+36 -6
View File
@@ -3,6 +3,7 @@
#include "FunctionTemplate.h"
#include "NelderMead.h"
#include "GradientDescent.h"
using namespace std;
using namespace optim;
@@ -27,20 +28,22 @@ void testNelderMead(fx_module* module) {
cout << "Nelder-Mead test ..." << endl;
int n_seed = fx_param_int(module, "seed", 10);
LengthEuclidianSquare f(2);
NelderMead<optim::LengthEuclidianSquare> algo(f);
// Seeding
ArrayList<Vector> vX;
vX.Init();
for (int i = 0; i < n_seed; i++) {
Vector x; x.Init(2);
Vector x;
f.Init(&x);
x[0] = 10+10*(double)rand()/RAND_MAX;
x[1] = 20+10*(double)rand()/RAND_MAX;
vX.PushBackCopy(x);
}
LengthEuclidianSquare f;
NelderMead<Vector, optim::LengthEuclidianSquare> algo(f);
algo.addSeed(vX);
// Optimization
Vector sol;
sol.Init(2);
double v = algo.optimize(sol);
@@ -51,11 +54,38 @@ void testNelderMead(fx_module* module) {
cout << "Nelder-Mead test succeeded." << endl;
}
void testGradientDescent(fx_module* module) {
cout << "GradientDescent test ..." << endl;
double param[] = {100, 0.00001, 0.001, 1e-4, 0.9, 0.4};
LengthEuclidianSquare f(2);
GradientDescent<optim::LengthEuclidianSquare> algo(f, param);
// Seeding
Vector x0;
f.Init(&x0);
x0[0] = 10+10*(double)rand()/RAND_MAX;
x0[1] = 20+10*(double)rand()/RAND_MAX;
// Optimization
Vector sol;
f.Init(&sol);
algo.setX0(x0);
double v = algo.optimize(sol);
cout << "Best value = " << v << endl;
ot::Print(sol, "Solution", stdout);
//algo.printHistory();
cout << "GradientDescent test succeeded." << endl;
}
int main(int argc, char** argv) {
fx_module* root = fx_init(argc, argv, &optimization_doc);
cout << "Optimization tests" << endl;
testNelderMead(root);
//testNelderMead(root);
testGradientDescent(root);
fx_done(root);
return 0;