基于二阶矩阵的优化问题（一）线搜索策略
非精确线搜索更新$X_{k+1}$
线搜索策略
最速下降+牛顿法
启发
修正牛顿法
特征值修正（eigenvalue modification）
谱分解
单位增量法
拟牛顿法

非精确线搜索更新$X_{k+1}$
优化算法的问题在于如何从$X_k$更新到$X_{k+1}$、确定步长公式详见 基于二阶矩阵的优化问题（二）和判定迭代终止的条件（详见 于二阶矩阵的优化问题（三））是我们需要在不精确搜索中解决的问题。

我们有两种方法来解决这个问题：

1.线搜索策略

2.信任域策略
线搜索策略
我们选择一个下降方向来使得lost function慢慢迭代至最小值。

$X_{k+1}=X_k+α_kp_k$

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

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

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

$p^N_K=-(\nabla^2f_k)^{-1}\nabla f_k$;

这边的$p^N_K$就是这里的局部二次型的顶点，在牛顿法中，我们可以一步就走到这个点。

ok，这个问题的理论情况是这样的，但在实际问题中，我们不可能去求一个hessian矩阵的逆矩阵，所以，我们把$1/(\nabla^2f_k)$移到等式的左边，变成：

$p^N_K*(\nabla^2f_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_kp_k=-\nabla f_k$

        $x_{k+1}=x_k+α_kp_k$（其中$α_k$在wolfe条件时可以取1）

end

好了，现在我们需要解决四件事，就可以实现修正牛顿法的运用：

1.判定$\nabla^2 f_k$是否正定；

2.解方程组$B_kp_k=-\nabla f_k$；

3.如何构建$E_k$；

4.如何使得$\vert\vert E_k\vert\vert$最小；
特征值修正（eigenvalue modification）
谱分解
首先我们先对对称矩阵$B_k$做正交分解（householder等），得到正交阵Q，和对角阵V，然后就可以得到$B_k$的所有特征值，我们知道，正定矩阵所有的特征值都是充分大的正数，即$λ_{min}>δ>>eps>0$，当出现负的特征值时，我们添加一个修正项$E_k$使得特征值大于等于δ。此时的$\vert\vert E_k\vert\vert$的Frobenius范数最小。

即用Matlab中的eig()函数求出特征值后，直接对特征值修改即可。
单位增量法
谱分解的问题在于，其真正改变的是特征值，所以对$B_k$的改动是比较大的，当我们不希望$B_k$发生很大的变动时，我们就使用单位增量法。单位增量法直接对$B_k$添加一个修正项$E_kI$，使得其正定。

首先，我们对$B_k$进行cholesky分解：

$B_k=LDL^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)*\frac{x_{k-1}-x_k}{f(x_{k-1})-f(x_k)}$

这在本质上是计算函数的数值微分，但这种方法在接近收敛时会出现数值不稳定的情况，当接近收敛时，舍入误差占比变大，数值会发生很大的起伏变化。

由此我们可以推导出拟牛顿法的公式

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

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

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

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

$B_{k+1}*s_k=y_k$

这里的$B_{k+1}$是$n^2$阶的矩阵，而我们只有n个方程，下面有SR1、BFGS等方法来求解$B_{k+1}$。

详见   拟牛顿法的下降方向计算（一）

此时，局部二次型的$p_K^N$方程变为

$p^N_K=-(B_k)^{-1}\nabla f_k$;

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

$p_k=-H_k*\nabla f_k$

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

1   二阶范数 
||A||2=tr(AAT)


2 矩阵求导基本公式（等式两边同时去除tr，公式不变）： 
tr(A)=tr(AT)

∂tr(AX)∂X=AT
∂tr(AX)∂XT=A
∂tr(AXB)∂X=(BA)T=ATBT


3   则简单公式均可以推导出来： 
∂tr(XTAX)∂X=AX+(XTA)T=(AT+A)X

∂tr(XTAX)∂XT=(AX)T+XTA=XT(AT+A)

∂tr(ABATC)∂A=(BATC)T+CAB=CTABT+CAB

学习阶段：高中数学。
前置知识：等比数列。
1. 斐波那契数列的定义
斐波那契数列大部分人都耳熟能详，就是 
 ，每一项都是前两项之和。我们用数列的递推公式严格定义一下：
2. 从零推导斐波那契数列的通项公式
相信也有不少的人见过斐波那契数列的通项公式，感觉还是比较复杂的，有什么 
 啊， 
 次方等等的复杂结构。成熟的推导通项公式的方法有
特征方程法、生成函数/母函数法、矩阵对角化法等等，但这些方法“集成度”过高，需要很多其他知识的铺垫，理解起来不是很清晰。现在我们忘记这些方法，从零开始推导一下它的通项公式，便于大家记（zhuang）忆（bi）。
等差数列 
 和等比数列 
 我们已经研究透彻了，我们一般的思路就是把待求数列化归为这两种数列。
2.1 二阶递推化为一阶递推
首先， 
 这个鬼二阶递推公式，似乎没办法化归啊？没事，我们暴力假设它能变形为 
 这种等比数列的形式（有人问我这个暴力假设怎么来的。呃，这要一点凑的数学直觉，就像凑不等式和凑微分。把递推式子凑出等比数列的形式，还要每项都能对应上，差不多只能设出这种形式），然后看看能不能把常数 
 与 
 求出来。对应系数，有：
成功解出 
 ，说明 
 和 
 是可以解出来的嘛！
然后利用等比数列的性质，有
成功地把二阶的递推式化成了一阶递推式！
2.2 一阶递推化为通项公式
现在故技重施，还是暴力假设把它化为等比数列的形式，巧妙地设计待定的系数：
因为 
设 
得 
哇哦！这下子 
 的通项已经被彻底地解出来了，剩下的只是繁琐的计算：
现在要把 
 的值全部带进去。 
 的值有两个，随便选一个就行了（结果是一样的），我选 
 ，则 
 ，立刻得到：
其实，这个过程正是特征方程法的由来。待定系数法是如此强大！

我有一个损失值/函数,我想计算所有关于张量f(大小为n)的二阶导数.我设法使用了两次tf.gradients,但是当第二次应用它时,它会对第一个输入中的导数求和(请参阅我的代码中的second_derivatives).
我还设法检索Hessian矩阵,但我想只计算它的对角线以避免额外的计算.
import tensorflow as tf
import numpy as np
f = tf.Variable(np.array([[1., 2., 0]]).T)
loss = tf.reduce_prod(f ** 2 - 3 * f + 1)
first_derivatives = tf.gradients(loss, f)[0]
second_derivatives = tf.gradients(first_derivatives, f)[0]
hessian = [tf.gradients(first_derivatives[i,0], f)[0][:,0] for i in range(3)]
model = tf.initialize_all_variables()
with tf.Session() as sess:
sess.run(model)
print "\nloss\n", sess.run(loss)
print "\nloss'\n", sess.run(first_derivatives)
print "\nloss''\n", sess.run(second_derivatives)
hessian_value = np.array(map(list, sess.run(hessian)))
print "\nHessian\n", hessian_value
我的想法是tf.gradients(first_derivatives,f [0,0])[0]可以检索例如关于f_0的二阶导数,但似乎张量流不允许从张量的张量中导出.