专栏名称: 码农经济学
Something Anything Nothi…
目录
51好读  ›  专栏  ›  码农经济学

使用NLopt做简单优化

码农经济学  · 知乎专栏  ·  · 2016-12-13 21:05

正文

工作中经常遇到优化的问题,比如多数的统计方法最终都可以归结为一个优化问题。一般的统计和数值计算软件都会包含最优化算法的函数,比如Matlab中的fminsearch等,不过对于C等其他语言,很多时候可以选择的并不多。

之前工作中经常用到的是GSL中的优化算法,不过由于最近更经常使用MKL,所以希望能摆脱对GSL的依赖。这个时候找到了NLopt这个包。

NLopt使用起来非常简单,且支持非常多的语言,常见的语言比如C/ C++/ Julia/ Python/ R/ Fortran/ Lua/ OCaml/ Octave等都支持,所以算是一个“一招鲜,吃遍天”的包。除此之外,NLopt还有很多其他优点,比如:

  1. 可以方便的在不更改代码的情况下更改算法进行尝试,
  2. 支持一些常见的全局最优算法
  3. 支持一些高维的问题
  4. 支持常见的需要导数和不需要导数的多种算法
  5. 支持等式、不等式、有界等不同的约束极值问题

总之,很好用!

NLopt的安装也非常简单,与一般的源码安装过程一样,从官网 NLopt - AbInitio 下载源码,然后:

./configure
make
sudo make install

就可以了,默认安装在/usr/local/lib/下。windows下可以用MinGW。

NLopt支持的算法可以从 NLopt Algorithms 查询,包括:

下面是如下问题的一个实例代码:

\max_{x_1,x_2}\ln x_1+\ln x_2
s.t.
 p_1 \cdot x_1+p_2\cdot x_2=5
x_1\leq x_2, x_1\geq 0, x_2\geq 0

注意其中有一个等式约束和一个不等式约束。我们使用Sequential Least-Squares Quadratic Programming (SLSQP) 算法来求解这个优化问题。

代码:

#include <stdio.h>
#include <math.h>
#include "nlopt.h"
#define INF (1.0/0.0)

double utility(unsigned n, const double *x, double *grad, void *data){
  grad[0]=1.0/x[0];
  grad[1]=1.0/x[1];
  printf("%f, %f, %f ", x[0],x[1],log(x[0])+log(x[1]));
  return log(x[0])+log(x[1]);
}

double constraint(unsigned n, const double *x, double *grad, void *data){
  double *p=(double *)data;
  grad[0]=*p;
  grad[1]=*(p+1);
  printf("Constraint: %f\n", x[0]*(*p)+x[1]*(*(p+1))-5);
  return x[0]*(*p)+x[1]*(*(p+1))-5;
}

double inconstraint(unsigned n, const double *x, double *grad, void *data){
  grad[0]=1;
  grad[1]=-1;
  return x[0]-x[1];
}

int main(int argc, char const *argv[]) {
  double p[2]={1,2};
  double tol=1e-8;
  double lb[2]={0,0};
  double ub[2]={INF,INF};
  double x[2]={1,1};
  double f_max=-INF;
  // set up optimizer
  nlopt_opt opter=nlopt_create(NLOPT_LD_SLSQP, 2);
  // lower and upper bound
  nlopt_set_lower_bounds(opter, lb);
  nlopt_set_upper_bounds(opter, ub);
  // objective function
  nlopt_set_max_objective(opter, utility, NULL);
  // equality constraint
  nlopt_add_equality_constraint(opter, constraint, p, tol);
  // inequality constraint
  nlopt_add_inequality_constraint(opter, inconstraint, NULL, tol);
  // stopping criterion
  nlopt_set_xtol_rel(opter, tol);
  nlopt_set_ftol_abs(opter, tol);
  nlopt_set_force_stop(opter, tol);
  // optimize
  nlopt_result result=nlopt_optimize(opter, x, &f_max);
  if (result)
    printf("Maximum utility=%f, x=(%f,%f)\n", f_max, x[0], x[1]);
  // free
  nlopt_destroy(opter);
  return 0;
}






请到「今天看啥」查看全文