函数优化方法详解与实践

# 函数优化方法详解与实践 ## 1 优化问题概述 在优化问题中，起始位置的选择至关重要。若起始点处于包含最小值的感兴趣区域之外，其斜率（一阶导数）和曲率（二阶导数）接近零，函数值近乎恒定。例如： ```mathematica D[function,x]/.x -> xStart 0.0000214326 D[function,{x,2}]/.x -> xStart −0.0000250037 ``` 这种情况下，迭代算法难以找到通往最小值的路径，搜索算法可能陷入停滞，无法收敛到最小值。而且在实际中，优化失败时往往难以找出问题所在。虽然有众多参数可用于调整局部和全局优化方法，但这并不能保证总能解决所有问题。此时，基于理论或实践经验的最优位置先验知识就显得尤为关键。 ## 2 迭代局部优化 ### 2.1 基本原理 迭代局部优化（求函数最大值等同于求 -f 或 f -1 的最小值）原则上很简单，即从给定起始位置，通过适当步骤尽可能快地向山下移动，直到在所需精度内达到局部最小值。不同的局部优化方法主要区别在于评估的函数信息数量，以此来确定沿所选下山方向的步长。评估部分决定每次迭代的计算成本，方向部分决定向局部最小值的收敛速度，这两部分通常相互制约：评估的函数信息越多，单次迭代执行越慢，但由于步长和方向更合适，迭代步数可能减少。 ### 2.2 常见方法 - **仅使用函数值评估的方法**：如单纯形法，仅利用不同位置的函数值来识别或智能或简单的下山路径，并采用自适应步长。 - **使用一阶导数信息的方法**：包括梯度法、共轭梯度法和拟牛顿法。这些方法除了函数值外，还使用（一阶导数）斜率/梯度信息，实现最速下降方向。共轭梯度法和拟牛顿法最多用 M 步就能找到 M 维抛物超曲面的（唯一且全局）最小值。 - **使用二阶导数信息的方法**：牛顿法利用（二阶导数）曲率信息，在局部最小值附近实现更快收敛。对于抛物超曲面，牛顿法一步就能直接到达最小值。不过，牛顿法每次迭代需要评估大量函数信息，耗时较长，因此仅在最小值附近有效。 ### 2.3 特殊函数的优化方法 对于特殊类型的函数，如平方和函数，Levenberg - Marquardt 等组合方法很有帮助，它们能在远离最小值时采用梯度步，接近最小值时采用牛顿步。 ### 2.4 方法比较与选择 一般来说，没有绝对最好的迭代局部优化方法。对于许多实际最小化问题，共轭梯度法和拟牛顿法在计算成本（函数和一阶导数评估）和局部最小值收敛速度之间取得了较好的平衡。 ### 2.5 实践示例 以下是使用 Mathematica 进行局部优化的示例： ```mathematica function=1.0-Cos[x]/(1.0+0.01*x^2); pureFunction=Function[argument,function/.x -> argument]; argumentRange={-10.0,10.0}; functionValueRange={-0.2,2.2}; startPosition=8.0; startPoint={startPosition,function/.x -> startPosition}; points2D={startPoint}; labels={"x","y","Function with multiple optima"}; CIP‘Graphics‘Plot2dPointsAboveFunction[points2D,pureFunction,labels, GraphicsOptionArgumentRange2D -> argumentRange, GraphicsOptionFunctionValueRange2D -> functionValueRange] localMinimum=FindMinimum[function,{x,startPosition}] {0.28015,{x →6.19389}} minimumPoint={x/.localMinimum[[2]],localMinimum[[1]]}; points2D={startPoint,minimumPoint}; labels，{"x","y","Local minimization"}; arrowGraphics=Graphics[{Thick,Red,{Arrowheads[Medium], Arrow[{startPoint,minimumPoint}]}}]; functionGraphics=CIP‘Graphics‘Plot2dPointsAboveFunction[points2D, pureFunction,labels, GraphicsOptionArgumentRange2D -> argumentRange, GraphicsOptionFunctionValueRange2D -> functionValueRange]; Show[functionGraphics,arrowGraphics] ``` 从不同的起始位置开始，可能会找到不同的最小值： ```mathematica startPosition=2.0; localMinimum=FindMinimum[function,{x,startPosition}] {0.,{x →9.64816×10−12}} startPoint={startPosition,function/.x -> startPosition}; minimumPoint={x/.localMinimum[[2]],localMinimum[[1]]}; points2D={startPoint,minimumPoint}; arrowGraphics=Graphics[{Thick,Red,{Arrowheads[Medium], Arrow[{startPoint,minimumPoint}]}}]; functionGraphics=CIP‘Graphics‘Plot2dPointsAboveFunction[points2D, pureFunction,labels, GraphicsOptionArgumentRange2D -> argumentRange, GraphicsOptionFunctionValueRange2D -> functionValueRange]; Show[functionGraphics,arrowGraphics] ``` 在某些情况下，近似最小值可能恰好是全局最小值，但一般来说，局部优化通常只能找到局部最优值。 ### 2.6 局部优化方法总结 | 方法类型 | 优点 | 缺点 | 适用场景 | | ---- | ---- | ---- | ---- | | 仅使用函数值评估的方法（如单纯形法） | 实现简单，无需导数计算 | 收敛速度慢 | 函数导数难以计算或导数不存在的情况 | | 使用一阶导数信息的方法（梯度法、共轭梯度法、拟牛顿法） | 收敛速度较快，能利用梯度信息 | 需要计算导数 | 函数导数容易计算的情况 | | 使用二阶导数信息的方法（牛顿法） | 收敛速度快，在最小值附近效果好 | 计算量大，需要二阶导数信息 | 接近最小值且函数二阶导数容易计算的情况 | ### 2.7 局部优化流程 ```mermaid graph TD; A[选择起始位置] --> B[评估函数信息]; B --> C[确定步长和方向]; C --> D[更新当前位置]; D --> E{是否达到收敛条件}; E -- 是 --> F[输出局部最小值]; E -- 否 --> B; ``` ## 3 迭代全局优化 ### 3.1 全局优化的目标与挑战 优化函数通常旨在找到科学相关参数空间的全局最优值。迭代局部搜索虽有可能找到全局最优值，但往往会陷入起始位置附近的局部最优值。全局优化策略试图通过对先验定义的整个搜索空间进行采样来解决这个问题。这需要为待全局优化函数 \( f(x_1,x_2,\cdots,x_M) \) 的每个参数 \( x_1,x_2,\cdots,x_M \) 设定一组最小/最大值，假定全局最优值位于由这些 \( M \) 个最小/最大区间 \( [x_{1,min}, x_{1,max}] \) 到 \( [x_{M,min}, x_{M,max}] \) 所构成的搜索空间内。 ### 3.2 系统网格搜索 最直接的全局优化方法是系统网格搜索，即在预先定义的参数搜索空间内，对等间距网格点的函数值进行评估，然后相互比较以找出最优值。以下是一个示例，用于近似之前绘制的曲面 \( f(x,y) \) 的全局最大值： ```mathematica function=1.9*(1.35+Exp[x]*Sin[13.0*(x-0.6)^2]*Exp[-y]* Sin[7.0*y]); pureFunction= Function[{argument1,argument2}, function/.{x -> argument1,y -> argument2}]; xMinBorderOfSearchSpace=0.0; xMaxBorderOfSearchSpace=1.0; yMinBorderOfSearchSpace=0.0; yMaxBorderOfSearchSpace=1.0; numberOfGridPointsPerDimension=10.0; gridPoints3D={}; Do[ Do[ AppendTo[gridPoints3D,{x,y,0.0}], {x,xMinBorderOfSearchSpace,xMaxBorderOfSearchSpace, (xMaxBorderOfSearchSpace-xMinBorderOfSearchSpace)/ (numberOfGridPointsPerDimension-1.0)} ], {y,yMinBorderOfSearchSpace,yMaxBorderOfSearchSpace, (yMaxBorderOfSearchSpace-yMinBorderOfSearchSpace)/ (numberOfGridPointsPerDimension-1.0)} ]; xRange={-0.1,1.1}; yRange={-0.1,1.1}; labels={"x","y","z"}; viewPoint3D={3.5,-2.4,1.8}; CIP‘Graphics‘Plot3dPointsWithFunction[gridPoints3D,pureFunction, labels, GraphicsOptionArgument1Range3D -> xRange, GraphicsOptionArgument2Range3D -> yRange, GraphicsOptionViewPoint3D -> viewPoint3D] ``` 接着，评估这些网格点的函数值并进行比较： ```mathematica winnerGridPoint3D={}; maximumFunctionValue=-Infinity; Do[ functionValue=pureFunction[gridPoints3D[[i, 1]], gridPoints3D[[i, 2]]]; If[functionValue>maximumFunctionValue, maximumFunctionValue=functionValue; winnerGridPoint3D={gridPoints3D[[i, 1]],gridPoints3D[[i, 2]], maximumFunctionValue} ], {i,Length[gridPoints3D]} ]; ``` 得到获胜的网格点和对应的最大函数值： ```mathematica winnerGridPoint3D {1.,0.222222,6.17551} maximumFunctionValue 6.17551 ``` 可以对结果进行可视化： ```mathematica Do[ If[gridPoints3D[[i,1]] == winnerGridPoint3D[[1]] && gridPoints3D[[i,2]] == winnerGridPoint3D[[2]], gridPoints3D[[i]] = winnerGridPoint3D ], {i,Length[gridPoints3D]} ]; arrowStartPoint={winnerGridPoint3D[[1]],winnerGridPoint3D[[2]],0.0}; arrowGraphics3D=Graphics3D[{Thick,Red,{Arrowheads[Medium], Arrow[{arrowStartPoint,winnerGridPoint3D}]}}]; plotStyle3D=Directive[Green,Specularity[White,40],Opacity[0.4]]; functionGraphics3D=CIP‘Graphics‘Plot3dPointsWithFunction[ gridPoints3D,pureFunction,labels, GraphicsOptionArgument1Range3D -> xRange, GraphicsOptionArgument2Range3D -> yRange, GraphicsOptionViewPoint3D -> viewPoint3D, GraphicsOptionPlotStyle3D -> plotStyle3D]; Show[functionGraphics3D,arrowGraphics3D] ``` 系统网格搜索得到的全局最优值只是近似值，误差与定义的网格间距有关。为了细化这个近似值，可以将其作为后续局部搜索的起点： ```mathematica globalMaximum=FindMaximum[function,{{x,winnerGridPoint3D[[1]]}, {y,winnerGridPoint3D[[2]]}}] {6.54443,{x →0.959215,y →0.204128}} ``` 最后对局部细化的结果进行可视化： ```mathematica globalMaximumPoint3D={x/.globalMaximum[[2,1]], y/.globalMaximum[[2,2]],globalMaximum[[1]]}; xRange={0.90,1.005}; yRange={0.145,0.26}; arrowGraphics3D=Graphics3D[{Thick,Red,{Arrowheads[Medium], Arrow[{winnerGridPoint3D,globalMaximumPoint3D}]}}]; points3D={winnerGridPoint3D,globalMaximumPoint3D}; functionGraphics3D=CIP‘Graphics‘Plot3dPointsWithFunction[points3D, pureFunction,labels, GraphicsOptionArgument1Range3D -> xRange, GraphicsOptionArgument2Range3D -> yRange, GraphicsOptionViewPoint3D -> viewPoint3D, GraphicsOptionPlotStyle3D -> plotStyle3D]; Show[functionGraphics3D,arrowGraphics3D] ``` 系统网格搜索仅适用于低维网格的全局优化问题，因为随着参数数量的增加，需要评估的网格点数量会呈指数级增长。例如，对于 12 个参数，每个参数有 10 个网格点，网格将包含一万亿（\( 10^{12} \)）个点，函数值评估的计算量会变得非常大。 ### 3.3 随机搜索 作为替代方案，可以将搜索空间中要测试的参数值数量限制在可管理的范围内。随机选择测试点是一种合理的选择，因为事先对搜索空间的任何偏好部分都没有先验知识。以下是一个使用 20 个随机测试点的示例： ```mathematica SeedRandom[1]; randomPoints3D= Table[ {RandomReal[{xMinBorderOfSearchSpace,xMaxBorderOfSearchSpace}], RandomReal[{yMinBorderOfSearchSpace,yMaxBorderOfSearchSpace}], 0.0}, {20} ]; CIP‘Graphics‘Plot3dPointsWithFunction[randomPoints3D,pureFunction, labels, GraphicsOptionArgument1Range3D -> xRange, GraphicsOptionArgument2Range3D -> yRange, GraphicsOptionViewPoint3D -> viewPoint3D] ``` 评估随机点的函数值并找出获胜点： ```mathematica winnerRandomPoint3D={}; maximumFunctionValue=-Infinity; Do[ functionValue=pureFunction[randomPoints3D[[i, 1]], randomPoints3D[[i, 2]]]; If[functionValue>maximumFunctionValue, maximumFunctionValue=functionValue; winnerRandomPoint3D={randomPoints3D[[i, 1]],randomPoints3D[[i, 2]], maximumFunctionValue} ], {i,Length[randomPoints3D]} ]; ``` 对结果进行可视化： ```mathematica Do[ If[randomPoints3D[[i,1]] == winnerRandomPoint3D[[1]] && randomPoints3D[[i,2]] == winnerRandomPoint3D[[2]], randomPoints3D[[i]] = winnerRandomPoint3D ], {i,Length[randomPoints3D]} ]; arrowStartPoint={winnerRandomPoint3D[[1]],winnerRandomPoint3D[[2]], 0.0}; arrowGraphics3D=Graphics3D[{Thick,Red,{Arrowheads[Medium], Arrow[{arrowStartPoint,winnerRandomPoint3D}]}}]; plotStyle3D=Directive[Green,Specularity[White,40],Opacity[0.4]]; functionGraphics3D=CIP‘Graphics‘Plot3dPointsWithFunction[ randomPoints3D,pureFunction,labels, GraphicsOptionArgument1Range3D -> xRange, GraphicsOptionArgument2Range3D -> yRange, GraphicsOptionViewPoint3D -> viewPoint3D, GraphicsOptionPlotStyle3D -> plotStyle3D]; Show[functionGraphics3D,arrowGraphics3D] ``` 如果对这个全局优化结果进行后续局部最大值搜索： ```mathematica globalMaximum=FindMaximum[function, {{x,winnerRandomPoint3D[[1]]},{y,winnerRandomPoint3D[[2]]}}] {4.55146,{x →0.265291,y →0.204128}} ``` 可能只能找到局部最大值，而错过全局最大值。这是因为随机测试点数量不足，增加随机测试点数量可以提高找到接近全局最大值测试点的概率，但随着参数数量（维度）的增加，随机搜索也会变得非常困难。 ### 3.4 进化算法 在面对大规模搜索空间的困难搜索问题时，进化算法应运而生。进化算法以基本随机的方式运行，类似于纯随机搜索，但还借鉴了生物进化中经过验证的优化策略，包括变异（随机变化）、交叉或重组（一种随机混合，引导向有前景的搜索空间区域跳跃）以及适者生存（放大到目前为止找到的最优解）。进化周期试图通过逐步组合最优解的部分（模式）来加速向全局最优值的搜索。Mathematica 通过 `NMinimize` 和 `NMaximize` 命令提供了基于进化算法的全局优化程序，使用 `DifferentialEvolution` 方法选项。以下是一个使用进化算法进行全局最大值搜索的示例： ```mathematica globalMaximum=NMaximize[{function, {xMinBorderOfSearchSpace<x<xMaxBorderOfSearchSpace, yMinBorderOfSearchSpace<y<yMaxBorderOfSearchSpace}}, {x,y}, Method -> {"DifferentialEvolution","PostProcess" -> False}] {6.54443,{x →0.959215,y →0.204128}} ``` 需要注意的是，虽然进化算法很受欢迎，但它们属于最后的手段，因为计算成本可能极高，非常耗时。以下是另一个使用进化算法对具有多个最优值的函数进行全局最小值搜索的示例： ```mathematica function=1.0-Cos[x]/(1.0+0.01*x^2); pureFunction=Function[argument,function/.x -> argument]; xMinBorderOfSearchSpace=-10.0; xMaxBorderOfSearchSpace=15.0; globalMinimum=NMinimize[{function, xMinBorderOfSearchSpace<x<xMaxBorderOfSearchSpace},x, Method -> {"DifferentialEvolution","PostProcess" -> False}] {5.16341×10−10,{x →−0.0000318188}} ``` 然而，如果搜索空间选择不当，可能无法在默认的最大迭代次数内找到全局最小值： ```mathematica xMinBorderOfSearchSpace=50.0; xMaxBorderOfSearchSpace=60.0; globalMinimum=NMinimize[{function, xMinBorderOfSearchSpace<x<xMaxBorderOfSearchSpace},x, Method -> {"DifferentialEvolution","PostProcess" -> False}] {0.9619,{x →50.2272}} xMinBorderOfSearchSpace=-100000.0; xMaxBorderOfSearchSpace=100000.0; globalMinimum=NMinimize[{function, xMinBorderOfSearchSpace<x<xMaxBorderOfSearchSpace},x, Method -> {"DifferentialEvolution","PostProcess" -> False}] {0.805681,{x →19.2638}} ``` ### 3.5 全局优化方法比较 | 方法 | 优点 | 缺点 | 适用场景 | | ---- | ---- | ---- | ---- | | 系统网格搜索 | 全面搜索，结果准确 | 计算量大，适用于低维问题 | 低维全局优化问题 | | 随机搜索 | 实现简单，可控制计算量 | 可能错过全局最优值 | 对计算资源有限且对精度要求不是极高的情况 | | 进化算法 | 能处理大规模搜索空间 | 计算成本高，耗时 | 复杂的高维全局优化问题 | ### 3.6 全局优化流程 ```mermaid graph TD; A[定义搜索空间] --> B[选择搜索方法]; B --> C[生成测试点]; C --> D[评估函数值]; D --> E[选择最优解]; E --> F{是否满足终止条件}; F -- 是 --> G[输出全局最优值]; F -- 否 --> C; ``` ## 4 约束迭代优化 ### 4.1 约束优化的概念 在之前的全局优化示例中，已经涉及到约束优化，因为预先定义的搜索空间是搜索的一个约束。一般来说，如果优化问题没有任何额外限制，则称为无约束优化；如果优化受到一个或多个约束，则进入约束优化领域。 ### 4.2 约束优化示例 以下是一个对函数进行全局最小化，同时将 \( x \) 值限制在定义区间内的示例： ```mathematica function=1.0-Cos[x]/(1.0+0.01*x^2); pureFunction=Function[argument,function/.x -> argument]; xMinConstraint=2.0; xMaxConstraint=11.0; constraint=xMinConstraint<x<xMaxConstraint; constrainedGlobalMinimum=NMinimize[{function,constraint},x, Method -> {"DifferentialEvolution","PostProcess" -> False}] {0.28015,{x →6.19386}} ``` 对约束全局最小值进行可视化： ```mathematica constrainedMinimumPoint={x/.constrainedGlobalMinimum[[2]], constrainedGlobalMinimum[[1]]}; points2D={constrainedMinimumPoint}; argumentRange={-12.0,17.0}; functionValueRange={-0.2,2.2}; labels={"x","y","Constrained global minimization"}; functionGraphics=CIP‘Graphics‘Plot2dPointsAboveFunction[points2D, pureFunction,labels, GraphicsOptionArgumentRange2D -> argumentRange, GraphicsOptionFunctionValueRange2D -> functionValueRange]; constraintGraphics=Graphics[{RGBColor[1,0,0,0.1], Rectangle[{xMinConstraint,functionValueRange[[1]]}, {xMaxConstraint,functionValueRange[[2]]}]}]; Show[functionGraphics,constraintGraphics] ``` 无约束和约束全局最优值可能不同，约束全局最优值甚至可能不是无约束优化问题的最优值。以下是一个来自 Mathematica 教程的 3D 表面示例： ```mathematica function= -1.0/((x+1.0)^2+(y+2.0)^2+1)-2.0/((x-1.0)^2+(y-1.0)^2+1)+2.0; pureFunction=Function[{argument1,argument2}, function/.{x -> argument1,y -> argument2}]; xRange={-3.0,3.0}; yRange={-3.0,3.0}; labels={"x","y","z"}; CIP‘Graphics‘Plot3dFunction[pureFunction,xRange,yRange,labels] ``` 这个表面包含一个局部最小值和一个全局最小值。使用 `FindMinimum` 命令进行迭代局部最小值搜索，根据不同的起始位置可以找到不同的最小值： ```mathematica startPosition={-2.5,-1.5}; localMinimum=FindMinimum[function,{{x,startPosition[[1]]}, {y,startPosition[[2]]}}] {0.855748,{x →−0.978937,y →−1.96841}} startPoint={startPosition[[1]],startPosition[[2]], function/.{x -> startPosition[[1]],y -> startPosition[[2]]}}; minimumPoint={x/.localMinimum[[2,1]],y/.localMinimum[[2,2]], localMinimum[[1]]}; points3D={startPoint,minimumPoint}; arrowGraphics=Graphics3D[{Thick,Red,{Arrowheads[Medium], Arrow[{startPoint,minimumPoint}]}}]; plotStyle3D=Directive[Green,Specularity[White,40],Opacity[0.4]]; functionGraphics=CIP‘Graphics‘Plot3dPointsWithFunction[points3D, pureFunction,labels, GraphicsOptionArgument1Range3D -> xRange, GraphicsOptionArgument2Range3D -> yRange, GraphicsOptionPlotStyle3D -> plotStyle3D]; Show[functionGraphics,arrowGraphics] startPosition={-0.5,2.5}; localMinimum=FindMinimum[function,{{x,startPosition[[1]]}, {y,startPosition[[2]]}}] {−0.071599,{x →0.994861,y →0.992292}} startPoint={startPosition[[1]],startPosition[[2]], function/.{x -> startPosition[[1]],y -> startPosition[[2]]}}; minimumPoint={x/.localMinimum[[2,1]],y/.localMinimum[[2,2]], localMinimum[[1]]}; points3D={startPoint,minimumPoint}; arrowGraphics=Graphics3D[{Thick,Red,{Arrowheads[Medium], Arrow[{startPoint,minimumPoint}]}}]; functionGraphics=CIP‘Graphics‘Plot3dPointsWithFunction[points3D, pureFunction,labels, GraphicsOptionArgument1Range3D -> xRange, GraphicsOptionArgument2Range3D -> yRange, GraphicsOptionPlotStyle3D -> plotStyle3D]; Show[functionGraphics,arrowGraphics] ``` ### 4.3 约束优化流程 ```mermaid graph TD; A[定义目标函数和约束条件] --> B[选择优化方法]; B --> C[初始化参数]; C --> D[迭代更新参数]; D --> E{是否满足收敛条件}; E -- 是 --> F[输出约束最优值]; E -- 否 --> D; ``` 综上所述，不同的优化方法适用于不同的场景，在实际应用中需要根据问题的特点选择合适的方法。通过对局部优化、全局优化和约束优化的深入了解，可以更好地解决各种复杂的优化问题。