• GSL 求多元函数极值:simplex 算法


    GSL 手册上有关于 gsl_multimin_fminimizer_nmsimplex2 的说明。它只需要被优化函数 \(f(\vec{x})\),而不需要偏导数。它自己有个 Nelder-Mead 算法(我不懂这个,还没有研究过,只是试用过),会自己去取点判断下降方向。

    手册上有测试代码:

    
    #include<iostream>
    using namespace std;
    
    #include<gsl/gsl_multimin.h>
    #include<gsl/gsl_sf_bessel.h>
    
    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 = "<<count_myf<<endl;
    
            gsl_vector_free(x);
            gsl_vector_free(ss);
            gsl_multimin_fminimizer_free(s);
    
            return status;
    }
    

    编译+运行:

    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 次迭代,就会收敛:https://www.cnblogs.com/luyi07/p/14560924.html
    所以共轭梯度法收敛更快,虽然它需要偏导数。

    但如果偏导数很昂贵,不妨拿 nmsimplex2 试一试。

  • 相关阅读:
    组合模式
    MySQL8.0 下载安装启动(Windows10)
    OI如逆旅,我亦是行人——省选
    闲话—江湖痴情浅,信步余生。平剑红烛,青丝微绾,却话奁中。
    此时彼方
    CSP 2019游记 & 退役记
    西狂 杨过
    SDOI 2019 Round1 游记
    NOIP2018游记
    未来可期,不知所终
  • 原文地址:https://www.cnblogs.com/luyi07/p/16280519.html
Copyright © 2020-2023  润新知