机器学习第二课：无约束优化问题（局部极小值的几种解法）（梯度下降法与拟牛顿法）

上一篇讲了一些微积分的概括。（http://blog.csdn.net/dajiabudongdao/article/details/52397343） 并引出了梯度下降法与拟牛顿法。当计算机求函数的局部极小值时，我们选用这两种方法（或衍生方法）来解。求极值的问题很复杂，我么、们大致分为两种：无约束极值问题，与有约束极值问题。这一章主要讨论无约束极值问题。有约束的极值在后面会提到。（http://blog.csdn.net/dajiabudongdao/article/details/52462942）

一、无约束优化问题

无约束优化问题，就是我们中学常见的求函数极值的方法。一个函数求极值无需外部控制。一旦有外部条件控制了，我们可以用拉格朗日乘数法求解，那就是条件极值的事情了。但是问题出现了，极值不一定是最值。我们求极值容易，最值很难。这TM就令人尴尬了。

我们工业界的同志发明了“局部最小值”一词的含义就是想要说“我的极值运算就是某种意义上的最值运算”。

1.局部极小值的解法

局部极小值的通用解法就是迭代。

例如当前自变量X_k，对应F(X)_k 。那么下次迭代的X_k+1  = X_k   + a*d_k其中a代表挪动步长，d_k表示挪动方向。迭代是个趋近过程，若使结果收敛需要遵循一个基本要求：下一步的结果要比这一步好，即X_k+1比X_k更接近值极值。当求解极小值时，我们自然希望X_k+1的结果比X_k小。按照泰勒级数展开，我们可以写成

F(X_k+1 ) =F( X_k )  + d_kF^'(X_k)+...         (其中d = X_k+1-X_k)

并且X_k+1 < X_k。忽略后面的部分可得d_kF^'(X_k) <0。那么我们可知每次迭代，d的方向为负。。但d的值是什么呢？

d_k = -F^'(X_k)

因为一点的梯度的本质是一种方向。当两个方向相同时，相乘正结果最大，相反时，相乘负结果最大。

下面是我用python模拟的，深度下降的程序

# coding=utf-8import numpyimport theanoimport theano.tensor as T###多项式的梯度下降####### y# x# theta：容忍的误差，因为编程有误，这里的误差仅控制循环次数# alpha：步长# x:用来标识对谁求导的标量，x=T.dvector('x')如'x'标识对x求导# y:带上面标量的公式# dimen:变量的维数# y=x**2# 使用：steepest(x,y,3,theta=0.1,alpha=0.1)def steepest(x, y, dimen, theta=0.1, alpha=0.1):    theta = theta * theta    J, updates = theano.scan(lambda i, y, x: T.grad(y[i], x), sequences=T.arange(y.shape[0]), non_sequences=[y, x])    Gradientfunction = theano.function([x], J, updates=updates)    x_0 = numpy.array([0.] * dimen)    x_k = x_0    wuchav = numpy.diag(numpy.array(Gradientfunction(x_k)))    wucha = numpy.dot(wuchav, wuchav)    while (wucha > theta):        x_k_1 = x_k - alpha * numpy.diag(Gradientfunction(x_k))        wuchav = x_k_1 - x_k        wucha = numpy.dot(wuchav, wuchav)        x_k = x_k_1    return x_kif __name__ == '__main__':    x = T.dvector('x')    y = x ** 2 + 2 * x + 1    res = steepest(x, y, 1, theta=0.01, alpha=0.1)    print "极小值的参数位置为"    print res

运行结果如下：

2.进一步探讨局部极小值的解法

有个问题，上面的泰勒级数展开只到了一级。如果展开到二级就变成牛顿法。

F(X_k+1 ) =F( X_k )  + dF^'(X_k)+1/2 d^TH(X)d + ...      (其中d = Xk+1-Xk)

其中H为Hassion矩阵，求二阶导使用。因为当X趋近与极值时，F(X_k+1)与F(X_k)很接近。F对d(d = X_k+1-X_k)求导为0，得的

d_k = -H^-1(X_k)F^'(X_k)

_{后面的算法与上面一样。但说实话，矩阵求逆是件很难的事情。很多情况，还无法求逆，这就造成牛顿法不太常用。}

  _{所以研究者提出了很多使用方法来近似Hessian矩阵,，这些方法都称作准牛顿算法,，BFGS就是其中的一种，以其发明者Broyden, Fletcher, Goldfarb和Shanno命名。}

_{我们用BFGS方法直接构造出Hassion矩阵的逆。具体步骤很复杂这里仅 帖公式，步骤如下：}

S_k = H^-1(X_k),

r_k=F^'(X_k+1)-F^'(X_k)

δ_k = X_k+1-X_k

同样的我用python模拟，BFGS的程序如下：

# coding=utf-8import numpy as npimport theanoimport theano.tensor as T###多项式的BFGS####### theta：容忍的误差# alpha：步长# x:用来标识对谁求导的标量，x=T.dvector('x')如'x'标识对x求导# y:带上面标量的公式# dimen:变量的维数# y=x**2# 使用：BFGS(x,y,3,theta=0.1,alpha=0.1)def BFGS(x, y, dimen, theta=0.1, alpha=0.1):    theta = theta * theta    S_0 = np.eye(dimen)    S_k = S_0    J, updates = theano.scan(lambda i, y, x: T.grad(y[i], x), sequences=T.arange(y.shape[0]), non_sequences=[y, x])    Gradientfunction = theano.function([x], J, updates=updates)    x_0 = np.array([0.] * dimen)    x_k = np.mat(x_0)    wuchav = np.diag(np.array(Gradientfunction(x_0)))    wucha = np.dot(wuchav, wuchav)    while (wucha > theta):        x_k_1 = x_k - alpha * np.mat(np.diag(Gradientfunction(np.array(x_k)[0]))) * S_k.T        theta_k = x_k_1 - x_k        gamma_k = np.mat(np.diag(Gradientfunction(np.array(x_k_1)[0])) - np.diag(Gradientfunction(np.array(x_k)[0])))        S_k_1 = S_k + (1 + (gamma_k * S_k.T * gamma_k.T) / (theta_k * gamma_k.T)) * (        (theta_k.T * theta_k) / (theta_k * gamma_k.T))        S_k_1 = S_k_1 - (theta_k.T * gamma_k * S_k + S_k * gamma_k.T * theta_k) / (theta_k * gamma_k.T)        S_k = S_k_1        wuchav = np.array(theta_k)[0]        wucha = np.dot(wuchav, wuchav)        x_k = x_k_1    return x_kif __name__ == '__main__':    x = T.dvector('x')    y = x ** 2 + 2 * x + 1    res = BFGS(x, y, 1, theta=0.001, alpha=0.1)    print "极小值的参数位置为"    print res