基于二阶矩阵的优化问题（一）线搜索策略（附matlab代码）

非精确线搜索更新$X_{k+1}$

优化算法的问题在于如何从$X_k$更新到$X_{k+1}$、确定步长公式详见 基于二阶矩阵的优化问题（二）和判定迭代终止的条件（详见 于二阶矩阵的优化问题（三））是我们需要在不精确搜索中解决的问题。
我们有两种方法来解决这个问题：
1.线搜索策略
2.信任域策略

线搜索策略

我们选择一个下降方向来使得lost function慢慢迭代至最小值。

$X_{k+1}=X_k+{\alpha }_k{p}_k$

这是$X_{k+1}$的更新公式，可以看到，我们需要知道其${\alpha }_k$和$p_k$的值才能进行迭代（${\alpha }_k$代表了步长、$p_k$代表下降方向），第一种方法就是$p_k$走一个最速下降方向（负梯度），而${\alpha }_k$走一个极其猥琐的距离，但这种方法在函数比较诡异时效果较差（如rosenbrock函数），但在某些复杂问题时，我们还是使用最速下降来解决问题，详见 大规模的优化问题（一）。

最速下降+牛顿法

让我们来列一下牛顿方向的式子：

$f(x_k+p)\approx f_k+p^T\nabla f_k+\frac{1}{2}{p}^T{\nabla }^2{f}_{k}p$

在这里，我们要假设这里的二次型（${\nabla }^2{f}_k$）是正定的，这是函数的局部二次型构型就是碗状的，我们可以求出这里的pk（下降方向）：

$p_K^N=-({\nabla }^2{f}_k{)}^{-1}\nabla f_k$;

这边的$p_K^N$就是这里的局部二次型的顶点，在牛顿法中，我们可以一步就走到这个点。
ok，这个问题的理论情况是这样的，但在实际问题中，我们不可能去求一个hessian矩阵的逆矩阵，所以，我们把$1/({\nabla }^2{f}_k)$移到等式的左边，变成：

$p_K^N\ast ({\nabla }^2{f}_k)=-\nabla f_k$，

这个等式实际上就等于Ax=b的形式，这边我们只需要求解一个线性方程组即可。
以上为牛顿法的基础理论，下面看一下牛顿法的实现步骤：

可以看到，牛顿法的意义在于足够靠近真解时的快速收敛，我们先走几部最速下降法，使得足够靠近阈值，后再启动牛顿法，会有较好的效果,即保证了算法的效率，有保证了算法的鲁棒性。

下面是最速下降法+牛顿法matlab实现的代码：

function [x1] = min_hybrid(func, dfunc, d2func ,x, MAX_IT ,eps, threshold)
    err = 10.0;
    steps = 0;
    alpha = 1.0;
    //最速下降法
    while (err > threshold)//此处阈值需要调参
        f = func(x);
        g = dfunc(x);
        p = -g;
        step=step+1;
        alpha = backtracking(@func, x, f, g, p, 1e-4, 0.8, alpha * 2.0);
        x = x + alpha * p;
        if (steps > MAX_IT)
            break;
        end
        error = norm(g);
    end
    //牛顿法
    while (error > eps)
        f = func(x);
        g = dfunc(x);
        G = d2func(x);
        p = -G\g;
        x = x + p;
        step=step+1;
        if (steps > MAX_IT)
            break;
        end
        error = norm(g);
    end
	x1=x;
end

function alpha = backtracking(func, x, f, df, p, c, rho, alpha0)
    alpha = alpha0;
    while (func(x(1) + alpha * p(1), x(2) + alpha * p(2)) > f + c * alpha * df' * p)
        alpha = alpha * rho;
    end
end

其中func、gfunc和hess需要自行解析计算，func的接口只有x，gfunc和hess没有输入参数。（这边要注意，使用wolfe条件时，alpha可以恒取1）

启发

牛顿法实际上是在改进真解，用少量次数的迭代即可实现一阶方法上千次的迭代效果，首先系数矩阵必须要正定，否则牛顿法不能保证收敛。

修正牛顿法

牛顿法最大的问题在于，系数矩阵必须要正定，那么如果系数矩阵不正定，则将$B_k$改成正定（modification），下面是修改系数矩阵的算法框架：

初值$x_0$
for k=0,1,2,…
        令 $B_k={\nabla }^2{f}_k+E_k$,其中$E_k$为修正项，使得$B_k$充分正定。
        求解 $B_k{p}_k=-\nabla f_k$
        $x_{k+1}=x_k+{\alpha }_k{p}_k$（其中${\alpha }_k$在wolfe条件时可以取1）
end

好了，现在我们需要解决四件事，就可以实现修正牛顿法的运用：
1.判定${\nabla }^2{f}_k$是否正定；
2.解方程组$B_k{p}_k=-\nabla f_k$；
3.如何构建$E_k$；
4.如何使得$||E_k||$最小；

特征值修正（eigenvalue modification）

谱分解

首先我们先对对称矩阵$B_k$做正交分解（householder等），得到正交阵Q，和对角阵V，然后就可以得到$B_k$的所有特征值，我们知道，正定矩阵所有的特征值都是充分大的正数，即${\lambda }_{min}>\delta >>eps>0$，当出现负的特征值时，我们添加一个修正项$E_k$使得特征值大于等于δ。此时的$||E_k||$的Frobenius范数最小。
即用Matlab中的eig()函数求出特征值后，直接对特征值修改即可。

单位增量法

谱分解的问题在于，其真正改变的是特征值，所以对$B_k$的改动是比较大的，当我们不希望$B_k$发生很大的变动时，我们就使用单位增量法。单位增量法直接对$B_k$添加一个修正项$E_{k}I$，使得其正定。
首先，我们对$B_k$进行cholesky分解：

$B_k=LD{L}^T$

L是对角线均为1的下三角阵，D为对角阵，若$B_k$不定，D的对角元$d_{ii}$会过大。
下面是单位增量法的Matlab代码：

function[L,D]=modification(A,delta,beta)
	if(norm(A-A')>eps
		return;
	end
	n=size(A,1);
	d=zero(n,1);
	L=zero(n);
	C=zero(n);
	theta=zeros(n,1);
	for j = 1:n
		C(j,j) = A(j,j)-sum(d(1:j-1)'.*L(j,1:j-1).^2);
		for i = j+1:n
			C(i,j) = A(i,j)-sum(d(1:j-1)'.*L(i,1:j-1).*L(j,1:j-1));
			absjj=abs(C(i,j));
			if theta < abs(C(i,j))
				theta = abs(C(i,j));
			end
		end
		d(j) = max([abs(C(j,j)), (theta/beta)^2, \delta]);
		for i = j+1:n
			L(i,j) = C(i,j)/d(j);
		end
		L(j,j) = 1.0;
	end
	D=diag(d);
end

单位增量法得到的$B_k$和$B_k+E_k$在形式上很接近，但其特征值完全不同，所以两种方法各有优势，可根据情况自行选择。

拟牛顿法

如果你的数据维数过大或者无法承受计算二阶矩阵的消耗，这边也可以使用拟牛顿法来计算下降方向，当然拟牛顿法也有其缺陷，下面我们来分析一下。
首先我们介绍一下割线法，用弦的斜率近似代替目标函数的切线斜率，下面给出割线法的公式

$x_{k+1}=x_k-f(x_k)\ast \frac{x_{k-1}-x_k}{f(x_{k-1})-f(x_k)}$

这在本质上是计算函数的数值微分，但这种方法在接近收敛时会出现数值不稳定的情况，当接近收敛时，舍入误差占比变大，数值会发生很大的起伏变化。
由此我们可以推导出拟牛顿法的公式

$x_{k+1}=x_k-\nabla f(x_k)\ast \frac{x_{k-1}-x_k}{\nabla f(x_{k-1})-\nabla f(x_k)}$

这里不在需要计算其hessian矩阵，由于hessian矩阵是$n^2$级别的，拟牛顿法对于大规模问题来说是非常节约资源的方法。

于是我们知道了

${\nabla }^{2}f(x_k)\ast (x_{k+1}-x_k)\approx \nabla f_{x+1}-\nabla f_k$

现在的情况与牛顿发不同，$\nabla f(x_k)$并不知道，这是一个不定方程组:

$B_{k+1}\ast s_k=y_k$

这里的$B_{k+1}$是$n^2$阶的矩阵，而我们只有n个方程，下面有SR1、BFGS等方法来求解$B_{k+1}$。
详见 拟牛顿法的下降方向计算（一）
此时，局部二次型的$p_K^N$方程变为

$p_K^N=-(B_k{)}^{-1}\nabla f_k$;

这里还有一条路，就是去构建$(B_k{)}^{-1}$,这时我们不需要再去求解线性方程组（这里设$H_k=(B_k{)}^{-1}$），这里$B_k$=I时，就是最速下降。

$p_k=-H_k\ast \nabla f_k$

这里有拟BFGS等方法来求解$H_k$,详见 拟牛顿法的下降方向计算（二）

github代码地址