计算物理基础-浅谈函数极值和局域优化
我们知道，除了微分方程，不少物理规律都可以用极值的形式来描述，比如光传播时的光程极值原理、运动体系的最小作用量原理、材料计算中的密度泛函原理、稳定结构的自由能极小等。而在这些问题中，我们都需要面对多变量问题的极值求解。
在数值计算中，优化相关的内容非常丰富，这里先介绍常见的局域优化，想了解更多的朋友可以参考《最优化导论》，作者: Edwin K. P. Chong / Stanislaw H. Zak，An Introduction to Optimization,Foulth Edition，译者: 孙志强 / 白圣建 / 郑永斌 / 刘伟。
（1） 函数极值和解方程（组）的关联
对于多变量函数，函数极值和之前提到的解方程（组）是互相关联的。因此，函数极值的计算（尤其是局域优化的思路）通常和初始值的选择有关，而程序通常在梯度达到阈值时停止。如果要得到更好的极值，通常需要全局搜索的算法（比如模拟退火、遗传-进化算法、微粒群算法等），这个在后续的内容再介绍。
（2） 库函数fminbnd和fminunc的使用
考虑一元函数 cos(3x).exp(-x)-1)在[-2,4]区间，从图像曲线看，有3个极小值点。在定义好目标函数fun = @(x) cos(3*x).*exp(-x)-1 后，如果用fminbnd函数，语法是fminbnd(fun,-2,4)，这里的-2和4就是指定的区间，返回结果是 3.0343，是其中一个极小值对应的自变量位置。如果用fminunc函数，fminunc(fun,-2) 这里的-2是初始点，取-2，1，4等作为初始点，分别得到极值位置-1.1544，0.9399，3.0343，从图像上也是很好理解的，就是初始点接近，而且中间梯度单调，自然落入附近的极值，这是局域优化常见的现象。
此外，fminunc 可以用于多变量函数。比如函数 fun = @(x) x(1).^2+x(2).^2-x(1)-x(2); 用 fminunc (fun, [1 1]) 计算极值，这里[1 1]是x_1和x_2的初值，结果返回 ： 0.5000 0.5000，和我们解析推导结果一致。如果处理带约束的问题，需要用fmincon函数，在matlab里是自带的，在octave中需要通过 pkg load optim 提前加载。
（3） 牛顿法、拟牛顿法
假定在点x_k附近进行二阶泰勒展开，得到q(x)是原函数f(x)的近似，进一步对q(x)求导，可以得到x_{k+1}和x_k的迭代公式。
注意，拓展到高维时，f ' 是梯度，f '' 对应海森矩阵。和之前的解方程是类似的。由于该方法每次要计算海森矩阵，20世纪60年代，Davidon-Fletcher-Powell提出DFP方法来简化，70年代，Broyden-Fletcher-Goldfard-Shanno等人提出改进方法，这类方法统称拟牛顿法。
此外，牛顿法对体系的二阶导数有要求，如果类似抛物线等情况，迭代很快指向体系的平衡位置，但有时会偏离平衡位置越来越远（有兴趣的朋友可以测试一下李纳-琼斯势的情况）。
为求函数（单变量、多变量）极值，在局域优化方法中，关键的是搜索方向和移动步长。在牛顿法中，可以认为方向和步长都乘在一起了，该方法最大的问题就是每次要求二阶导数。最速下降法的搜索方向是负梯度方向，步长通过一元函数极值来确定，假设前进步长是alpha，新的位置是当前位置加上alpha乘以负梯度，求一元函数的极小值确定alpha，然后得到新的位置；之后不断迭代，在达到局域极小附近，梯度向量对应的模逐渐减小，可以通过设置阈值使程序停止迭代。
附录程序是用最速下降法求函数的最小值。这里选择初始x=2.5,计算梯度gx和函数值y，如果gx向量的模没有达到阈值，进入循环，这里每次步长alpha从2开始，或者更大的数，之后不断减半，直到对应的函数值小于当前数值。这里当前函数值是y，更新后是ny。在初始点附近的3个点做插值（这3个点满足第二点比另外两点数值小），找到使函数值最小的步长，更新自变量x，然后判断该点梯度是否达到阈值。
在高维问题中（上图右边是二维情况），最速下降法在接近极小值附近会出现锯齿形的迭代路径，影响计算效率。共轭梯度法可以很好解决该问题，搜索方向在梯度方向上加了修正，保证每次搜到方向关于海森矩阵共轭；通过Hestenes-Stiefel公式消去迭代公式中的海森矩阵，可以使效率极大提高。此外，对于最速下降法，沿着梯度做一维搜索步长也是耗时的，可以根据梯度分量的大小适当变步长即可，在多数问题时，评估目标函数是很慢的。
附录：
A. 牛顿法求函数极值
clear
x(1)=2;
for i=2:10
f1=4*x(i-1)^3-3;
f2=12*x(i-1)^2;
x(i)=x(i-1)-f1/f2;
end
y=x.^4-3*x-10;
plot(x,y,'o-')
hold on
clear
x=linspace(-2,2);
y=x.^4-3*x-10;
plot(x,y)
B. 最速下降法求函数极值
clear
x=0.8;y=1/x^12-2/x^6;i=1;gx=-12/x^13+12/x^7;
f(i)=1/x^12-2/x^6;
d(i)=x;
while abs(gx)>0.001&&i<100 %梯度向量模的阈值
alpha=0.01;
tp=[];ys=[];
tp(1)=0;ys=y;
x=x-gx*alpha; %计算更新后的x和y
gx=-12/x^13+12/x^7;
i=i+1;
f(i)=1/x^12-2/x^6;
d(i)=x;
end
plot(d,f,'o-')
hold on
clear
x=linspace(0.8,2);
y=1./x.^12-2./x.^6;
plot(x,y)
