梯度下降（学习率优化，以F(w)=w^4为例）和拟牛顿

'''
学习率α如何确定 
    使用固定学习率还是变化学习率?
    学习率设置多大比较好?
下降方向  
    处理梯度方向，其他方向是否可以？ 
    可行方向和梯度方向有何关系？

'''
#----------------------------------------------
#固定学习率的梯度下降:以h_w=w^4为例

def fix_xr():
    w0=1.5
    w1=1.5
    h_x=w0**4
    a=0.01
    score1=[]
    score=[]
    for i in range(1000):
        g_w=4*(w0**3)#h_x的导数
        w0-=g_w*a
        g_w1=2*w1
        w1-=g_w1*a
        score.append(w0)
        score1.append(w1)
    print('s:',score[199],score[999])
    print('s1:',score1[199],score1[999])
    fig=plt.figure()
    ax=fig.add_subplot(1,1,1)
    ax.plot(np.arange(1000),score,'bo')
    ax.plot(np.arange(1000),score1,'ro')
    ax.set_title('蓝：y=x^4 | 红：y=x^2')
    plt.legend()
    plt.show()
fix_xr()#可以看出固定学习率在h_w函数不同时收敛的速率差别很大，函数越复杂可能收敛越慢
#----------------------------------------------------------------------
'''
优化学习率
调整学习率：
    在斜率(方向导数)大的地方，使用小学习率
    在斜率(方向导数)小的地方，使用大学习率
对梯度下降的运行过程进行分析，

1.wk=a，沿着负梯度方向，移动到wk+1=b，有：
    b=a-α*▽F(a) ==> f(a)>=f(b)
2.从w0为出发点，每次沿着当前函数梯度反方向移动 一定距离αk，得到序列：
    w0,w1,...,wn
3.对应的各点函数值序列之间的关系为：
   f(w0)>=f(w1)>=...>=f(wn)
4.当n 达到一定值时，函数f(w)收敛到局部最小值

#转换视角
记当前点为wk，当前搜索方向为dk(如：负梯度方向)，因为学习率α是待考察的对象，因此，将下列函数f(wk+αdk)看做是关于α的函数h(α)。
    h(α)= f(wk+αdk),α>0
    当α=0时，h(0)=f(wk)
    导数▽h(α)=▽f(wk+αdk) T dk  #T即转置
   
学习率的计算标准
    因为梯度下降是寻找f(w)的最小值，那么， 在wk和dk给定的前提下，即寻找函数 f(wk+αdk)的最小值。即：
    ? 因为梯度下降是寻找f(x)的最小值，那么， 在xk和dk给定的前提下，即寻找函数 f(xk+αdk)的最小值。即：
        α=arg min(h(α))=arg min(f(xk+αdk)),α>0
    进一步，如果h(α)可导，局部最小值处的α满足：
        h'(α)=▽f(wk+αdk) T dk=0
将α=0带入上式，得h'(0)=▽f(wk) T dk 
下降方向dk可以选负梯度方向(或者选择与负梯度夹角小于90°的方向)
从而可得：h'(0)<0

##如果能够找到足够大的α，使得h'(α)>0 则必存在某α，使得h'(α*)=0
            α*即为要寻找的学习率
##然后可使用 
1.二分线性搜索(Bisection Line Search) 
    不断将区间[α1, α2]分成两半，选择端点异号的一侧，直到区间足够小或者找到当前最优学习率。
2.回溯线性搜索(Backing Line Search)
    基于Armijo准则计算搜素方向上的最大步长，其基本思想是沿着搜索方向移动一
个较大的步长估计值，然后以迭代形式不断缩减 步长，直到该步长使得函数值
f(wk+αdk)相对与当前函数值f(wk)的减小程度大于预设的期望值(即满足Armijo准则
)为止。
        f(wk+αdk)<=f(wk)+c1*α▽f(wk) T dk ,c1=(0,1)
        
##回溯与二分线性搜索的异同        
    二分线性搜索的目标是求得满足h'(α)≈0的最优步长近似值，而回溯线性搜索放松了对步长的约束，只要步长能使函数值有足够大 的变化即可。
    二分线性搜索可以减少下降次数，但在计算最优步长上花费了不少代价；回溯线性搜索找到一个差不多的步长即可
'''
import sys
sys.setrecursionlimit(100000)
class LearningRate(object):
    def __init__(self):
        self.at=0
        self.a0=0
        self.w0=1.5
        self.i=0
        self.score=[]
    def F(self,w):
        return w**4
    def f(self,w):
        return 4*(w**3)
    def ft(self,w,a,dk):
        return 4*((w+a*dk)**3)*dk
    def search_max_a(self,a1,w0):

        dk = -self.f(w0)
        max_a=self.ft(w0,a1,dk)
        while max_a<0:
            a1*=2
            a1+=1
            max_a = self.ft(w0, a1, dk)

        return a1

    def BLS(self,a0,a1,w):
        dk=-self.f(w)
        at=(a0+a1)/2
        optimal = self.ft(w, at, dk)
        if  abs(optimal)<=0.001:
            self.at=at
            return
        if optimal>0.001 :
            a1=at
            self.BLS(a0,a1,w)
        else :
            a0=at
            self.BLS(a0, a1, w)
    def dynamic_lr(self,w,at):
        t=self.f(w)
        w-=t * at
        return w
    def curve_plot(self,score):
        print('s:',score[10])
        plt.figure()
        plt.plot(np.arange(self.i+1),score,'bo')
        plt.show()
    def start(self,w0,n):
        self.score.append(w0)
        if self.i ==n:
            self.curve_plot(self.score)
            return
        else:
            a1 = self.search_max_a(self.a0,w0)
            print(a1)
            self.BLS(self.a0, a1,w0)
            a2=self.at
            w0=self.dynamic_lr(w0,a2)
            self.i+=1
            self.start(w0,n)
lr=LearningRate()
lr.start(1.5,20)#可以得到大概在第十次迭代就已经是-3.1215834180041617e-05
#但是正如前面所说，在搜索最优步长时花费了很大的计算代价
#-------------------------------------------------------
#回溯线性搜索
class GetA_Armijo(object):
    def __init__(self):
        self.c1=0.3
        self.score=[]
        self.i=0

    def F(self,w):
        return w**4
    def f(self,w):
        return 4*(w**3)
    def ft(self,w,a,dk):
        return 4*((w+a*dk)**3)*dk
    def dynamic_lr(self,w,at):
        t=self.f(w)
        w-=t * at
        return w
    def BLS(self,w,a,d):
        now=self.F(w)
        next=self.F(w-a*d)
        print(now,next)
        c=30
        while(next<now):
            a*=2
            next=self.F(w-a*d)
            c-=1
            if c==0:
                break
        c=50
        while(next>now-self.c1*a*d*d):
            a /=2
            next=self.F(w-a*d)
            c-=1
            if c==0:
                break
        return a
    def curve_plot(self,score):
        print('s:',score[12],score[100])
        plt.figure()
        plt.plot(np.arange(self.i+1),score,'bo')
        plt.show()
    def start(self,w0,n):
        self.score.append(w0)
        if self.i ==n:
            self.curve_plot(self.score)
            return
        else:
            at=self.BLS(w0,0.02,self.f(w0,))
            w0=self.dynamic_lr(w0,at)
            self.i+=1
            self.start(w0,n)
tt=GetA_Armijo()
tt.start(1.5,120)
#第12次迭代: 0.0003698107138697752
#第100次迭代：2.554384621316504e-06

#----------------------------------------------------------

'''
#采用多项式插值法拟合简单函数，
    然后根据该简单函数估计函数的极值点，这样选择合适步长的效率会高很多
已知数据为wk处的函数值f(wk)及其导数f'(wk),再加上第一次尝试的随机步长
α0.
   如果α满足条件，显然算法退出；
   若α不满足条件，则根据上述信息可以构造一个二次近似函数 
                h(α0)-h'(0)*α0 -h(0)
        h(α)=  ----------------------- α^2 +h'(0)α+h(0)
                        α0^2
    由来：
        设二次函数f(x)=bx^2+cx+d,已知f(0),f'(0),f(a)
        f(0)=d
        f'(0)=c    =======>f(x)
        f(a)=b*a^2+c*a+d 
                                              h'(0)*α0^2
    显然，导数为0的最优值[-b/2a]==>α1= ---------------------------
                                         2[h'(0)*α0 + h(0)-h(α0)]
   若α1满足 Armijo准则，则输出该学习率；否则，继续迭代    
'''
class Quadratic_Interpolation(object):
    def __init__(self):
        self.c1=0.3
        self.score=[]
        self.i=0

    def F(self,w):
        return w**4
    def f(self,w):
        return 4*(w**3)
    def ft(self,w,a,dk):
        return 4*((w+a*dk)**3)*dk
    def dynamic_lr(self,w,at):
        t=self.f(w)
        w-=t * at
        return w
    def BLS(self,w,a,d):#d为F(w)的导数
        now=self.F(w)
        next=self.F(w-a*d)
        print(now,next)
        c=30
        while(next<now):
            a*=2
            next=self.F(w-a*d)
            c-=1
            if c==0:
                break
        c=50
        while(next>now-self.c1*a*d*d):
            b=d*a*a/(now+d*a-next)#二次插值
            b /=2#二次插值
            if b<0:#如果b<0则使用回溯法求解
                a/=2
            else:#否则继续使用二次插值
                a=b
            next=self.F(w-a*d)
            c-=1
            if c==0:
                break
        return a
    def curve_plot(self,score):
        print('s:',score[12],score[100])
        plt.figure()
        plt.plot(np.arange(self.i+1),score,'bo')
        plt.show()
    def start(self,w0,n):
        self.score.append(w0)
        if self.i ==n:
            self.curve_plot(self.score)
            return
        else:
            at=self.BLS(w0,0.02,self.f(w0,))
            w0=self.dynamic_lr(w0,at)
            self.i+=1
            self.start(w0,n)
qi=Quadratic_Interpolation()
qi.start(1.5,120)

'''
通过使用线性搜索的方式，能够比较好的解决学习率问题
一般的说，回溯线性搜索和二次插值线性搜索能够 基本满足实践中的需要
问题：
    可否在搜索过程中，随着信息的增多，使用三次或者更高次的函数曲线，从而
得到更快的学习率收敛速度？
    为避免高次产生的震荡，可否使用三次Hermite(特征值均为实数的方阵)多项式，在端点保证函数值和一阶导都相等，从而构造更光顺的简单低次函数？
'''

'''
#搜索方向
因为函数二阶导数反应了函数的凸凹性；二阶导越大，一阶导的变化越大。在搜索中， 可否用二阶导做些“修正”？如：二者相除？
x(k+1) =xk - f'(xk)/f''(xk)
由来：
   若f(x)二阶导连续，将f(x)在xk处Taylor展开：
        h(x)=f(xk)+f'(xk)(x-xk)+1/2 * f''(xk)(x-xk)^2 + R(x)
        h'(x)≈f'(xk)+f''(xk)(x-xk)=0
        => x =xk - f'(xk)/f''(x)  
    上式即，牛顿法
    该方法可以直接推广到多维：用方向导数代替一阶导，用Hessian矩阵代替二阶导

经典牛顿法虽然具有二次收敛性，但是要求初始点 需要尽量靠近极小点，否则有可能不收敛。
? 计算过程中需要计算目标函数的二阶偏导数，难度较大。
? 目标函数的Hessian矩阵无法保持正定，会导致算法产生的方向不能保证是f在xk处的下降方向，从而令牛顿法失效；
? 如果Hessian矩阵奇异[方阵行列式=0，有无穷解]，牛顿方向可能根本是不存在的。

BFGS / LBFGS ------拟牛顿
Broyden – Fletcher – Goldfarb - Shanno
    求Hessian矩阵的逆影响算法效率，同时，搜索方向只要和负梯度的夹角小于90°即可， 因此，可以用近似矩阵代替Hessian矩阵，只要满足该矩阵正定、容易求逆，或者可以通过若干步递推公式计算得到
    矩阵迭代公式：
    
               yk(yk)T        (Bksk)(Bksk)T
    Bk+1=Bk + ------------ -  --------------
               (yk)T sk         (sk)T Bksk
        sk=xk - x(k-1)
        yk=gk - g(k-1) #k-1只表示迭代次序,g为▽f(xk)
    B0=I 单位阵
L-BFGS    
    BFGS需要存储n×n的矩阵Hk用于近似 Hessian矩阵的逆矩阵；
    而L-BFGS仅需要存 储最近m(m一般小于10，m=20足够)对n维的更新数据(x,f'(xi))即可。L-BFGS的空间复杂度O(mn)，若将m看做常数，则为线性，特别适用于变量非常多的优化问题。
    
'''

 Hessian非正定时情况