函数最小化与最大化方法全解析

### 函数最小化与最大化方法全解析 #### 1. 引言 在实际应用中，我们常常会遇到这样的问题：给定一个依赖于一个或多个自变量的函数 \(f\)，需要找出这些自变量取何值时，函数 \(f\) 达到最大值或最小值，进而计算出该最大值或最小值。其实，最大化和最小化问题紧密相关，因为一个人的函数 \(f\) 可能是另一个人的 \(-f\)。在计算时，我们通常希望快速、低成本且占用较少内存。 函数的极值可以分为全局极值（函数在整个定义域内的最高或最低值）和局部极值（在有限邻域内的最高或最低值，且不在该邻域的边界上）。寻找全局极值通常是一个极具挑战性的问题，不过有两种常用的启发式方法：一是从不同的自变量初始值出发寻找局部极值，然后挑选出其中最极端的值；二是对局部极值进行有限幅度的扰动，观察是否能得到更好的点。 近年来，模拟退火方法在解决各种全局极值问题上取得了显著成效。此外，经济学家和一些工程师特别关注约束优化问题，其中线性规划是约束优化中发展较为成熟的一个领域，它涉及的函数和约束条件恰好都是自变量的线性函数。 #### 2. 一维黄金分割搜索 在一维情况下，我们可以回顾二分法求函数根的过程。类似地，在函数最小化问题中，我们需要明确“包围”最小值的含义。对于函数的根，当两个点处函数值异号时，我们知道根被这两个点包围；而对于最小值，需要三个点 \(a < b < c\)（或 \(c < b < a\)），使得 \(f(b)\) 小于 \(f(a)\) 和 \(f(c)\)，此时我们知道函数在区间 \((a, c)\) 内有最小值。 接下来，我们要选择一个新的点 \(x\)，可以在 \(a\) 和 \(b\) 之间，也可以在 \(b\) 和 \(c\) 之间。通过评估 \(f(x)\) 的值，我们可以得到一个新的包围三元组。这个过程会一直持续，直到包围区间足够小。 那么，多小才算“足够小”呢？一般来说，由于函数在最小值附近的形状可以用泰勒定理近似表示，我们很难要求包围区间的宽度小于其中心值的 \(\sqrt{\epsilon}\) 倍（\(\epsilon\) 是计算机的浮点精度）。因此，在实际应用中，我们通常将用户提供的参数 \(tol\) 设置为不小于机器浮点精度的平方根。 为了选择新的点 \(x\)，我们希望最小化最坏情况下的可能性。通过一系列的推导，我们发现最优的选择是让 \(x\) 位于较大区间中，且与 \(b\) 的距离满足黄金分割比例（约为 0.38197）。这种方法被称为黄金分割搜索，它保证每次新的函数评估都会将最小值包围在一个是前一个区间 0.61803 倍大小的区间内。 以下是初始包围最小值的代码： ```c #include <math.h> #include "nrutil.h" #define GOLD 1.618034 #define GLIMIT 100.0 #define TINY 1.0e-20 #define SHFT(a,b,c,d) (a)=(b);(b)=(c);(c)=(d); void mnbrak(float *ax, float *bx, float *cx, float *fa, float *fb, float *fc, float (*func)(float)) { float ulim,u,r,q,fu,dum; *fa=(*func)(*ax); *fb=(*func)(*bx); if (*fb > *fa) { SHFT(dum,*ax,*bx,dum) SHFT(dum,*fb,*fa,dum) } *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(FMAX(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) { SHFT(*bx,*cx,u,*cx+GOLD*(*cx-*bx)) SHFT(*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); } SHFT(*ax,*bx,*cx,u) SHFT(*fa,*fb,*fc,fu) } } ``` 以下是黄金分割搜索的代码： ```c #include <math.h> #define R 0.61803399 #define C (1.0-R) #define SHFT2(a,b,c) (a)=(b);(b)=(c); #define SHFT3(a,b,c,d) (a)=(b);(b)=(c);(c)=(d); float golden(float ax, float bx, float cx, float (*f)(float), float tol, float *xmin) { float f1,f2,x0,x1,x2,x3; x0=ax; x3=cx; if (fabs(cx-bx) > fabs(bx-ax)) { x1=bx; x2=bx+C*(cx-bx); } else { x2=bx; x1=bx-C*(bx-ax); } f1=(*f)(x1); f2=(*f)(x2); while (fabs(x3-x0) > tol*(fabs(x1)+fabs(x2))) { if (f2 < f1) { SHFT3(x0,x1,x2,R*x1+C*x3) SHFT2(f1,f2,(*f)(x2)) } else { SHFT3(x3,x2,x1,R*x2+C*x0) SHFT2(f2,f1,(*f)(x1)) } } if (f1 < f2) { *xmin=x1; return f1; } else { *xmin=x2; return f2; } } ``` #### 3. 一维抛物线插值与布伦特方法 黄金分割搜索主要是为了应对函数最小化问题中最不利的情况。然而，如果函数在最小值附近呈现出良好的抛物线形状，那么通过三个点拟合的抛物线应该能够让我们一步到位，或者至少非常接近最小值。这种通过抛物线拟合来寻找最小值横坐标的方法被称为逆抛物线插值。 逆抛物线插值的公式为： \[x = b - \frac{1}{2} \frac{(b - a)^2[f(b) - f(c)] - (b - c)^2[f(b) - f(a)]}{(b - a)[f(b) - f(c)] - (b - c)[f(b) - f(a)]}\] 不过，这个公式在三个点共线时会失效，而且它既可能跳到抛物线的最大值，也可能跳到最小值，因此单纯依赖这个公式的最小化方案在实际中很难成功。 布伦特方法巧妙地解决了这个问题。在任何特定阶段，它会跟踪六个函数点 \(a\)、\(b\)、\(u\)、\(v\)、\(w\) 和 \(x\)。其中，最小值被 \(a\) 和 \(b\) 包围，\(x\) 是到目前为止函数值最小的点，\(w\) 是函数值第二小的点，\(v\) 是 \(w\) 的前一个值，\(u\) 是最近一次评估函数的点。 布伦特方法会尝试进行抛物线插值，通过 \(x\)、\(v\) 和 \(w\) 这三个点拟合抛物线。为了使抛物线步骤可接受，它必须满足两个条件：一是落在包围区间 \((a, b)\) 内；二是从当前最佳值 \(x\) 移动的距离小于上一步移动距离的一半。在最坏的情况下，该方法会在抛物线步骤和黄金分割步骤之间交替进行，最终通过黄金分割步骤收敛。 以下是布伦特方法的代码： ```c #include <math.h> #include "nrutil.h" #define ITMAX 100 #define CGOLD 0.3819660 #define ZEPS 1.0e-10 #define SHFT(a,b,c,d) (a)=(b);(b)=(c);(c)=(d); float brent(float ax, float bx, float cx, float (*f)(float), float tol, float *xmin) { int iter; float a,b,d,etemp,fu,fv,fw,fx,p,q,r,tol1,tol2,u,v,w,x,xm; float e=0.0; a=(ax < cx ? ax : cx); b=(ax > cx ? ax : cx); x=w=v=bx; fw=fv=fx=(*f)(x); for (iter=1;iter<=ITMAX;iter++) { 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=x; else b=x; SHFT(v,w,x,u) SHFT(fv,fw,fx,fu) } else { if (u < x) a=u; else 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; } } } nrerror("Too many iterations in brent"); *xmin=x; return fx; } ``` #### 4. 带一阶导数的一维搜索 在这部分，我们的目标与前面相同，即找出被三个横坐标 \((a, b, c)\) 包围的函数最小值，但我们可以利用函数的一阶导数信息。 我们不能简单地通过寻找导数的零点来确定最小值，因为这样无法区分最大值和最小值，也无法处理初始条件下导数指示“下坡”方向在包围区间之外的情况。因此，导数的作用只能是帮助我们在包围区间内选择新的试验点。 一种保守的方法是根据包围三元组中心点处导数的符号，确定下一个试验点应在区间 \((a, b)\) 还是 \((b, c)\) 内。然后，通过割线法（逆线性插值）将该点处的导数和第二好点处的导数外推到零。同时，我们会对这个新的试验点施加与布伦特方法类似的限制，如果试验点不符合要求，我们会对当前考察的区间进行二分。 以下是带一阶导数的一维搜索代码： ```c #include <math.h> #include "nrutil.h" #define ITMAX 100 #define ZEPS 1.0e-10 #define MOV3(a,b,c, d,e,f) (a)=(d);(b)=(e);(c)=(f); float dbrent(float ax, float bx, float cx, float (*f)(float), float (*df)(float), float tol, float *xmin) { int iter,ok1,ok2; float a,b,d,d1,d2,du,dv,dw,dx,e=0.0; float fu,fv,fw,fx,olde,tol1,tol2,u,u1,u2,v,w,x,xm; a=(ax < cx ? ax : cx); b=(ax > cx ? ax : cx); x=w=v=bx; fw=fv=fx=(*f)(x); dw=dv=dx=(*df)(x); for (iter=1;iter<=ITMAX;iter++) { xm=0.5*(a+b); tol1=tol*fabs(x)+ZEPS; tol2=2.0*tol1; if (fabs(x-xm) <= (tol2-0.5*(b-a))) { *xmin=x; return fx; } if (fabs(e) > tol1) { d1=2.0*(b-a); d2=d1; if (dw != dx) d1=(w-x)*dx/(dx-dw); if (dv != dx) d2=(v-x)*dx/(dx-dv); u1=x+d1; u2=x+d2; ok1 = (a-u1)*(u1-b) > 0.0 && dx*d1 <= 0.0; ok2 = (a-u2)*(u2-b) > 0.0 && dx*d2 <= 0.0; olde=e; e=d; if (ok1 || ok2) { if (ok1 && ok2) d=(fabs(d1) < fabs(d2) ? d1 : d2); else if (ok1) d=d1; else d=d2; if (fabs(d) <= fabs(0.5*olde)) { u=x+d; if (u-a < tol2 || b-u < tol2) d=SIGN(tol1,xm-x); } else { d=0.5*(e=(dx >= 0.0 ? a-x : b-x)); } } else { d=0.5*(e=(dx >= 0.0 ? a-x : b-x)); } } else { d=0.5*(e=(dx >= 0.0 ? a-x : b-x)); } if (fabs(d) >= tol1) { u=x+d; fu=(*f)(u); } else { u=x+SIGN(tol1,d); fu=(*f)(u); if (fu > fx) { *xmin=x; return fx; } } du=(*df)(u); if (fu <= fx) { if (u >= x) a=x; else b=x; MOV3(v,fv,dv, w,fw,dw) MOV3(w,fw,dw, x,fx,dx) MOV3(x,fx,dx, u,fu,du) } else { if (u < x) a=u; else b=u; if (fu <= fw || w == x) { MOV3(v,fv,dv, w,fw,dw) MOV3(w,fw,dw, u,fu,du) } else if (fu < fv || v == x || v == w) { MOV3(v,fv,dv, u,fu,du) } } } nrerror("Too many iterations in routine dbrent"); return 0.0; } ``` #### 5. 多维下山单纯形法 从这部分开始，我们将考虑多维函数的最小化问题。下山单纯形法由 Nelder 和 Mead 提出，该方法只需要函数评估，不需要导数信息。虽然它在函数评估次数方面效率不高，但在计算负担较小的问题中，它可能是快速解决问题的最佳方法。 在 \(N\) 维空间中，单纯形是由 \(N + 1\) 个点（或顶点）及其所有相互连接的线段、多边形面等组成的几何图形。例如，在二维空间中，单纯形是一个三角形；在三维空间中，它是一个四面体。 下山单纯形法需要从 \(N + 1\) 个点开始，定义一个初始单纯形。通常，我们可以选择一个初始点 \(P_0\)，然后将其他 \(N\) 个点定义为 \(P_i = P_0 + \lambda e_i\)，其中 \(e_i\) 是 \(N\) 个单位向量，\(\lambda\) 是我们对问题特征长度尺度的猜测。 该方法通过一系列步骤来寻找最小值。大多数步骤是将单纯形中函数值最大的点（“最高点”）通过对面移动到一个较低的点，这种操作称为反射。当条件允许时，方法会在一个或多个方向上扩展单纯形以采取更大的步骤；当到达“谷底”时，方法会在横向方向上收缩并尝试沿着山谷向下移动；如果单纯形试图“穿过针眼”，它会在所有方向上收缩，围绕其最低点（最佳点）收紧。 以下是下山单纯形法的代码： ```c #include <math.h> #include "nrutil.h" #define TINY 1.0e-10 #define NMAX 5000 #define GET_PSUM \ for (j=1;j<=ndim;j++) {\ for (sum=0.0,i=1;i<=mpts;i++) sum += p[i][j];\ psum[j]=sum;} #define SWAP(a,b) {swap=(a);(a)=(b);(b)=swap;} void amoeba(float **p, float y[], int ndim, float ftol, float (*funk)(float []), int *nfunk) { float amotry(float **p, float y[], float psum[], int ndim, float (*funk)(float []), i ```