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 试一试。

相关