GSL 求多元函数极值:simplex 算法
GSL 手册上有关于 gsl_multimin_fminimizer_nmsimplex2 的说明。它只需要被优化函数 \(f(\vec{x})\),而不需要偏导数。它自己有个 Nelder-Mead 算法(我不懂这个,还没有研究过,只是试用过),会自己去取点判断下降方向。
手册上有测试代码:
#include
using namespace std;
#include
#include
int count_myf = 0;
/* Paraboloid centered on (p[0],p[1]), with
scale factors (p[2],p[3]) and minimum p[4] */
double my_f (const gsl_vector *v, void *params)
{
count_myf ++;
double x, y;
double *p = (double *)params;
x = gsl_vector_get(v, 0);
y = gsl_vector_get(v, 1);
return p[2] * (x - p[0]) * (x - p[0]) +
p[3] * (y - p[1]) * (y - p[1]) + p[4];
}
int main(){
// define the multi-var function to be optimized
gsl_multimin_function minex_func;
/* Paraboloid center at (1,2), scale factors (10, 20), minimum value 30 */
double par[5] = { 1.0, 2.0, 10.0, 20.0, 30.0 };
minex_func.n = 2; /* number of function components */
minex_func.f = my_f;
minex_func.params = par;
// set initial position x
gsl_vector *x;
x = gsl_vector_alloc(2);
gsl_vector_set( x, 0, 5.0 );
gsl_vector_set( x, 1, 7.0 );
// set initial step size to 1
gsl_vector *ss = gsl_vector_alloc(2);
gsl_vector_set_all( ss, 1.0 );
// set minimizer
const gsl_multimin_fminimizer_type *T = gsl_multimin_fminimizer_nmsimplex2;
gsl_multimin_fminimizer *s = gsl_multimin_fminimizer_alloc( T, 2 );
gsl_multimin_fminimizer_set( s, & minex_func, x, ss );
int iter = 0, status; double size;
do{
iter++;
status = gsl_multimin_fminimizer_iterate( s );
if( status ) break;
size = gsl_multimin_fminimizer_size (s);
status = gsl_multimin_test_size( size, 1E-8);
if( status == GSL_SUCCESS )
printf( " Converged to minimum at :\n" );
printf( "%5d %.10f %.10f f() = %7.3f size = %.10f\n",
iter,
gsl_vector_get( s->x, 0 ),
gsl_vector_get( s->x, 1 ),
s->fval, size );
}
while ( status == GSL_CONTINUE && iter < 100 );
cout<<" count_myf = "<
编译+运行:
g++ test_simplex.cpp -lgsl -lgslcblas
./a.out
得到结果
1 6.5000000000 5.0000000000 f() = 512.500 size = 1.1303883305
2 5.2500000000 4.0000000000 f() = 290.625 size = 1.4092945438
3 5.2500000000 4.0000000000 f() = 290.625 size = 1.4092945438
4 5.5000000000 1.0000000000 f() = 252.500 size = 1.4092945438
5 2.6250000000 3.5000000000 f() = 101.406 size = 1.8474832731
6 2.6250000000 3.5000000000 f() = 101.406 size = 1.8474832731
7 -0.0000000000 3.0000000000 f() = 60.000 size = 1.8474832731
8 2.0937500000 1.8750000000 f() = 42.275 size = 1.3213162892
9 0.2578125000 1.9062500000 f() = 35.684 size = 1.0689461829
10 0.5878906250 2.4453125000 f() = 35.664 size = 0.8409024133
11 1.2583007812 2.0253906250 f() = 30.680 size = 0.4761530328
12 1.2583007812 2.0253906250 f() = 30.680 size = 0.3672920466
13 1.0926208496 1.8494873047 f() = 30.539 size = 0.2995597113
14 0.8829574585 2.0041198730 f() = 30.137 size = 0.1724325444
15 0.8829574585 2.0041198730 f() = 30.137 size = 0.1261625606
16 0.9581913948 2.0604190826 f() = 30.090 size = 0.1062198609
17 1.0218096972 2.0041832924 f() = 30.005 size = 0.0626448959
18 1.0218096972 2.0041832924 f() = 30.005 size = 0.0433856521
19 1.0218096972 2.0041832924 f() = 30.005 size = 0.0433856521
20 1.0218096972 2.0041832924 f() = 30.005 size = 0.0274263092
21 1.0218096972 2.0041832924 f() = 30.005 size = 0.0218800866
22 0.9920430849 1.9969168063 f() = 30.001 size = 0.0156702845
23 0.9920430849 1.9969168063 f() = 30.001 size = 0.0133695954
24 0.9920430849 1.9969168063 f() = 30.001 size = 0.0079586570
25 0.9987625058 2.0047084907 f() = 30.000 size = 0.0079586570
26 1.0025252407 1.9999882754 f() = 30.000 size = 0.0053914447
27 1.0025252407 1.9999882754 f() = 30.000 size = 0.0034382764
28 1.0025252407 1.9999882754 f() = 30.000 size = 0.0027835269
29 0.9985776580 2.0003782319 f() = 30.000 size = 0.0020123788
30 0.9985776580 2.0003782319 f() = 30.000 size = 0.0017260798
31 1.0008632701 2.0003940352 f() = 30.000 size = 0.0010139823
32 0.9996159870 1.9995509089 f() = 30.000 size = 0.0010139823
33 0.9994086433 2.0001753520 f() = 30.000 size = 0.0007350933
34 1.0001877926 2.0001285828 f() = 30.000 size = 0.0004349777
35 1.0001877926 2.0001285828 f() = 30.000 size = 0.0003513674
36 1.0002168497 1.9998973397 f() = 30.000 size = 0.0002633416
37 0.9999547118 1.9999321997 f() = 30.000 size = 0.0001553283
38 0.9999547118 1.9999321997 f() = 30.000 size = 0.0001215448
39 0.9999601990 2.0000167371 f() = 30.000 size = 0.0000940103
40 0.9999601990 2.0000167371 f() = 30.000 size = 0.0000557365
41 0.9999601990 2.0000167371 f() = 30.000 size = 0.0000420072
42 0.9999914230 1.9999886035 f() = 30.000 size = 0.0000378037
43 1.0000114660 2.0000003713 f() = 30.000 size = 0.0000240432
44 1.0000114660 2.0000003713 f() = 30.000 size = 0.0000145614
45 0.9999911331 2.0000000498 f() = 30.000 size = 0.0000109784
46 0.9999963613 1.9999944070 f() = 30.000 size = 0.0000090450
47 1.0000026066 1.9999987999 f() = 30.000 size = 0.0000052764
48 1.0000026066 1.9999987999 f() = 30.000 size = 0.0000037733
49 1.0000026066 1.9999987999 f() = 30.000 size = 0.0000037733
50 0.9999986944 1.9999995432 f() = 30.000 size = 0.0000023683
51 0.9999986944 1.9999995432 f() = 30.000 size = 0.0000018371
52 1.0000012525 1.9999995221 f() = 30.000 size = 0.0000013433
53 1.0000005378 2.0000002391 f() = 30.000 size = 0.0000011223
54 0.9999997948 1.9999997119 f() = 30.000 size = 0.0000006582
55 0.9999997948 1.9999997119 f() = 30.000 size = 0.0000004499
56 1.0000001234 2.0000000980 f() = 30.000 size = 0.0000002732
57 1.0000001234 2.0000000980 f() = 30.000 size = 0.0000002028
58 1.0000001234 2.0000000980 f() = 30.000 size = 0.0000001210
59 0.9999998954 2.0000000247 f() = 30.000 size = 0.0000000827
60 0.9999999427 1.9999999776 f() = 30.000 size = 0.0000001100
61 0.9999999427 1.9999999776 f() = 30.000 size = 0.0000000599
62 0.9999999427 1.9999999776 f() = 30.000 size = 0.0000000599
63 1.0000000134 2.0000000198 f() = 30.000 size = 0.0000000543
64 1.0000000233 2.0000000006 f() = 30.000 size = 0.0000000398
65 0.9999999805 1.9999999939 f() = 30.000 size = 0.0000000213
66 0.9999999805 1.9999999939 f() = 30.000 size = 0.0000000187
67 1.0000000087 2.0000000009 f() = 30.000 size = 0.0000000143
Converged to minimum at :
68 0.9999999944 1.9999999993 f() = 30.000 size = 0.0000000077
count_myf = 134
可见,nmsimplex2 需要 68 次迭代,才真正收敛。
之前的帖子里记录了,共轭梯度法只需要 13 次迭代,就会收敛:
所以共轭梯度法收敛更快,虽然它需要偏导数。
但如果偏导数很昂贵,不妨拿 nmsimplex2 试一试。