Optimization of Machine Learning

2018-09-18 本文已影响8人冒绿光的盒子

机器学习就是需要找到模型的鞍点，也就是最优点。因为模型很多时候并不是完全的凸函数，所以如果没有好的优化方法可能会跑不到极值点，或者是局部极值，甚至是偏离。所以选择一个良好的优化方法是至关重要的。首先是比较常规的优化方法：梯度下降。以下介绍的这些算法都不是用于当个算法，可以试用于能可微的所有算法。

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$
这个就是当前点的梯度，而函数的更新方向就是向着梯度的负方向进行的，所以最后是要减：
$\theta_j = \theta_j - \alpha\frac{1}{m}\sum_{i=0}^{m}(h(x_i)-y_i)x_i^j$

最优点很明显是在左边，而待优化点是

左边的就是成功拟合的，右边的那个就是方向错误，拟合没有成功的例子。所以牛顿法不一定会按照正确的方向拟合。上面牛顿法的式子是对于单变量的，如果是对变量，那么下面的二阶导要用到Hession矩阵。所以对于多变量的牛顿法：

如果是随机值的话，多试几次有可能方向是错误的。这里能跑这么快是因为初始值设置的好。
对于牛顿法的这个问题还是有改进方法的。对于

可以看到接受这个条件有两个区间，有时候会选择到第一个区间的内容，也就是第一个区间的内容，所以第二条式子的作用就是舍得步长不要太小了。

为什么 $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^i*x_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准则也有两个表达式，第一个表达式就是

e和f后面的点就是满足的了，

直接限定在了(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}$
递推一下：
$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代码：https://github.com/GreenArrow2017/MachineLearning/tree/master/MachineLearning/Linear%20Model/LogosticRegression/LinearRegaression

Optimization of Machine Learning

Gradient Descent

为什么 $c_1 \in (0, \frac{1}{2})$

阻尼牛顿法

Wolfe-Powell准则

线搜索的总结

DFP

推导证明

BFGS

推导证明

代码实现

L-BFGS

这是怎么想出来的，说实话第一眼看根本没有get到这个算法就是实现了LBFGS，所以如果有大神知道麻烦私信我！渣渣感激不尽。

代码实现

总结

GitHub代码：https://github.com/GreenArrow2017/MachineLearning/tree/master/MachineLearning/Linear%20Model/LogosticRegression/LinearRegaression

猜你喜欢

热点阅读