#include<cmath>

#include<ctime>

using namespace std

double f(double x)

double g(double x)

double gd(double xs,double s)

{

    double x=xs

    double y

    for(int i=0i!=20++i)

    {

        double grad= -1*g(x)

        x+=grad*s

        y=f(x)

        cout <<"EPOCH "<<i<<" grad = " <<setprecision(15)<<  grad<< " x = "<< x <<"  y = " << y<< endl

        if(abs(grad) < 1e-6)

            break

    }

    return x

}

int main()

{

    double xk=-5,ak=0.1

    gd(xk,ak)

}

double f(double x)

{

    return x*x-2*x+1

}

double g(double x)

{

return 2*x-2

}

minf = 0

for x1 = -2 : 0.1 : 1

for x2 = -2 : 0.1 : 1

f = 4*x1+2*x2+10*x1*x2+6*x1^2+5*x2^2

if f <minf

minf = f

mx1 = x1

mx2 = x2

end

end

end

[minf mx1 mx2]

ans =

-1.2000 -1.00000.8000

当x1 = -1,x2 = 0.8时，取得最小值-1.2

一书及所附程序，有完整程序。不过我封装了它的C++版本，可以对但参数或多参数求极值，完整的头文件为：

#ifndef __OPTIMIZATION_H__

#define __OPTIMIZATION_H__

//////////////////////////////////////////////////////////////////////////

// class TOptimization

//

// $ 求函数一个或多个参数的最小值

//

// 该类默认对一个参数优化，一般只要输入优化参数数目

// 优化目标函数就可以使用。

//

// ...主要代码来自：

// Numerical Recipes in C++

// The Art of Scientific Computing

// Second Edition

// William H. Press Saul A. Teukolsky

// William T. Vetterling Brian P. Flannery

//

// 中译本：

// C++ 数值算法(第二版) 胡健伟 赵志勇 薛运华 等译

// 电子工业出版社 北京 (2005)

//

// Author: Jian Feng

// Email: [email protected]

// Dec. 9, 2006

//

//////////////////////////////////////////////////////////////////////////

//

// 输入函数：

//

// @MaxIterationStep: 最大迭代次数，默认 1000

// @ParameterNumbers: 优化参数数目，默认 1

// @InitMatrix:初始化矩阵参数(N*N)， 默认

// @Tolerance: 容许公差，默认 1E-7

//

// 执行函数：

//

// @ExecutePowell: 利用powell方法进行多参数优化

// @ExecuteBrent: 利用brent方法进行单参数优化

//

// 输出函数：

//

// @OptimizatedParameters: 优化结果数据

// @ObjectiveFunctionValue: 目标函数在优化值处的值

//

// 使用示例：

//

// 1. 单参数

// double objfun(double a){

// double sum = 0

// for(int i = 0i <DataPoints++i)

// sum += SQR(Exps[i] - Theo(a))

// }

// double value

// TOptimization opt

// if(opt.ExecuteBrent(objfun, -10, -1)) opt.OptimizatedParameters(&value)

//

// 2. 多参数

// double objfun(double *a){

// double sum = 0

// for(int i = 0i <DataPoints++i)

// sum += SQR(Exps[i] - Theo(a))

// }

// double value[3]

// TOptimization opt(3)

// double ival[3] = {-1, 0, 1}

// if(opt.ExecutePowell(objfun, ival)) opt.OptimizatedParameters(value)

//

namespace{

static int ncom //公用变量

static double *pcom_p//公用变量

static double *xicom_p //公用变量

static double (*nrfunc)(double*) //公用函数指针

}

class TOptimization

{

private:

typedef double (*Reff)(double *)

typedef double (*Ptrf)(double )

public:

TOptimization(int n = 1)

~TOptimization(){ FreeMemory()}

//主要方法

void ParameterNumbers(int n){ FreeMemory()num = nAllocateMemory()}

//利用powell方法对一个或多个参数优化

bool ExecutePowell(Reff obj, double *a = 0)

//利用brent方法对一个参数优化，需给出参数所在的区间

bool ExecuteBrent(Ptrf obj, double vFrom = 0, double vTo = 1)

void OptimizatedParameters(double *a){ for(int i=0i<num++i) a[i]=coef[i]}

void OptimizatedParameters(double &a){ a = vmin }

//void OptimizatedParameters(double *a){

// if(method) for(int i=0i<num++i) a[i]=coef[i]

// else *a = vmin

//}

//其它方法

void InitMatrix(double **m)

{

for(int i=0i<num++i)

for(int j = 0j<num++j)

matx[i][j]=m[i][j]

setm = true

}

void MaxIterationStep(int s){ ITMAX = s}

void Tolerance(double eps){ ftol = eps}

double ObjectiveFunctionValue()const{ return fret }

private:

double brent(double ax, double bx, double cx, Ptrf f, double tol, double &xmin, int &flag)

void mnbrak(double &ax, double &bx, double &cx, double &fa, double &fb, double &fc, Ptrf func)

void linmin(double *p, double *xi, double &fret, Reff func)

bool powell(double *p, double **xi, double ftol, int &iter, double &fret, Reff func)

void shft2(double &a, double &b, const double c){ a=bb=c}

void shft3(double &a, double &b, double &c, const double d){ a=bb=cc=d}

double SQR(double x){ return x * x}

void SWAP(double &a, double &b){ double dum=aa=bb=dum}

double SIGN(const double &a, const double &b){return b >= 0?(a>=0?a:-a):(a>=0?-a:a)}

double MAX(const double &a, const double &b){return b >a ? (b) : (a)}

void AllocateMemory()

void FreeMemory()

static double f1dim(double x)

{

int j

double *xt = new double [ncom]

//Vec_Dp &pcom=*pcom_p,&xicom=*xicom_p

double *pcom = pcom_p, *xicom = xicom_p

for (j=0j<ncomj++)

xt[j]=pcom[j]+x*xicom[j]

//delete []xt

double val = nrfunc(xt)

delete []xt

return val

}

bool setm //是否设置优化方向初始矩阵

int num//优化参数

int ITMAX //最大迭代数

int iter //实际迭代步数

int method //优化方法 0: 1-D brent, 2: M-D Powell

double vmin//一维优化参数

double ftol//容许差

double fret//目标函数值

double *coef //多维优化参数值

double **matx //多维优化参数方向的初始值

}

//////////////////////////////////////////////////////////////////////////

inline TOptimization::TOptimization(int n )

{

num= n

ftol = 1e-7

ITMAX = 1000

iter = 0

fret = 0.

vmin = 0.

method = 0

setm = false

AllocateMemory()

}

inline void TOptimization::AllocateMemory()

{

pcom_p = new double [num]

xicom_p = new double [num]

coef= new double [num]

matx= new double *[num]

for(int i = 0i <num++i)

{

coef[i] = 0.

matx[i] = new double [num]

for(int j = 0j <num++j)

matx[i][j]=(i == j ? 1.0 : 0.0)

}

}

inline void TOptimization::FreeMemory()

{

for(int i = 0i <num++i)

{

delete []matx[i]

}

delete []matx

delete []pcom_p

delete []xicom_p

delete []coef

}

inline bool TOptimization::ExecutePowell(Reff obj, double *a)

{

method = 1

if(a)

for(int i = 0i <num++i) coef[i] = a[i]

return powell(coef, matx, ftol, iter, fret, obj)

}

inline bool TOptimization::ExecuteBrent(Ptrf obj, double vFrom, double vTo)

{

method = 0

int flag

double cx, fa, fb, fc

mnbrak(vFrom,vTo,cx,fa,fb,fc,obj)

fret = brent(vFrom,vTo,cx,obj, ftol,vmin, flag)

return flag ? true : false

}

inline void TOptimization::mnbrak(double &ax, double &bx, double &cx, double &fa,

double &fb, double &fc, Ptrf func)

{

const double GOLD=1.618034,GLIMIT=100.0,TINY=1.0e-20

double ulim,u,r,q,fu

fa=func(ax)

fb=func(bx)

if (fb >fa) {

SWAP(ax,bx)

SWAP(fb,fa)

}

cx=bx+GOLD*(bx-ax)

fc=func(cx)

while (fb >fc) {

r=(bx-ax)*(fb-fc)

q=(bx-cx)*(fb-fa)

u=bx-((bx-cx)*q-(bx-ax)*r)/

(2.0*SIGN(MAX(fabs(q-r),TINY),q-r))

ulim=bx+GLIMIT*(cx-bx)

if ((bx-u)*(u-cx) >0.0) {

fu=func(u)

if (fu <fc) {

ax=bx

bx=u

fa=fb

fb=fu

return

} else if (fu >fb) {

cx=u

fc=fu

return

}

u=cx+GOLD*(cx-bx)

fu=func(u)

} else if ((cx-u)*(u-ulim) >0.0) {

fu=func(u)

if (fu <fc) {

shft3(bx,cx,u,cx+GOLD*(cx-bx))

shft3(fb,fc,fu,func(u))

}

} else if ((u-ulim)*(ulim-cx) >= 0.0) {

u=ulim

fu=func(u)

} else {

u=cx+GOLD*(cx-bx)

fu=func(u)

}

shft3(ax,bx,cx,u)

shft3(fa,fb,fc,fu)

}

}

inline double TOptimization::brent(double ax, double bx, double cx,

Ptrf f, double tol, double &xmin, int &flag)

{

flag = 1

const double CGOLD=0.3819660

const double ZEPS=1.0e-20

int iter

double a,b,d=0.0,etemp,fu,fv,fw,fx

double p,q,r,tol1,tol2,u,v,w,x,xm

double e=0.0

a=(ax <cx ? ax : cx)

b=(ax >cx ? ax : cx)

x=w=v=bx

fw=fv=fx=f(x)

for (iter=0iter<ITMAXiter++) {

xm=0.5*(a+b)

tol2=2.0*(tol1=tol*fabs(x)+ZEPS)

if (fabs(x-xm) <= (tol2-0.5*(b-a))) {

xmin=x

return fx

}

if (fabs(e) >tol1) {

r=(x-w)*(fx-fv)

q=(x-v)*(fx-fw)

p=(x-v)*q-(x-w)*r

q=2.0*(q-r)

if (q >0.0) p = -p

q=fabs(q)

etemp=e

e=d

if (fabs(p) >= fabs(0.5*q*etemp) || p <= q*(a-x) || p >= q*(b-x))

d=CGOLD*(e=(x >= xm ? a-x : b-x))

else {

d=p/q

u=x+d

if (u-a <tol2 || b-u <tol2)

d=SIGN(tol1,xm-x)

}

} else {

d=CGOLD*(e=(x >= xm ? a-x : b-x))

}

u=(fabs(d) >= tol1 ? x+d : x+SIGN(tol1,d))

fu=f(u)

if (fu <= fx) {

if (u >= x) a=xelse b=x

shft3(v,w,x,u)

shft3(fv,fw,fx,fu)

} else {

if (u <x) a=uelse b=u

if (fu <= fw || w == x) {

v=w

w=u

fv=fw

fw=fu

} else if (fu <= fv || v == x || v == w) {

v=u

fv=fu

}

}

}

flag = 0

xmin=x

return fx

}

inline void TOptimization::linmin(double *p, double *xi, double &fret, Reff func)

{

int j, flag

const double TOL=1.0e-8

double xx,xmin,fx,fb,fa,bx,ax

int n=num

ncom=n

//pcom_p=new Vec_Dp(n)

//xicom_p=new Vec_Dp(n)

nrfunc=func

//Vec_Dp &pcom=*pcom_p,&xicom=*xicom_p

double *pcom = pcom_p, *xicom = xicom_p

for (j=0j<nj++) {

pcom[j]=p[j]

xicom[j]=xi[j]

}

ax=0.0

xx=1.0

mnbrak(ax,xx,bx,fa,fx,fb,f1dim)

fret=brent(ax,xx,bx,f1dim,TOL,xmin, flag)

for (j=0j<nj++) {

xi[j] *= xmin

p[j] += xi[j]

}

//delete xicom_p

//delete pcom_p

}

inline bool TOptimization::powell(double *p, double **xi, double ftol, int &iter,

double &fret, Reff func)

{

const int ITMAX=500

const double TINY=1.0e-20

int i,j,ibig

double del,fp,fptt,t

int n=num

//Vec_Dp pt(n),ptt(n),xit(n)

double *pt, *ptt, *xit

for(i = 0i <n++i)

{

pt = new double [n]

ptt = new double [n]

xit = new double [n]

}

fret=func(p)

for (j=0j<nj++) pt[j]=p[j]

for (iter=0++iter) {

fp=fret

ibig=0

del=0.0

for (i=0i<ni++) {

for (j=0j<nj++) xit[j]=xi[j][i]

fptt=fret

linmin(p,xit,fret,func)

if (fptt-fret >del) {

del=fptt-fret

ibig=i+1

}

}

if (2.0*(fp-fret) <= ftol*(fabs(fp)+fabs(fret))+TINY) {

delete []pt

delete []ptt

delete []xit

return true

}

if (iter == ITMAX)

{

delete []pt

delete []ptt

delete []xit

return false

//cerr<<"powell exceeding maximum iterations."

}

for (j=0j<nj++) {

ptt[j]=2.0*p[j]-pt[j]

xit[j]=p[j]-pt[j]

pt[j]=p[j]

}

fptt=func(ptt)

if (fptt <fp) {

t=2.0*(fp-2.0*fret+fptt)*SQR(fp-fret-del)-del*SQR(fp-fptt)

if (t <0.0) {

linmin(p,xit,fret,func)

for (j=0j<nj++) {

xi[j][ibig-1]=xi[j][n-1]

xi[j][n-1]=xit[j]

}

}

}

}

}

#endif