## 无约束最优化方法
最优化(Optimization)是工程技术、经济管理、 科学研究、社会生活中经常遇到的问题, 如：结构设计，资源分配，生产计划，运输方案等等。
最优化: 在一定条件下，寻求使目标最大(小)的决策解决优化问题的手段
• 经验积累，主观判断
• 作试验，比优劣
• 建立数学模型，求解最优策略
### 优化问题的数学模型
无约束优化问题标准形式：
求n维决策变量x=[x1,x2,…,xn]T在可行域Ω中
使得目标函数f(x)→min,x⊆Ω⊆Rn
⇒
**可行解**只满足(2)，**最优解**满足(1)(2)
**无约束优化**只有(1)，**约束优化**需要(1)(2)
### 无约束优化的基本方法
### 求解无约束优化的基本思路：[下山法](http://blog.csdn.net/u013007900/article/details/43525733)
#### 基本思想：
和爬山法基本思想一样，在Rn中某一点，确定一个**搜索放向**然后沿着该方向的**移动步长**，得到使目标函数下降的最新点。
#### 迭代步骤：
Step 1 初始化：初始点x0，终止准则等
Step 2 迭代改进：方向dk，步长αk
xk+1=xk+(αd)k 且 f(xk+1)<f(xk)
Step 3 终止检验：得到近优解或k+1⇒k转2> 选择dk,αk使得f下降更快⇒利用不同的算法
### 最速下降法（梯度法）
泰勒级数：
其中Δx可正可负，但必须充分接近于0。
X沿D方向移动步长α后，变为X+αD。由泰勒展开式：
目标函数：
α确定的情况下即最小化：
向量的内积何时最小？当然是两向量方向相反时。所以X移动的方向应该和梯度的方向相反。
**下降方向：**∇fT(xk)dk<0
**最速下降方向：**dk=−∇f(xk)负梯度方向
**迭代改进形式：**xk+1=xk−∇f(xk)
算法特点：在初始状态改进较快，在最优解附近改进缓慢。
求函数f(x1,x2)=x21+4x22在极小点，以初始点X0=(1,1)T。
~~~
#include"matrix.h"//大家自己YY一个矩阵类吧
#include<iostream>
#include<iomanip>
#include<cmath>
#include<limits>
#include<cassert>
using namespace std;
const int SIZE=2;
const double ZETA=0.001;
inline double func(Matrix<double> &X){
assert(SIZE==X.getRows());
double x1=X.get(0,0);
double x2=X.get(1,0);
return x1*x1+4*x2*x2;
}
inline Matrix<double> gradient(Matrix<double> &X){
assert(SIZE==X.getRows());
double x1=X.get(0,0);
double x2=X.get(1,0);
Matrix<double> rect(SIZE,1);
rect.put(0,0,2*x1);
rect.put(1,0,8*x2);
return rect;
}
inline Matrix<double> Hesse(Matrix<double> &X){
Matrix<double> rect(SIZE,SIZE);
rect.put(0,0,2);
rect.put(0,1,0);
rect.put(1,0,0);
rect.put(1,1,8);
return rect;
}
int main(int agrc,char *argv[]){
Matrix<double> X(SIZE,1);
X.put(0,0,1);
X.put(1,0,1);
int iteration=0;
double value=func(X);
double newValue=numeric_limits<double>::max();
while(++iteration){
Matrix<double> G=gradient(X);
double factor=((G.getTranspose()*G).get(0,0))/((G.getTranspose()*Hesse(X)*G).get(0,0));
for(int i=0;i<G.getRows();++i)
for(int j=0;j<G.getColumns();++j)
G.put(i,j,G.get(i,j)*factor);
Matrix<double> newX=X-G;
//cout<<"X=["<<newX.get(0,0)<<","<<newX.get(1,0)<<"]"<<endl;
newValue=func(newX);
if(fabs(newValue-value)<ZETA)
break;
else{
X=newX;
value=newValue;
}
//cout<<"本次迭代找到的极小值"<<value<<endl;
}
cout<<"迭代次数"<<iteration<<endl;
cout<<"极小值"<<value<<endl;
return 0;
}
~~~
> 这份算法使用的是利用Hesse矩阵和泰勒公式计算最优步长
若f(X)具有二阶连续偏导，由泰勒展开式可得：
H是f(X)的Hesse矩阵。
可得最优步长：
g是f(X)的梯度矩阵。
此时：
可见最速下降法中最优步长不仅与梯度有关，而且与Hesse矩阵有关。
### 牛顿迭代法
将f(xk+1)在xk点作泰勒展开至二阶项，用d替代dk，有
f(xk+1)=f(xk+d)=f(xk)+∇fT(xk)d+12dT∇2f(xk)d
求d使得f(xk+1)极小⇒右端对d的导数为0⇒∇f(xk)+∇2f(xk)d=0
可得到> **牛顿方程**∇2f(xk)d = −∇f(xk)
**牛顿方向**dk = −(∇2f(xk))−1∇f(xk)
**迭代格式**xk+1 = xk−(∇2f(xk))−1∇f(xk)
其中∇2f(xk)为[Hessen矩阵](http://baike.baidu.com/link?url=E2wdH1BhTKFos1JYd_NhHvkmGAyABsHJDIgSyygvzTuoAs41tvb4OGzzWThbUMwobdBAEPpS7nzia-GQUHHf7a)
可见x的搜索方向是−(∇2f(xk))−1∇f(xk)，函数值在此方向上下降，就需要它与梯度方向相反，即−∇f(xk)(∇2f(xk))−1∇f(xk)<0，所以要求在每一个迭代点上Hessian矩阵必须是正定的。
算法特点：
牛顿法是二次收敛的，并且收敛阶数是2。一般目标函数在最优点附近呈现为二次函数，于是可以想像最优点附近用牛顿迭代法收敛是比较快的。而在开始搜索的几步，我们用梯度下降法收敛是比较快的。将两个方法融合起来可以达到满意的效果。
收敛快是牛顿迭代法最大的优点，但也有致命的缺点：Hessian矩阵及其逆的求解计算量大，更何况在某个迭代点Xk处Hessian矩阵的逆可能根本就不存在（即Hessian
矩阵奇异），这样无法求得Xk+1。
求f(x1,x2)=(x1−2)4+(x1−2x2)2的极小点，初始点取X=(0,3)T
~~~
#include"matrix.h"
#include<iostream>
#include<iomanip>
#include<cmath>
#include<limits>
#include<cassert>
using namespace std;
const int SIZE=2;
const double ZETA=0.001;
inline double func(Matrix<double> &X){
assert(SIZE==X.getRows());
double x1=X.get(0,0);
double x2=X.get(1,0);
return pow(x1-2,4)+pow(x1-2*x2,2);
}
inline Matrix<double> gradient(Matrix<double> &X){
assert(SIZE==X.getRows());
double x1=X.get(0,0);
double x2=X.get(1,0);
Matrix<double> rect(SIZE,1);
rect.put(0,0,4*pow(x1-2,3)+2*(x1-2*x2));
rect.put(1,0,-4*(x1-2*x2));
return rect;
}
inline Matrix<double> Hesse(Matrix<double> &X){
Matrix<double> rect(SIZE,SIZE);
double x1=X.get(0,0);
double x2=X.get(1,0);
rect.put(0,0,12*pow(x1-2,2)+2);
rect.put(0,1,-4);
rect.put(1,0,-4);
rect.put(1,1,8);
return rect;
}
int main(int agrc,char *argv[]){
Matrix<double> X(SIZE,1);
X.put(0,0,0);
X.put(1,0,3);
int iteration=0;
double value=func(X);
double newValue=numeric_limits<double>::max();
while(++iteration){
Matrix<double> G=gradient(X);
Matrix<double> H=Hesse(X);
Matrix<double> newX=X-H.getInverse()*G;
cout<<"X=["<<newX.get(0,0)<<","<<newX.get(1,0)<<"]"<<endl;
newValue=func(newX);
if(fabs(newValue-value)<ZETA)
break;
else{
X=newX;
value=newValue;
}
cout<<"本次迭代找到的极小值"<<value<<endl;
}
cout<<"迭代次数"<<iteration<<endl;
cout<<"极小值"<<value<<endl;
return 0;
}
~~~
### 拟牛顿法
Hessian矩阵在拟牛顿法中是不计算的，拟牛顿法是构造与Hessian矩阵相似的正定矩阵，这个构造方法，使用了目标函数的梯度（一阶导数）信息和两个点的“位移”（Xk−Xk−1）来实现。
有人会说，是不是用Hessian矩阵的近似矩阵来代替Hessian矩阵，会导致求解效果变差呢？事实上，效果反而通常会变好。
拟牛顿法与牛顿法的迭代过程一样，仅仅是各个Hessian矩阵的求解方法不一样。
在远离极小值点处，Hessian矩阵一般不能保证正定，使得目标函数值不降反升。而拟牛顿法可以使目标函数值沿下降方向走下去，并且到了最后，在极小值点附近，可使构造出来的矩阵与Hessian矩阵“很像”了，这样，拟牛顿法也会具有牛顿法的二阶收敛性。
对目标函数f(X)做二阶泰勒展开：
f(x)=f(xk+1)+(x−xk+1)T∇f(xk+1)+12(x−xk+1)TG(xk+1)(x−xk+1)
其中∇2f(xk) ⇔G(xk+1)
上式两边对x求导：
∇f(x)=∇f(xk+1)+G(xk+1)(x−xk+1)
当x=xk的时候，有
G−1(xk+1)[∇f(xk+1)−∇f(xk)]=xk+1−xk
这里我们用Hk表示在点xk处的Hessian矩阵的**逆矩阵**
Hk+1[∇fd(xk+1)−∇f(xk)]=xk+1−xk ......(5)
(5)式就是**拟牛顿方程**
下面提供两种拟牛顿法
### 非线性最小二乘(Least Square)拟合
问题：
给定(ti,yi),i=1,…n, 拟合一个函数y=f(t,x)，
其中x为待定的参数向量, f 对x非线性。
记误差ri(x)=yi−f(ti,x) r(x)=(r1(x),…,rn(x))T
优化模型：
minxR(x)=rT(x)r(x)
记r(x)的雅克比行列式为J(x)=(∂ri∂xj)n×m
∇R=2J(x)Tr(x)∇2R=2J(x)TJ(x)+2SS=∑ni=1ri(x)∇2ri(x)∇2ri(x)=(∂2ri∂xk∂xl)n×m
牛顿法要计算Hessian矩阵，其中S计算量大
若f 对x线性, 则化为线性最小二乘拟合, 此时
S=0
特定算法考虑如何忽略或近似矩阵S
#### Gauss-Newton算法(GN)：忽略矩阵S
∇R=2J(x)Tr(x)∇2R=2J(x)TJ(x)
牛顿方程：∇2f(xk)d=−∇f(xk)
f用R代替，下降方向dk满足
J(x)TJ(x)dk=−J(x)Tr(x)
GN算法收敛性特点:
S=∑ni=1ri(x) 收敛性依赖f对x的线性程度，即偏差r的大小（”小残差问题“）
#### LM算法: Levenbery (1944), Marquardt (1963)
G-N算法修正
J(x)TJ(x)dk=−J(x)Tr(x)
⇓
(J(x)TJ(x)+αkI)dk=−J(x)Tr(x)
其中αk>0为修正参数，防止JTJ出现病态。
LM算法特点
dk位于GN方向（αk很小）和负梯度方向(αk很大)之间
LM算法的另一种观点
min∥r(xk)+J(xk)(x−xk)∥2
s.t. ∥x−xk∥2≤hk (∗)置信域
其解为
(J(x)TJ(x)+αkI)dk=−J(x)Tr(x)
若GN方向满足约束(*)，αk=0,否则ak>0
置信域算法特点：
先确定最长步长hk，再确定方向dk
