Optimization of Machine Learning

机器学习就是须要找到模型的鞍点,也就是最优势。由于模型不少时候并非彻底的凸函数,因此若是没有好的优化方法可能会跑不到极值点,或者是局部极值,甚至是偏离。因此选择一个良好的优化方法是相当重要的。首先是比较常规的优化方法:梯度降低。如下介绍的这些算法都不是用于当个算法,能够试用于能可微的全部算法。 ###Gradient Descent 常见会用在logistics regression或者是linear regression里面。好比logistics regression,这个模型直接等于0是求不出解析解的,全部只能用一些迭代方法求解。当样本的label是\{1, 0\} 每个样本正确几率:P(y|x;\theta) = (h(x_i)^{y_i}(1-h(x_i))^{1-y_i}) 最大似然函数log化简:l(\theta) = \sum_{i = 0}^{m}(y_ilogh(x_i) + (1-y_i)logh(x_i)) target function就是如上的l(\theta),要作的就是优化上面的损失函数。 \frac{d(l(\theta))}{d\theta} = -\frac{1}{m}\sum_{i=0}^{m}\{y_i\frac{1}{h(x_i)}\frac{d}{d\theta}h(x_i)-(1-y_i)\frac{1}{1-h(x_i)}\frac{d}{d\theta}h(x_i)\} = -\frac{1}{m}\sum_{i=0}^{m}\{y_i\frac{1}{g(\theta^Tx_i)}-(1-y_i)\frac{1}{1-g(\theta^Tx_i)}\}\frac{d}{d\theta}g(\theta^Tx_i) = -\frac{1}{m}\sum_{i=0}^{m}\{y_i\frac{1}{g(\theta^Tx_i)}-(1-y_i)\frac{1}{1-g(\theta^Tx_i)}\}g(\theta^Tx_i)(1-g(\theta^Tx_i))\frac{d}{d\theta}\theta^Tx_i =-\frac{1}{m}\sum_{i=0}^{m}\{y_i(1-g(\theta^Tx_i)) - (1-y_i)g(\theta^Tx_i)\}x_i^j = \frac{1}{m}\sum_{i=0}^{m}(h(x_i)-y_i)x_i^j 这个就是当前点的梯度,而函数的更新方向就是向着梯度的负方向进行的,因此最后是要减:git

\theta_j = \theta_j - \alpha\frac{1}{m}\sum_{i=0}^{m}(h(x_i)-y_i)x_i^j

最优势很明显是在左边,而待优化点是 P_0点,而在 P_0点的梯度是上升的,因此负方向才是正确的,最后获得的梯度就是要减去了。这个就是gradient descent。对于步长 \alpha,能够选择 \frac{\alpha}{
\sqrt{t}},t就是循环次数,随着迭代次数的增多,愈来愈靠近极值点,再相同的迭代步长是有可能越过极值点的。事实上步长的选择也是一个很重要的方法,后面会使用到。梯度降低是用给一个平面去拟合,也就是一次曲线去拟合,由于梯度降低是用泰勒一阶展开的,因此是一次曲线。这样形成的问题就是,gradient descent方法的视野很小,由于它只有一阶拟合,因此当前找到那个方向并不必定是全局最优的方向,多是有点偏离,因此常常出现之字形的走势,迭代的速度也不快。它只会看当前一阶导数哪一个最适合,而不会去管二次或者三次的。 **梯度降低有三种:全量梯度降低,全部是数据一块儿作迭代更新。数据量很大的话速度回很慢。小批量梯度降低,一部分一部分数据的就行迭代数据。随机梯度降低SGD,随机选取几个进行迭代,可能迭代的方向会有误差,可是随着时间流逝大方向仍是同样的。**代码实现前面的logistics regression中已经有了,再也不重复。 给予这种不足,给出了牛顿法的改进。 ###NewTon Method 既然梯度降低作的是一阶拟合,那么试一下二阶拟合。泰勒展开二阶: f(x) = f(x_0) + f'(x_0)(x - x_0) + \frac{f''(x_0)}{2!}(x-x_0)^2两边对 x求导数: f'(x) = f'(x_0) + f''(x_0)(x-x_0) = 0 x - x_0 = -\frac{f'(x_0)}{f''(x_0)} x = x_0 - \frac{f'(x_0)}{f''(x_0)}代换一下: x_{n+1} = x_n - \frac{f'(x_n)}{f''(x_n)} 牛顿法有点像坐标上升的方法,先迭代一个坐标,再迭代另外一个,和SMO相似。能够看到牛顿法这个时候就是二次曲线的拟合了:
左边的就是成功拟合的,右边的那个就是方向错误,拟合没有成功的例子。因此牛顿法不必定会按照正确的方向拟合。上面牛顿法的式子是对于单变量的,若是是对变量,那么下面的二阶导要用到Hession矩阵。因此对于多变量的牛顿法:
d_k = -H_k^{-1}g_k称为是牛顿方向,刚刚提到,迭代方向要和一阶导数方向相反,也就是 g_k^Td_k < 0,这个想法和刚刚提出的沿着梯度方向的反方向降低的思路是同样的。分解一下这个公式: g_k^TH_k^{-1}g_k > 0,因此Hession矩阵要是正定矩阵并且非奇异,而正定矩阵这个条件是很强的,不是在特定条件下达不到,因此牛顿法虽然降低的速度能够很快,可是方向不必定是正确的,因此牛顿法要使用的话必定要要求是靠近极值点的状况下使用,不然容易偏离。从本质上去看,牛顿法是二阶收敛,梯度降低是一阶收敛,因此牛顿法就更快。若是更通俗地说的话,好比你想找一条最短的路径走到一个盆地的最底部,梯度降低法每次只从你当前所处位置选一个坡度最大的方向走一步,牛顿法在选择方向时,不只会考虑坡度是否够大,还会考虑你走了一步以后,坡度是否会变得更大。因此,能够说牛顿法比梯度降低法看得更远一点,能更快地走到最底部。(牛顿法目光更加长远,因此少走弯路;相对而言,梯度降低法只考虑了局部的最优,没有全局思想。)这一段简单来讲就是梯度降低仅仅考虑了梯度,而牛顿法不只考虑了梯度,并且还考虑了梯度的变化,也就是梯度的梯度。 #####代码实现

def sigmoid(
        x1, x2, theta1, theta2, theta3):
    z = (theta1*x1+ theta2*x2+theta3).astype("float_")
    return 1.0 / (1.0 + np.exp(-z))

'''cross entropy loss function '''
def cost(
        x1, x2, y, theta_1, theta_2, theta_3):
    sigmoid_probs = sigmoid(x1,x2, theta_1, theta_2,theta_3)
    return -np.mean(y * np.log(sigmoid_probs)
                  + (1 - y) * np.log(1 - sigmoid_probs))
复制代码

sigmoid函数和cost函数,使用牛顿法代替梯度降低作逻辑回归。github

def gradient(
        x_1, x_2, y, theta_1, theta_2, theta_3):
    sigmoid_probs = sigmoid(x_1,x_2, theta_1, theta_2,theta_3)
    return np.array([np.sum((y - sigmoid_probs) * x_1),
                     np.sum((y - sigmoid_probs) * x_2),
                     np.sum(y - sigmoid_probs)])

复制代码

一阶导数,数据就是二维的,再加上一个偏置值就是三个导数了。算法

'''calculate hassion matrix'''
def Hassion(
        x1, x2, y, theta1, theta2, theta3):
    sigmoid_probs = sigmoid(x1,x2, theta1,theta2,theta3)
    d1 = np.sum((sigmoid_probs * (1 - sigmoid_probs)) * x1 * x1)
    d2 = np.sum((sigmoid_probs * (1 - sigmoid_probs)) * x1 * x2)
    d3 = np.sum((sigmoid_probs * (1 - sigmoid_probs)) * x1 )
    d4 = np.sum((sigmoid_probs * (1 - sigmoid_probs)) * x2 * x1)
    d5 = np.sum((sigmoid_probs * (1 - sigmoid_probs)) * x2 * x2)
    d6 = np.sum((sigmoid_probs * (1 - sigmoid_probs)) * x2 )
    d7 = np.sum((sigmoid_probs * (1 - sigmoid_probs)) * x1 )
    d8 = np.sum((sigmoid_probs * (1 - sigmoid_probs)) * x2 )
    d9 = np.sum((sigmoid_probs * (1 - sigmoid_probs)))
    H = np.array([[d1, d2,d3],[d4, d5,d6],[d7,d8,d9]])
    return H
复制代码

方法比较笨,直接一个一个迭代。bash

'''update parameter by newton method'''
def newton_method(x1, x2, y):
    #initialize parameter
    theta1 = 0.001
    theta2 = -0.4
    theta3 = 0.6
    delta_loss = np.Infinity
    loss = cost(x1, x2, y, theta1, theta2, theta3)
    kxi = 0.0000001
    max_iteration = 10000
    i = 0
    print('while i = ', i, ' Loss = ', loss, 'delta_loss = ', delta_loss)
    while np.abs(delta_loss) > kxi and i < max_iteration:
        i += 1
        g = gradient(x1, x2, y, theta1, theta2, theta3)
        hassion = Hassion(x1, x2, y, theta1, theta2, theta3)
        H_inverse = np.linalg.inv(hassion)
        delta = H_inverse @ g.T
        delta_theta_1 = delta[0]
        delta_theta_2 = delta[1]
        delta_theta_3 = delta[2]
        theta1 += delta_theta_1
        theta2 += delta_theta_2
        theta3 += delta_theta_3
        loss_new = cost(x1, x2, y, theta1, theta2, theta3)
        delta_loss = loss - loss_new
        loss = loss_new
        print('while i = ', i, ' Loss = ', loss, 'delta_loss = ', delta_loss)
    return np.array([theta1, theta2, theta3])
复制代码

运行函数。app

def get_param_target(self):
        np.random.seed(12)
        num_observations = 5000
        x1 = np.random.multivariate_normal([0, 0], [[1, .75], [.75, 1]], num_observations)
        x2 = np.random.multivariate_normal([1, 4], [[1, .75], [.75, 1]], num_observations)
        simulated_separableish_features = np.vstack((x1, x2)).astype(np.float32)
        simulated_labels = np.hstack((np.zeros(num_observations),
                                      np.ones(num_observations)))
        return (simulated_separableish_features[:, 0], simulated_separableish_features[:, 1], simulated_labels)

复制代码

生成数据。dom

若是是随机值的话,多试几回有可能方向是错误的。这里能跑这么快是由于初始值设置的好。 对于牛顿法的这个问题仍是有改进方法的。对于 \alpha步长的算法研究是能够解决这个问题。既然有时候牛顿法降低的方向或许是不对的,可使用一个阻力来阻挡一下,减慢速度,等到下一次能够从新选择。因此接下来就是对于 \alpha步长的研究。 ###线搜索 步长最小化意味着 f(\alpha) = minimizef(x_k + \alpha d_k)d_k就是方向了,这里的 x_k,d_k都是常数,因此能够看作是一个 \alpha的函数。 #####二分搜索算法 target function: \alpha = argmin_{\lambda>0}f(x+\lambda p) 假设 g(\alpha)=f(x_k + \alpha d_k),若是是一个凸函数,须要 g'(\alpha) = 0,由于 g'(\lambda) = \nabla f(x_k + \lambda p)p,p就当前的梯度降低的方向,因此向着梯度的反方向得: g'(0) < 0,假如咱们知道一个 g'(\lambda)>0,那么就能够知道确定有一个 x_k使得 g'(x_k) = 0,x_k \in (0, x_k),因此二分搜索的流程很简单: ①k = 0, \alpha_0 = 0\alpha_1 = am = \frac{\alpha_0 + \alpha_1}{2},若是 g'(m) > 0,\alpha_1 = m,k+=1 若是 g'(m) < 0,\alpha_0 = m,k+=1 若是 g'(m) = 0,end,不然就继续循环。 二分搜索也叫精确搜索,由于这样不断的迭代下去是能够很精确的搜索到一个合适的值,使用二分查找法来求步长的计算复杂度很高, 由于在最小化f(x)f(x)的每次迭代中咱们都须要执行一次线搜索, 而每次线搜索都要用上述的二分查找算法,由于这个精确的值在一些凸函数里面会是很小很小的,甚至后面多加一个小数,这样算起来复杂度很是高的。咱们能够在牺牲必定的精度的条件下来加快计算速度, 回溯线搜索是一种近似线搜索算法。 #####回溯线搜索 对于上诉二分法的问题,这里提出了改进,牺牲掉一些精度,求一个大概的范围就好,因而有:

f(x_k+\alpha p_k) <= f(x_k) + c_1 \alpha \nabla f_k^Tp_k$$上诉条件称为**充分降低条件**,$c_1 \in (0, 1)$
$$f(x_k+\alpha p_k) >= f(x_k) + c_1 (1-\alpha) \nabla f_k^Tp_k$$这两个准则就是Armijo-Goldstein准则,①目标函数值应该有足够的降低;②一维搜索的步长α不该该过小。这两个目标其实很明显,足够的降低其实就是式子1:![](https://upload-images.jianshu.io/upload_images/10624272-4e8a494e2452c6e8.png?imageMogr2/auto-orient/strip%7CimageView2/2/w/1240)能够看到接受这个条件有两个区间,有时候会选择到第一个区间的内容,也就是第一个区间的内容,因此第二条式子的做用就是舍得步长不要过小了。
>>####为何$c_1 \in (0, \frac{1}{2})$
 >> **①对于$c_1$的选择是必定要在(0,0.5)之间的,若是没有这个条件,会影响算法的超线性收敛性。这个怎么收敛的我就不知道了,还得看看其余的凸优化书籍。
②$f(x_k + \alpha p_k)$作泰勒展开,$f(x_k + \alpha p_k) = f(x_k) + \alpha g_k^T d_k+o(\alpha)$,去掉高阶无穷小,$f(x_k + \alpha p_k) = f(x_k) + \alpha g_k^T d_k$,和上式就是相差了一个$p$而已,这样就保证了小于的这个条件。
③因为$p \in (0, \frac{1}{2})$,因此$(1-p) \in (\frac{1}{2}, 1)$,并且$g_k^Td_k<0$,因此有$\alpha_k(1-p)g_k^Td_k < \alpha_k p g_k^Td_k < 0$,这就保证了a不能过小**

综上,这就获得了Armijo-Goldstein准则。后面用的Armijo搜索只是用了第一条式子进行搜索。
#####阻尼牛顿法
回到刚刚的牛顿法,提到有可能会存在Hession Matrix不是正定的状况,咱们能够加一点阻碍,也就是加上一个步长的限制,这个步长就是由Arimji搜索获得,这里只是会用到第一条准则,第二三条先不用。
```
class DampedNewton(object):

    def __init__(self, feature, label, iterMax, sigma, delta):
        self.feature = feature
        self.label = label
        self.iterMax = iterMax
        self.sigma = sigma
        self.delta = delta
        self.w = None

    def get_error(self, w):
        return (self.label - self.feature*w).T * (self.label - self.feature*w)/2

    def first_derivative(self):
        m, n = np.shape(self.feature)
        g = np.mat(np.zeros((n, 1)))
        for i in range(m):
            err = self.label[i, 0] - self.feature[i, ]*self.w
            for j in range(n):
                g[j, ] -= err*self.feature[i, j]
        return g

    def second_derivative(self):
        m, n = np.shape(self.feature)
        G = np.mat(np.zeros((n, n)))
        for i in range(m):
            x_left = self.feature[i, ].T
            x_right = self.feature[i, ]
            G += x_left * x_right
        return G

    def get_min_m(self, d, g):
        m = 0
        while True:
            w_new = self.w + pow(self.sigma, m)*d
            left = self.get_error(w_new)
            right = self.get_error(self.w) + self.delta*pow(self.sigma, m)*g.T*d
            if left <= right:
                break
            else:
                m += 1
        return m

    def newton(self):
        n = np.shape(self.feature)[1]
        self.w = np.mat(np.zeros((n, 1)))
        it = 0
        while it <= self.iterMax:
            g = self.first_derivative()
            G = self.second_derivative()
            d = -G.I * g
            m = self.get_min_m(d, g)
            self.w += pow(self.sigma, m) * d
            if it % 100 == 0:
                print('it: ', it, ' error: ', self.get_error(self.w)[0, 0])
            it += 1
        return self.w

```
求海塞矩阵换了一个写法,好看一点,可是效果都是同样的。
线性回归的损失函数:$$loss = \frac{1}{2}\sum_{i=1}^{m}(y^i-\sum_{j-0}^{n-1}w_j*x_j^i)^2

线性回归的一阶导数:$$\frac{dloss}{dw_j} = -\sum_{i=1}^{m}(y^i - \sum_{j=0}^{n-1}w_j*x_j^i)x_j^i$$ 线性回归的二阶导数:$$\frac{dloss}{dw_jdw_k} = \sum_{i=1}^m(x_j^ix_k^i)$$ 对照上面的公式是能够一一对应的。对于Armrji搜索只是用了第一条公式:$$f(x_k+\alpha p_k) <= f(x_k) + c_1 \alpha \nabla f_k^Tp_k$$一旦超过这个限制就能够退出了。 机器学习

效果其实仍是能够的。固然,也能够三条准则都使用上,只不过使用一条通常也足够了。 到这里线搜索尚未完。 #####Wolfe-Powell准则 有Armijo-Goldstein准则仍是不够的,以下图:
能够看到(a,b),(c,d)就是选择的区间,可是很明显这些区间已经避开了最低的点,固然这不是必定会,可是有这个可能,为了解决这个问题就出现了Wolfe-Powell准则。 Wolfe-Powell准则也有两个表达式,第一个表达式就是$$f(x_k+\alpha p_k) <= f(x_k) + c_1 \alpha \nabla f_k^Tp_k$$第二条式子就去掉了,由于它会很容易的把极值点排除在外,因此并不须要这个条件。添加的就是:$$\nabla f(x_k + \alpha_k d_k)^Td_k >= c_2 g_k^Td_k,c_2 \in (c_1, 1)$$这条式子的直观解释就是, 当前可选点的斜率要大于初始点斜率的c_2倍。,而初始点( α=0 的点)处的切线是比 e 点处的切线要“斜”的,因为 c_2∈(ρ,1) ,使得 e 点处的切线变得“不那么斜”了:
e和f后面的点就是知足的了, f(x_k+\alpha p_k) >= f(x_k) + c_1 (1-\alpha) \nabla f_k^Tp_k这个条件是没有的了。这个条件仍是不太好,有时候会给出更严谨的一个条件: |\nabla f(x_k + \alpha_k d_k)^Td_k| >= -c_2 g_k^Td_k,c_2 \in (c_1, 1),其实就是两边都给绝对值,可是右边的原本就是一个负数了,加上符号就是正数。这个条件就更强了:
直接限定在了(e,g)和(f,h)范围以内了。这就是整个Wolfe-Powell准则。 #####线搜索的总结 线搜索到这里基本就结束了。 寻找一个最佳步长就是要求优化以后的函数值是要最小的,因而给出的目标函数就是:\alpha = argmin_{\lambda>0}f(x+\lambda p)。解决这个问题有两种方法,一种就是二分法的查找,这种方法可能很是的精确搜索到一个最佳的值,可是他的计算复杂度有点高;有时候咱们其实不太须要太过精确的值,咱们只是须要一个大概模糊的就行了,因而出现了回溯搜索,这是属于模糊搜索,给出的就是Armijo-Goldstein准则,可是这两个准则也有很差的地方,有时候会把极值点排除在外,因而引入曲率条件,把Armijo-Goldstein准则的第二条式子换掉,就出现了Wolfe-Powell准则,条件再强一点两边加上绝对值,这样就出现了最终的Wolfe-Powell准则。 对于步长的研究并不针对某一个算法,对于优化降低的算法均可以使用,梯度降低,牛顿法,拟牛顿法均可以用到,具备很强的普适性。梯度降低的步长也是能够经过这种方式进行选择最优步长,牛顿法用Armijo搜索的方法是能够获得全局牛顿法,也叫阻尼牛顿法,这样可使得迭代方向能够避免向错误的方向进行,增长点阻力。 这篇文章主要的研究对象仍是牛顿法,因此下面的三个算法分别就是DFP,BFGS,LBFGS。 ###DFP 前面的阻尼牛顿法解决的就是对于迭代方向的问题,可是对于计算复杂度这个问题尚未获得解决,由于矩阵求拟是一个很大的工程量,若是维度一可能是计算复杂度是很大的,因此拟牛顿法基本上都是构造一个和Hession矩阵类似的矩阵来替代。可是问题来了,用近似矩阵替代,效果可能会很差,给出以下证实: 首先由牛顿法能够获得: -\frac{\nabla f(x_i)}{H} = x - x_i \nabla f(x_i) = -(x-x_i)H \nabla f(x_i)d = (x-x_i)H * H^{-1}g \nabla f(x_i)d = -(x-x_i)H(x-x_i) 因为 \nabla f(x_i) d < 0,因此右边也是要求小于0,那么天然就是大于0了,这也就证实了Hession矩阵的正定性。在逐步寻找最优解的过程当中要求函数值是降低的,可是有时候Hession矩阵不必定是半正定的,可能会使得函数值不降反升。 拟牛顿法可使目标函数值沿降低方向走下去,而且到了最后,在极小值点附近,可以使构造出来的矩阵与Hesse矩阵“很像”了,这样,拟牛顿法也会具备牛顿法的二阶收敛性。 #####推导证实 涉及到Hession矩阵,天然就涉及到二阶导数矩阵了,Hession矩阵只是对于多变量二阶导数的一种描述方式。Taylor展开:$$f(x) = f(x_{i+1})+(x-x_{i+1})^T\nabla f(x_{i+1}) + \frac{1}{2}(x-x_{i+1})^TH_{i+1}(x+ x_{i+1})+o(x+x_{i+1})....$$仍是用二次导数作近似,不少函数均可以用二次导数作近似。去掉高次项两边求导数:$$\nabla f(x_i) = \nabla f(x_{i+1}) + H_{i+1}(x_i-x_{i+1}) \ H^{-1}\nabla f(x_i)=H^{-1}\nabla f(x_{i+1}) + (x_i-x_{i+1}) \ H^{-1}(\nabla f(x_i) - \nabla f(x_{i+1})) = (x_i - x_{i+1})$$ 上面那个方程就是拟牛顿方程。上面那个矩阵就是一个近似矩阵。近似矩阵有不少种,若是按照上面那种方式进行计算,复杂度没有降多少,因此比较常见的方法就是迭代计算,比较常见的迭代方式:$$H_{i+1} = H_i + E_i$$ 那么接下来问题来了, E要这么设计。通常给出的设计就是$$E = mvv^T+nww^T$$由于但愿每次迭代h能有一个微小的变化而不是巨大的变化,这样才有可能收敛。并且这个结构设计的也很简单,也是对称的。 代入上面的拟牛顿方程:$$H_{i+1}(\nabla f(x_{i+1})-\nabla f(x_i)) = (H_i+E_i)(\nabla f(x_{i+1}) - \nabla f(x_i)) \ (H_i + mvv^T+nww^T)(\nabla f(x_{i+1}) - \nabla f(x_i))=x_{i+1} - x_i \ H_i(\nabla f(x_{i+1})-\nabla f(x_i)) + mvv^T(\nabla f(x_{i+1})-\nabla f(x_i)) + nww^T(\nabla f(x_{i+1})-\nabla f(x_i)) = x_{i+1}-x_i \ support: q_i = (\nabla f(x_{i+1})-\nabla f(x_i)) ,s_i = x_{i+1}-x_i \ then: H_iq_i+mvv^Tq_i+nww^Tq_i = s_i \ H_iq_i+v(mv^Tq_i)+w(nw^T)q_i = s_i$$仔细看一下这个式子,右边的 s_i是一个nx1的向量,而 mv^Tq_i,nw^Tq_i均为实数,也就是点积。能够直接假设 v = s_i,w = H_iq_i,因而有 mv^Tq_i = 1, nw^Tq_i = -1。 代入上面的式子:$$H_{i+1} = H_i+\frac{s_is_i^T}{s_i^Tq_i}-\frac{H_iq_iq_i^TH_i}{q_i^TH_iq_i}$$一开始我看到这个推导过程,我是一脸懵逼的。原来数学家也有猜的时候, 这个过程和前面假设H的形式在凸优化里面,是用"很显然"来描述的,嗯,很显然!我就是他妈的没看出来显然在哪了,后面的LBFGS也是,真的很显然,可是这个我是真的不知道。看到网上的不少解释都离不开特殊状况,也就是为0或者是直接假设v_i = ms_i的,都差很少的解释。 已知初始正定矩阵H,从一个初始点开始(迭代),用式子  d_k = -H_kg_k 来计算出下一个搜索方向,并在该方向上求出可以使目标函数极小化的步长α,而后用这个步长,将当前点挪到下一个点上,并检测是否达到了程序停止的条件,若是没有达到,则用上面所说的[13]式的方法计算出下一个修正矩阵H,并计算下一个搜索方向……周而复始,直到达到程序停止条件。 因此整一个流程: ①给定一个初值H_0 = I ②搜索方向:d_k = -H_k *g_k ③利用搜索获得步长\alpha,可使用上面提到的Armrji搜索或者等等的改进方法。 ④更新一波 ⑤计算y_k = g_{k+1}-g_{k}H_{i+1} = H_i+\frac{s_is_i^T}{s_i^Tq_i}-\frac{H_iq_iq_i^TH_i}{q_i^TH_iq_i},转回去继续更新。 DFP仍是有缺点的:
具体主要内容就是说若是Hession矩阵的是错误估计的,大家BFGS会在不多的迭代回到正确的方向,可是DFP在这方面的效果不明显。至于为何不明显,我从公式还不能直接看出来。既然提到了BFGS,下面就是BFGS了。 ###BFGS BFGS和DFP实际上是属于对偶的解法,一个直接求海赛矩阵逆矩阵(DFP),另外一个就是求海赛矩阵(BFGS)。仍是使用拟牛顿公式: q_i = H_is_i,而DFP使用的是 H_i^{-1}q_i = s_i,直接求出来的就是逆矩阵了。推导过程其实很简单,和上面DFP相似,差异不大,甚至除了换一下没什么差异的。 #####推导证实 假设都和上面DFP的同样。仍是 H_{i+1} = H_{i} + E_i,E就是矫正矩阵。E的形式和以前的是同样的 E = mvv^T+nww^T。代入上面的公式:$$q_i = H_{i+1}s_i \ q_i = (H_i + mvv^T + nww^T)s_i \ q_i = H_is_i+(mv^Ts_i)v+(nw^Ts_i)w$$一样假设 v = q_i,w = H_is_i,一样选取特殊状况, v = q_i,w = H_is_i,m = \frac{1}{v^Ts_i},n = -\frac{1}{w^Ts_i},代进矫正函数的表达式:$$H_{i+1} = H_i+\frac{vv^T}{v^Ts_i} - \frac{ww^T}{w^Ts_i} \ H_{i+1} = H_i+\frac{q_iq_i^T}{q_i^Ts_i} - \frac{H_is_is_i^TH_i^T}{s_i^TH_is_i} $$仔细看一些这个式子,和DFP那个式子是换了一下位置而已,因此这个过程和步骤真的没有什么好说的。按照常规套路,通常是这样: 因此整一个流程: ①给定一个初值H_0 = I ②搜索方向:d_k = -H_k^{-1} *g_k ③利用搜索获得步长\alpha,可使用上面提到的Armrji搜索或者等等的改进方法。 ④更新一波 ⑤计算y_k = g_{k+1}-g_{k}H_{i+1} = H_i+\frac{q_iq_i^T}{q_i^Ts_i} - \frac{H_is_is_i^TH_i^T}{s_i^TH_is_i},转回去继续更新。 然而,若是是这样,复杂度仍是存在的,仍是得求个导数啊。因此这样的作法天然是不行的,由于拟牛顿法的改进有一个很重要的诱因就是计算复杂度的问题了。 对于求逆矩阵,有一个很是重要的公式——Sherman-Morrison公式:$$(A+uv^T)^{-1}=A^{-1}-\frac{A^{-1}uv^TA^{-1}}{1+v^TA^{-1}u}$$ 用Sherman-Morrison公式改造一下:H_{k+1}^{-1}=(I-\frac{s_kq_k^T}{q_k^Ts_k})H_k^{-1}(I-\frac{q_ks_k^T}{q_k^Ts_k})+\frac{s_ks_k^T}{q_k^Ts_k}至因而怎么改造的,有兴趣本身看吧,原谅我并无看懂,懂了我也想不出来。 根据这个结果改造一下: 因此整一个流程: ①给定一个初值H_0 = I ②搜索方向:d_k = -H_k^{-1} *g_k ③利用搜索获得步长\alpha,可使用上面提到的Armrji搜索或者等等的改进方法。 ④更新一波 ⑤计算y_k = g_{k+1}-g_{k}H_{k+1}^{-1}=(I-\frac{s_kq_k^T}{q_k^Ts_k})H_k^{-1}(I-\frac{q_ks_k^T}{q_k^Ts_k})+\frac{s_ks_k^T}{q_k^Ts_k},转回去继续更新。 完美的达到了下降复杂度的目标。可是我仍是不知道BFGS为何相对DFP来讲对于纠正方向要更快一点。这两个是对偶式子,难到是由于用了Sherman-Morrison公式吗? #####代码实现

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from DataProcesser import DataProcesser

def bfgs(feature, label, lam, maxCycle):
    n = np.shape(feature)[1]
    w0 = np.mat(np.zeros((n, 1)))
    rho = 0.55
    sigma = 0.4
    Bk = np.eye(n)
    k = 1
    while (k < maxCycle):
        print('Iterator: ', k, ' error: ', get_error(feature, label, w0))
        gk = get_gradient(feature, label, w0, lam)
        dk = np.mat(-np.linalg.solve(Bk, gk))
        m = 0
        mk = 0
        while (m < 20):
            newf = get_result(feature, label, (w0 + rho**m*dk), lam)
            oldf = get_result(feature, label, w0, lam)
            if (newf < oldf + sigma * (rho ** m)*(gk.T*dk)[0, 0]):
                mk = m
                break
            m += 1
        #BFGS
        w = w0 + rho**mk * dk
        sk = w-w0
        yk = get_gradient(feature, label, w, lam) - gk
        if (yk.T * sk > 0):
            Bk = Bk - (Bk * sk * sk.T * Bk) / (sk.T * Bk * sk) + (yk * yk.T) / (yk.T * sk)
        k = k + 1
        w0 = w
    return w0

def get_error(feature, label, w):
    return (label - feature * w).T*(label - feature * w)/2

def get_gradient(feature, label, w, lam):
    err = (label - feature * w).T
    left = err * (-1) * feature
    return left.T + lam * w

def get_result(feature, label, w, lam):
    left = (label - feature * w).T * (label - feature * w)
    right = lam * w.T * w
    return (left + right)/2

复制代码

最主要的部分也就是公式那部分了,没有什么好讲的。实现效果: 函数

###L-BFGS L-BFGS的出现仍是为了节省,但此次是为了节省内存。每一次存储高纬度的数据须要的空间太多了,因此LBFGS的主要做用就是减小内存的开销。再也不完整的存储整一个矩阵了,它对BFGS算法作了近似,而是存储计算时所须要的\{s_k\},\{q_k\},而后利用哪两个向量的计算来近似代替。对于\{s_k\},\{q_k\}也不是存储全部的,而是只保存那么后面的m个,也就是最新的m个就好。这样就把空间复杂度O(N^2)降到O(mN)。 毫无疑问,LBFGS是从BFGS演化而来的。对于逆矩阵的公式:$$H_{k+1}^{-1}=(I-\frac{s_ky_k^T}{y_k^Ts_k})H_k^{-1}(I-\frac{y_ks_k^T}{y_k^Ts_k})+\frac{s_ks_k^T}{y_k^Ts_k}$$ 假设p_k = \frac{1}{y_k^Ts_k},v_k = 1-p_ky_ks_k^T,改造一下上式:$$H_{k+1}^{-1} = v_k^TH_k^{-1}v_k+p_ks_ks_k^T$$ 假设:D_k = H_k^{-1} 递推一下:oop

D_1 = v_0^TD_0v_0+p_0s_0s_0^T\\
D_2 = v_1^TD_1v_1+p_1s_1s_1^T\\
   = v_1^T(v_0^TD_0v_0+p_0s_0s_0^T)v_1+p_1s_1s_1^T\\
= v_1^Tv_0^TD_0v_0v_1+v_1s_0s_0^Tv_1+p_1s_1s_1^T
D_3 = v_2^Tv_1^Tv_0^TD_0v_0v_1v_2+v_2^Tv_1^Tp_0s_0s_0^Tv_1v_2+v_2^Tp_1s_1s_1^Tv_2+p_2s_2s_2^T

根据推论,通常的有:$$D_{k+1} = (v_k^Tv_{k-1}^T...v_1^Tv_0^T)D_0(v_0v_1...v_{k-1})\ +(v_k^Tv_{k-1}^T...v_2^Tv_1^T)(p_0s_0s_0^T)(v_1v_2...v_{k-1}v_k)\ +(v_k^Tv_{k-1}^T...v_3^Tv_2^T)(p_1s_1s_1^T)(v_2v_3...v_k-1v_k)\ +...\ +(v_k^Tv_{k-1}^T)(p_{k-2}s_{k-2}s_{k-2}^T)(v_{k-1}v_{k})\ +v_k^T(p_{k-1}s_{k-1}s_{k-1}^T)v_k\ +p_ks_ks_k^T$$ 能够看的出,当前的海赛矩阵的逆矩阵是能够由\{s_i,q_i\}给出的,因此直接存储对应的s_i.q_i便可,上面提到过为了压缩内存,能够只是存储后面的m个,因此更新一下公式:$$D_{k+1} = (v_k^Tv_{k-1}^T...v_{k-m+1}^Tv_{k-m}^T)D_0(v_{k-m}v_{k-m+1}...v_{k-1})\ +(v_k^Tv_{k-1}^T...v_{k-m+2}^Tv_{k-m+1}^T)(p_0s_0s_0^T)(v_{k-m+1}v_{k-m+2}...v_{k-1}v_k)\ +(v_k^Tv_{k-1}^T...v_{k-m+3}^Tv_{k-m+2}^T)(p_1s_1s_1^T)(v_{k-m+2}v_{k-m+3}...v_k-1v_k)\ +...\ +(v_k^Tv_{k-1}^T)(p_{k-2}s_{k-2}s_{k-2}^T)(v_{k-1}v_{k})\ +v_k^T(p_{k-1}s_{k-1}s_{k-1}^T)v_k\ +p_ks_ks_k^T$$ 公式长是有点长,可是已经有大神给出了一个很好的算法解决这个问题,two-loop-recursion算法: y_k=\nabla f(x_k)-\nabla f(x_{k-1}) q = \nabla f(x_k) for (i=1...m)do \alpha_i = p_{k-i}s_{k-i}^Tq q = q-\alpha_iy_{k-i} end for r = H_{k}^0q for(i=m...1)do \beta = p_{k-1} y_{k-1}^T r r = r+s_{k-i}(\alpha_i - \beta) return r 这个式子到底行不行呢?证实一下理论: 学习

这样就证实这个算法的正确性。然而其实我根本不关心这个算法正确性,我只是想知道 ###这是怎么想出来的,说实话第一眼看根本没有get到这个算法就是实现了LBFGS,因此若是有大神知道麻烦私信我!渣渣感激涕零。 按照上面算法获得方向以后就可使用线搜索获得合适的步长最后结合更新了。 #####代码实现 前面求梯度都同样的,就是后面的更新有不一样:

def lbfgs(feature, label, lam, maxCycle, m = 10):
    n = np.shape(feature)[1]
    w0 = np.mat(np.zeros((n, 1)))
    rho = 0.55
    sigma = 0.4

    H0 = np.eye(n)
    s = []
    y = []

    k = 1
    gk = get_gradient(feature, label, w0, lam)
    dk = -H0 * gk

    while (k < maxCycle):
        print('iter: ', k, ' error:', get_error(feature, label, w0))
        m1 = 0
        mk = 0
        gk = get_gradient(feature, label, w0, lam)
        while (m1 < 20):
            newf = get_result(feature, label, (w0 + rho ** m1 * dk), lam)
            oldf = get_result(feature, label, w0, lam)
            if newf < oldf + sigma * (rho ** m1) * (gk.T * dk)[0, 0]:
                mk = m1
                break
            m1 = m1 + 1
        w = w0 + rho ** mk * dk
        if k > m:
            s.pop(0)
            y.pop(0)
        sk = w - w0
        qk = get_gradient(feature, label, w, lam)
        yk = qk - gk
        s.append(sk)
        y.append(yk)

        t = len(s)
        a = []

        for i in range(t):
            alpha = (s[t - i -1].T * qk) / (y[t - i - 1].T * s[t - i - 1])
            qk = qk - alpha[0, 0] * y[t - i -1]
            a.append(alpha[0, 0])
        r = H0 * qk

        for i in range(t):
            beta = (y[i].T * r) / (y[i].T * s[i])
            r = r + s[i] * (a[t - i - 1] - beta[0, 0])
        if yk.T * sk > 0:
            dk = -r
        k = k + 1
        w0 = w
    return w0

复制代码

中间还要有一个判断降低方向的过程。

数据比较少,因此效果差异都不大。 ###总结 首先提到的是梯度降低,梯度降低算法虽然很简单,可是降低的方向会有所误差,可能胡局部不稳定,速度不会特别快,可是最终是会到达终点。因而改进一下,梯度降低是一阶拟合,那么换牛顿法二阶拟合,可是牛顿法问题来了,迭代的方向有多是错误的,因此改进一下,加点阻力,就算是不许确的,用linear search也能够调整一下。可是对于阻尼牛顿法仍是存在了计算复杂度的问题,因而改进一下,用DFP直接作近似逆矩阵,达到了下降复杂度的问题,可是对于方向梯度的问题仍是不是特别好,因而又出来了BFGS,用了比较牛逼的Sherman-Marrion公式求出了逆矩阵。问题又来了,每次都存储这么大的矩阵,有没有办法下降一些呢?因而改进了一下,就出现了LBFGS,用两个nx1的矩阵来表示逆矩阵,而且只存储后m个更新的内容,改进一下出现了two-loop-recursion算法。后面又继续改进了一下。这部分的知识比较靠近数值优化我心脏已经受不了了。

####GitHub代码:github.com/GreenArrow2…

相关文章
相关标签/搜索