线性搜索

精确线性搜索——黄金分割法

单峰函数

设f(x)是[a,b]上的一元函数，x_m是f(x)在[a,b]上的极小值点，且对任意的x1,x2∈[a,b]，x1<x2，有：

1.当x2<x_m时，f(x1)>f(x2)；

2.当x1>x_m时，f(x1)<f(x2)。

则f(x)是单峰函数。

通过计算区间内两个不同点的函数值，就能确定包含极小值点的子区间。

算法步骤

非精确线性搜索

精确搜索方法缺点:

1.精确一维搜索方法一般需要花费很大的工作量, 特别是当迭代点远离问题的解时, 精确的求解一个一维子问题通常不是有效的.

2.有些最优化方法,其收敛速度并不依赖于精确一维搜索过程.

不精确一维搜索方法：

总体希望收敛快，每一步不要求达到精确最小，速度快，虽然步数增加，则整个收敛达到快速.

 

Armijo准则

Armijo(1966)和Goldstein(1965)分别了提出了不精确一维搜索过程.

设

是一个区间， 区间为[0,a]。

  下降准则: 目标函数单调下降, 同时要求下降不是太小(否则下降序列的极限值可能不是极小值).

  必须避免所选择的太靠近区间的端点. 从而要求:

其中

满足上式要求的[a_k]构成区间J=[0,c],这排除了区间右端点附近的点.

Armijo算法步骤

（1）给定下降方向d_k ，

（2）若

则停止，返回a；否则，转（3）

（3）a=a*β，转（2）；

Wolfe 准则

在Goldstein 准则中, 步长因子的极小值可能被排除在可接受域之外。

为了克服这一缺点, Wolfe  给出了如下的条件

几何解释:

可接受点处的切线的斜率大于或等于初始斜率的倍(曲率条件)。

Wolfe准则如下

测试

在线性回归的梯度下降和正规方程组求解中，用梯度下降拟合线性回归模型，迭代次数达到上万次。

采用精确线性搜索，8次即能达到指定精度，比普通梯度下降更接近准确结果。

采用Wolfe非精确线性搜索，5次即能达到指定精度。

采用Armijo法，所需仍然是一万多次，可能是由于上文中提到的“步长因子的极小值可能被排除在可接受域之外”。

输出两种有效搜索方式得到的步长，可见步长在1e-6和1e-2两个数量级变化；另外两种方法得到的步长都在1e-6数量级。

# coding:utf-8
import numpy as np


def g618(x,y,gradient,theta,alpha=1,iters=50):  # 精确线性搜索-黄金分割法
    a = 0
    b=alpha
    l = 0.618
    u = b - l * (b - a)
    v = a + l * (b - a)

    def loss(alpha):
        return np.sum((np.dot(x, theta-alpha*gradient)- y)**2)

    for i in range(iters):
        if b-a < 1e-8:
            break
        if loss(u) == loss(v):
            a = u
            b = v
        elif loss(u) < loss(v):
            b = v
            v = u
            u = b - l * (b - a)
        else:
            a = u
            u = v
            v = a + l * (b - a)
    print (a+b)/2,i
    return (a+b)/2

def armijo(x,y,gradient,theta,alpha=0.1,iters=50):  # 线性搜索-armijo法
    def loss(alpha):
        return np.sum((np.dot(x, theta-alpha*gradient)- y)**2)

    def phi(a):
        return -np.sum(2 * np.dot(x, gradient) * (np.dot(x, theta- a * gradient) - y))

    sigma=0.4
    rho=0.8
    loss_0=loss(0)
    # print phi(0)   #  -41319209434.1
    # print (loss(1e-8)-loss_0)*1e8   #  -41251360122.9

    def check(alpha):
        return loss(alpha)>loss_0+sigma*alpha*phi(0)

    for i in range(iters):
        if check(alpha):
            alpha*=rho
        else:
            break
    return alpha


def wolfe(x, y, gradient, theta, iters=50):  # 线性搜索-wolfe法
    def loss(alpha):
        return np.sum((np.dot(x, theta - alpha * gradient) - y) ** 2)

    def phi(a):
        return -np.sum(2 * np.dot(x, gradient) * (np.dot(x, theta- a * gradient) - y))

    sigma1 = 0.1
    sigma2 = 0.7
    loss_0 = loss(0)
    a1=0
    phi_0=phi(0)
    alpha=0.1

    for i in range(iters):
        if loss(alpha)>loss_0+ sigma1 * alpha * phi_0:
            temp=a1+(alpha-a1)/2/(1+(loss_0-loss(alpha))/(alpha-a1)/phi_0)
            alpha=temp

        elif phi(alpha) < sigma2 * phi_0:
            temp = a1 + (alpha - a1) / 2 / (1 + (loss_0 - loss(alpha)) / (alpha - a1) / phi_0)
            a1 = alpha
            alpha = temp
            loss_0=loss(alpha)
            phi_0=phi(alpha)
        else:
            break

    print alpha,i
    return alpha

def dataN(length):  # 生成数据
    x = np.zeros(shape = (length,2))
    y = np.zeros(shape = length)
    for i in range(0,length):
      x[i][0] = 1
      x[i][1] = i
      y[i] = (i + 25) + np.random.uniform(0,1) *10
    return x,y


def alphA(x,y):  # 选取前20次迭代cost最小的alpha
    c=float("inf")
    for k in range(1,1000):
            a=1.0/k**3
            f=gD(x,y,20,a)[1][-1]
            if f>c:
                break
            c=f
            alpha=a
    return alpha


def gD(x,y,iters,alpha):  # 梯度下降
    theta=np.ones(2)
    cost=[]
    for i in range(iters):
        hypothesis = np.dot(x,theta)
        loss = hypothesis - y
        cost.append(np.sum(loss ** 2))
        gradient = np.dot(x.transpose(),loss)
        theta -=alpha * gradient
        if cost[-1]<epsilon:
            break
    return theta,cost,i


def sD(x,y,iters,fun_num=1):  # 梯度下降，线性搜索
    funs=[g618,armijo,wolfe]
    theta=np.ones(2)
    cost=[]
    for i in range(iters):
        hypothesis = np.dot(x,theta)
        loss = hypothesis - y
        cost.append(np.sum(loss ** 2))
        gradient = np.dot(x.transpose(),loss)
        alpha=funs[fun_num](x,y,gradient,theta)
        theta -=alpha * gradient
        if cost[-1]<epsilon:
            break
    return theta,i


def eQ(x,y):  # 正则方程组
   x=np.matrix(x)
   y=np.matrix(y).T
   a=np.dot(x.T,x).I
   b=np.dot(a,x.T)
   c=np.dot(b,y)
   return c

length=100
iters=50000
epsilon=1000

x,y=dataN(length)
theta1,_,its1=gD(x,y,iters,alphA(x,y))
print theta1,its1
print sD(x,y,iters,0)  # g618
print sD(x,y,iters,1)  # armijo
print sD(x,y,iters,2)  # wolfe
print eQ(x,y)

"""
[ 27.4752872    1.04342315] 16921
3.04503223176e-06 39
0.0131747301801 35
3.04503223176e-06 39
0.0294140419714 37
3.04503223176e-06 39
0.0131718000286 37
3.04503223176e-06 39
0.0294143276007 37
3.04503223176e-06 39
(array([ 29.52746209,   1.01248425]), 8)
(array([ 27.47554186,   1.04368584]), 14038)
3.0449186251e-06 1
0.0294131338763 1
3.0449186251e-06 1
0.0294131338764 1
3.0449186251e-06 1
0.0294131338747 1
(array([ 29.88521493,   0.99985976]), 5)
[[ 30.36492265]
 [  0.99985744]]
"""