这是Matlab代码的一个简单翻译,原作者有关信息参见:
http://www.kyb.tuebingen.mpg.de/bs/people/carl/code/minimize/
mlminimize.h
mlminimize.cpp
使用样例
http://www.kyb.tuebingen.mpg.de/bs/people/carl/code/minimize/
mlminimize.h
- #ifndef GUARD_mlminimize_h
- #define GUARD_mlminimize_h
- #include <ml.h>
- struct CV_EXPORTS CvMinimizeParams
- {
- double INT;
- double EXT;
- int ITMAX;
- double RATIO;
- double SIG;
- double RHO;
- CvMinimizeParams()
- : INT(.1),
- EXT(3.),
- ITMAX(20),
- RATIO(10.),
- SIG(.1),
- RHO(.05)
- {}
- CvMinimizeParams( double _INT, double _EXT, int _ITMAX, double _RATIO, double _SIG, double _RHO )
- : INT(_INT),
- EXT(_EXT),
- ITMAX(_ITMAX),
- RATIO(_RATIO),
- SIG(_SIG),
- RHO(_RHO)
- {}
- };
- class CV_EXPORTS CvMinimize
- {
- private:
- CvMinimizeParams* params;
- virtual bool function( CvMat* x, double& result ) = 0;
- virtual bool derivative( CvMat* x, CvMat* result ) = 0;
- public:
- CvMinimize( CvMinimizeParams* _params )
- : params(_params)
- {}
- CvMat* minimize( CvMat* x, int length, double red = 1. );
- };
- #endif
- #include "mlminimize.h"
- CvMat* CvMinimize::minimize( CvMat* x,
- int length,
- double red )
- {
- CvMat *df0, *df3, *dF0;
- CvMat *s;
- CvMat *X, *X0, *Xn;
- X = cvCreateMat( x->rows, x->cols, CV_64FC1 );
- cvCopy( x, X );
- dF0 = cvCreateMat( X->rows, X->cols, CV_64FC1 );
- df0 = cvCreateMat( X->rows, X->cols, CV_64FC1 );
- df3 = cvCreateMat( X->rows, X->cols, CV_64FC1 );
- s = cvCreateMat( X->rows, X->cols, CV_64FC1 );
- X0 = cvCreateMat( X->rows, X->cols, CV_64FC1 );
- Xn = cvCreateMat( X->rows, X->cols, CV_64FC1 );
- cvZero( dF0 );
- cvZero( df0 );
- cvZero( df3 );
- cvZero( s );
- cvZero( X0 );
- cvZero( Xn );
- double F0 = 0, f0 = 0, f1 = 0, f2 = 0, f3 = 0, f4 = 0;
- double x1 = 0, x2 = 0, x3 = 0, x4 = 0;
- double d0 = 0, d1 = 0, d2 = 0, d3 = 0, d4 = 0;
- double A = 0, B = 0;
- bool ls_failed = 0;
- function( X, f0 );
- derivative( X, df0 );
- cvSubRS( df0, cvScalar(0), s );
- d0 = -cvDotProduct( s, s );
- x3 = red/(1.-d0);
- int i = 0;
- int l = ( length > 0 ) ? length : -length;
- int ls = ( length > 0 ) ? 1 : 0;
- int eh = ( length > 0 ) ? 0 : 1;
- while ( i < l )
- {
- i+=ls;
- cvCopy( X, X0 );
- F0 = f0;
- cvCopy( df0, dF0 );
- int m = ( length > 0 ) ? params->ITMAX : l-i;
- if ( params->ITMAX < m )
- m = params->ITMAX;
- for ( ; ; )
- {
- x2 = 0;
- f2 = f0;
- d2 = d0;
- f3 = f0;
- cvCopy( df0, df3 );
- while ( m > 0 )
- {
- m--;
- i+=eh;
- cvScaleAdd( s, cvScalar(x3), X, Xn );
- cvZero( df3 );
- if ((function( Xn, f3 ))&&(derivative( Xn, df3 )))
- break;
- else
- x3 = (x2+x3)*.5;
- }
- if ( f3 < F0 )
- {
- cvCopy( Xn, X0 );
- F0 = f3;
- cvCopy( df3, dF0 );
- }
- if ( (d3 > params->SIG*d0)||(f3 > f0+x3*params->RHO*d0)||(m <= 0) )
- break;
- x1 = x2;
- f1 = f2;
- d1 = d2;
- x2 = x3;
- f2 = f3;
- d2 = d3;
- A = 6.*(f1-f2)+3.*(d2+d1)*(x2-x1);
- B = 3.*(f2-f1)-(2.*d1+d2)*(x2-x1);
- x3 = B*B-A*d1*(x2-x1);
- if ( x3 < 0 )
- x3 = x2*params->EXT;
- else {
- x3 = x1-d1*(x2-x1)*(x2-x1)/(B+sqrt(x3));
- if ( x3 < 0 )
- x3 = x2*params->EXT;
- else {
- if ( x3 > x2*params->EXT )
- x3 = x2*params->EXT;
- else if ( x3 < x2+params->INT*(x2-x1) )
- x3 = x2+params->INT*(x2-x1);
- }
- }
- }
- while ( ((fabs(d3) > -params->SIG*d0)||(f3 > f0+x3*params->RHO*d0 ))&&(m > 0) )
- {
- if ( (d3 > 1e-8)||(f3 > f0+x3*params->RHO*d0) )
- {
- x4 = x3;
- f4 = f3;
- d4 = d3;
- } else {
- x2 = x3;
- f2 = f3;
- d2 = d3;
- }
- if ( f4 > f0 )
- x3 = x2-(.5*d2*(x4-x2)*(x4-x2))/(f4-f2-d2*(x4-x2));
- else {
- A = 6.*(f2-f4)/(x4-x2)+3.*(d4+d2);
- B = 3.*(f4-f2)-(2.*d2+d4)*(x4-x2);
- x3 = B*B-A*d2*(x4-x2)*(x4-x2);
- if ( x3 < 0 )
- x3 = (x2+x4)*.5;
- else
- x3 = x2+(sqrt(x3)-B)/A;
- }
- if ( x3 > x4-params->INT*(x4-x2) )
- x3 = x4-params->INT*(x4-x2);
- if ( x3 < x2+params->INT*(x4-x2) )
- x3 = x2+params->INT*(x4-x2);
- cvScaleAdd( s, cvScalar(x3), X, Xn );
- function( Xn, f3 );
- cvZero( df3 );
- derivative( Xn, df3 );
- if ( f3 < F0 )
- {
- cvCopy( Xn, X0 );
- F0 = f3;
- cvCopy( df3, dF0 );
- }
- m--;
- i+=eh;
- d3 = cvDotProduct( df3, s );
- }
- if ( (fabs(d3) < -params->SIG*d0)&&(f3 < f0+x3*params->RHO*d0) )
- {
- cvCopy( Xn, X );
- f0 = f3;
- cvScaleAdd( s, cvScalar((cvDotProduct( df0, df3 )-cvDotProduct( df3, df3 ))/cvDotProduct( df0, df0 )), df3, s );
- cvSubRS( s, cvScalar(0), s );
- cvCopy( df3, df0 );
- d3 = d0;
- d0 = cvDotProduct( df0, s );
- if ( d0 > 0 )
- {
- cvSubRS( df0, cvScalar(0), s );
- d0 = -cvDotProduct( s, s );
- }
- x3 = x3*(params->RATIO < d3/(d0-1e-8) ? params->RATIO:d3/(d0-1e-8));
- ls_failed = 0;
- } else {
- cvCopy( X0, X );
- f0 = F0;
- cvCopy( dF0, df0 );
- if ( ls_failed )
- break;
- cvSubRS( df0, cvScalar(0), s );
- d0 = -cvDotProduct( s, s );
- x3 = red/(1.-d0);
- ls_failed = 1;
- }
- }
- cvReleaseMat( &s );
- cvReleaseMat( &X0 );
- cvReleaseMat( &Xn );
- cvReleaseMat( &dF0 );
- cvReleaseMat( &df0 );
- cvReleaseMat( &df3 );
- return X;
- }
- #include "mlminimize.h"
- #include <iostream>
- class CV_EXPORTS Rosenbrock : public CvMinimize
- {
- private:
- virtual bool function( CvMat* x, double& result )
- {
- double* x_vec = x->data.db;
- result = 0;
- for ( int i = 0; i < x->cols-1; i++ )
- result += 100*(x_vec[i+1]-x_vec[i]*x_vec[i])*(x_vec[i+1]-x_vec[i]*x_vec[i])+(1-x_vec[i])*(1-x_vec[i]);
- return 1;
- }
- virtual bool derivative( CvMat* x, CvMat* result )
- {
- double* x_vec = x->data.db;
- double* result_vec = result->data.db;
- for ( int i = 0; i < x->cols-1; i++ )
- result_vec[i] = -400*x_vec[i]*(x_vec[i+1]-x_vec[i]*x_vec[i])-2*(1-x_vec[i]);
- for ( int i = 1; i < x->cols; i++ )
- result_vec[i] += 200*(x_vec[i]-x_vec[i-1]*x_vec[i-1]);
- return 1;
- }
- public:
- Rosenbrock( CvMinimizeParams* params )
- : CvMinimize(params) {}
- bool derivative1( CvMat* x, CvMat* result )
- {
- double* x_vec = x->data.db;
- double* result_vec = result->data.db;
- for ( int i = 0; i < x->cols-1; i++ )
- result_vec[i] = -400*x_vec[i]*(x_vec[i+1]-x_vec[i]*x_vec[i])-2*(1-x_vec[i]);
- for ( int i = 1; i < x->cols; i++ )
- result_vec[i] += 200*(x_vec[i]-x_vec[i-1]*x_vec[i-1]);
- return 1;
- }
- };
- int main()
- {
- CvMinimizeParams* params = new CvMinimizeParams();
- Rosenbrock* rb = new Rosenbrock( params );
- CvMat* X = cvCreateMat( 1, 2, CV_64FC1 );
- cvZero( X );
- CvMat* result = rb->minimize( X, 25 );
- for ( int k = 0; k < 2; k++ )
- printf("%f/n", result->data.db[k]);
- }