第二十一节 牛顿法和L-BFGS求函数最优解
这一节中，我们讲解一个新的求函数最优化的方法就是L-BFGS。以下是本节目录。
目录
1-L-BFGS算法简介
我们知道算法在计算机中运行的时候是需要很大的内存空间的.就像我们解决函数最优化问题常用的梯度下降,它背后的原理就是依据了泰勒一次展开式.泰勒展开式展开的次数越多,结果越精确,没有使用三阶四阶或者更高阶展开式的原因就是目前硬件内存不足以存储计算过程中演变出来更复杂体积更庞大的矩阵.L-BFGS算法翻译过来就是有限内存中进行BFGS算法,L是limited memory的意思.那算法为什么叫BFGS呢,请看下图:
上图中从左到右依次是Broyden,Fletcher,Goldfarb,Shanno.四位数学家名字的首字母是BFGS,所以算法的名字就叫做BFGS算法.接下来我们就一起来学习BFGS算法的内容。
2-牛顿法求根问题
我们先来回顾下牛顿法求根问题，比如求1元2次方程的根公式为，我们通常管这种形式的根叫解析根。所谓解析解就是你不用给我具体的值，就是一个公式。3次方程也是有解析解的，但是当函数达到5次方以上，就不好找解析解了，对于这种复杂的函数，很遗憾我们不能找到它的全部根，但是至少有办法找到它的一个根。
我们看一个对于一元函数的例子：
对于一个无法解出解析解的函数来说，现在是一元函数，在只有一个x的情况下，我最终想找到x令y=0，即函数的根，怎么找到它？牛顿和另外一个人同时分别发现，假如这个函数是连续可导的，我随机出来一个x1，总能求它在这点x1的导数，导数是一个实数，它是能代表切线的斜率，那么我们就在这个x1点上画一个原函数切线，这个切线一定会与x轴相交，除非特别倒霉，你一步就随机到它的驻点上了，也不用求了，你找的就是驻点。不倒霉的情况下一定会使x1点的切线和x轴有个交点，我们命名为x2。 点完之后，又可以找一下x2对应在函数上的点，再画一条切线找到了x3，这时候相比x1，距离跟的位置来说比较近，然后发现经过有限次数的迭代之后，怎么迭代都不会变化。它还是会产生震荡，经过震荡之后，最终会收敛在某一个点上，而这个点就是函数的根。所以说牛顿法利用几何直觉就是在找到某一个函数与x轴的交点。
我们分析下几何原理：
函数本身f(x)也是已知的，那么f(x1)就可以计算。x2是x1切线的位置，我们建立一个关于x2，x1的解析式，三角形这条边是x1-x2，另外一条边是f(x1)，如果用f(x1)/ x1-x2就应该是这条线的斜率，斜率又等于x1这点导数的值。所以x2，x1建立一个关系，x2 = x1-f(x1)/ f`(x1)，此时x2就可求了。我们大的思路想从x1找到x2，从x2找到x3，最终找着发现xn，xn+1它俩没区别的时候xn就是函数的根。因此建立了一个xn，xn-1的之间关系的迭代公式：
xn=xn-1- f(xn-1)/ f`(xn-1)
我们在求数值解的时候，虽然我们求根法是求损失函数最小值，当最小值只差一点点的时候，对模型的预测结果影响不大，但是它在底部会震荡很久。 对于这种在底部频繁震荡的情况下，通常会设置超参tolerance，当xk，xk+1的变化已经小于你设的阈值的时候，我就认为它已经收敛了，而不去非得接近它，驻点数本身就是不那么可控的。设置超参tolerance，它几乎在不牺牲精度的情况下来对我们这个函数最优化算法进行。
总结下牛顿法求根的流程：
.1、已知函数f(x)的情况下随机产生x0。
2、由已知的x0按照 xn=xn-1- f(xn-1)/ f`(xn-1)公式进行n次迭代。
3、当迭代结果xn与上一次迭代结果xn-1相同或小于一定阈值时，本次的结果即为函数f(x)的根。
3-牛顿法求驻点问题
既然利用上述办法找到任意函数的根，能不能借这个找到函数最小值呢？我们可以利用驻点的知识，虽然驻点不一定是极值点，极值点也不一定是驻点，但大部分情况下两个值相等，驻点也就是导数为零的点，即原函数最小值的点。 所以我们需要找到导函数为零的点，实际上也就是导函数的根。导函数是一个函数，它是用来计算一个函数导数的工具，给我一个x，我就算这个原函数在这个点的导数是多少的这么一个函数。如果我们能把原函数的导函数写出来，进而求导函数的根，就解决了求原函数最小值的问题。
既然都是求根问题，我们可不可以利用上面的方法呢？肯定是可以的，只不过现在的函数是我们要求的导函数而已。我们看下刚才原函数求根的迭代公式：xn+1 = xn-f(xn)/ f`(xn)。对于导函数来说，唯一的区别就是f(xn)表达式不同，所以我们此时的原函数就是f`(xn)，此时f`(xn)的导函数就是f``(xn)。最终的迭代公式就是：
xk = xk-1-f'(xk-1)/ f''( xk-1
通过这个东西，我最终找到xn不再是f(x)的根，而是f(x)导数的根，也就是驻点。
总结下利用牛顿法求函数的驻点过程：
1、当函数f(x) 的一阶导数 f’(x) = 0 时点(x,f(x))为函数f(x)的驻点
2、求某函数的驻点即为求该函数的导函数的根,同样可以利用牛顿法进行求解
3、对于f(x) 函数来说 迭代公式为 xk = xk-1 - f’(xk-1)/f’’(xk-1)
4-牛顿法求驻点的本质
牛顿法求驻点本质实际上是二阶泰勒展开公式。我们先来回顾下什么是泰勒展开。所谓泰勒展开就是把任意一个复杂函数在某个点附近，用一个有限长度的多项式来拟合，那么通常多项式在xk点附近展开是φ（x）=f(x)+ f’(x-xk)+ 1/2！f''(xk)*(x-xk)^2+1/3！f'''(xk)(x-xn)^3...，一直往上加，越加越像原函数。 我们看下一阶展开就是φ（x）=f(x)+ f’(x-xk)。它本质就是y=ax+b，一条直线，只不过a和b里面杂糅了很多关于xk函数的数值，a和b是通过xk这个点的函数还有一阶导函数算出来的，所以一阶泰勒展开就是在xk附近用一条直线尽量的去拟合原函数。什么叫做xk尽量的去拟合原函数？意思是在xk点附近，我不关注其它地方，我就在你周围画一条直线，你能跟我原先xk点的附近重合的越多越好。而二阶泰勒展开式是ax2+bx+c，只不过a，b，c是通过这些导数什么杂糅在一起算出来，无论这些东西等于多少，它始终是一个实数，它一定就是这个形势，这个形式是一个抛物线。
我们看下任意函数在 xk点附近的二阶泰勒展开公式为：
该公式表达的函数φ（x） 的几何意义为：通过2次函数对于原函数的最佳拟合。当φ'（x）=0时：
解释下这个公式怎么来的：对x求导，没有对xk求导的，xk是个已知数，第一项f(xk)没有x，求导等于0。第二项对x的求导，把它展开变成x*f’（xk）就剩这么一个东西，f’（xk）它是个数，它不是x函数。第三项对x求导结果是1/2 f''(xk).2(x-xk) 。综合以上结果就得到φ'（x）=0的解析式，也就是函数二阶泰勒展开所拟合出来的抛物线的最小值。
把上述φ'（x）=0的解析式展开成x= xk -f’(xk)/f’’(xk)。发现刚好和之前牛顿法求驻点xk = xk-1 - f’(xk-1)/f’’(xk-1)的迭代公式一致。所以下一次xxk的迭代值就是以它拟合出来的抛物线的最小值对应的x来作为我原函数根的下一次迭代值。
我们画图示意牛顿法求驻点的本质：
因为一阶泰勒展开是一条直线它没有最小值，而抛物线有最小值，实际上它就是在拿拟合出来的抛物线的最小值当成更好一点的x作为下一次的起始点，下次到这新的点，就又找了一个抛物线，以此类推，类似于我们每次拿一个碗形的曲线去套原先的曲线，所以它能更好的拟合在某一点附近的真实曲线情况，正因为如此，这个东西带来的好处就是牛顿法比梯度下降在底部震荡次数小很多，它能够让我们在求函数最优解的过程中有更少的收敛次数。
5-多元函数利用牛顿法求驻点
上面所说的都是针对一元的，而机器学习里面都是多元的，所谓的x在机器学习里面是谁？你要优化损失函数，损失函数是跟着w来的，我们一定要了然一个事，在机器学习的训练集里，x是已知数，你唯一未知数就是w。x都是w的系数，它在损失函数里面无论什么情况，最终都会变成w的系数。那么这里面的w实际上是机器学习里面的未知数x，w有多少个，往往不止一个。对于多元函数求驻点怎么求？实际在多元函数中，一阶导数映射到多元函数里就是梯度：即一阶导数f'(x) ----->梯度
梯度就是它的所有一阶导数给它写一块去。比如你原来就一个x，求个导就完事了。你现在有十个x，你求谁都不合适，你就干脆把十个逐个求个叫偏导，然后把这第一个x求偏导的函数放在这个向量的第一位，第二个求偏导的函数放在向量第二位
，最后一个求偏导放在最后一位，组成一个向量。这个向量里面存的目前来讲还是一个函数，这个叫做梯度函数向量，而梯度是不是指的是某个点的梯度，就像导数指的是某一个点的导数一样，那么梯度里边存的就不再是函数，而是把某组x带到这一组函数里面，去得到了那一组实数的向量值，所有你只要说梯度，虽然它里面写的是字母，但实际上知道这是一个真真切切的实数。 它是一个实数向量。 那么梯度是什么情况？假设你有十个w，梯度就有10个元素，因为要对10个未知的w逐个求偏导放到向量里面。
对第一个w求偏导，它得的结果是一个导函数，它们通常我们在写梯度的时候，有时候会这么写，这个x代表一个向量，它并不是向量中的某个元素，这个x就已经包含了若干个x了，x1到xn，k就代表第几代x。 这个梯度的意思就是说你先把这些偏导函数搁在这之后，然后再对当前这一代xk的具体的值带到这里来求出的结果，这个东西叫梯度。涉及到梯度，就意味着有一定的方向。
多元函数中，我们之前的求驻点的迭代公式xk = xk-1 - f’(xk-1)/f’’(xk-1)，f’(x)在多元函数中就是它对应的的梯度g（x） ，f’’( x)叫做Hessian矩阵对应关系为：一阶导数f''(x) ----->Hessian矩阵
这个矩阵是怎么求？就是它要求两次导，假如说原来对x1求导，求了之后还对x1求一次导，它还需要再对x1到xn分别求一次导，这样就形成了一个巨大的偏导数的矩阵，这个矩阵的情况应该是N乘N的，它的对角线元素是对同一个自变量求两次，而其它是各种各样的排列组合。这个就是二阶导数Hessian矩阵。
所以我们的迭代公式演变为：从到
其中gk就是多元函数一阶导，Hk就是多元函数二阶导。gk原来在分子中，Hessian矩阵应该是除以二阶导数，所以变成了矩阵的逆。我们可以看一下Hk矩阵，逆矩阵还是一个N*N的，gk是一个N*1的矩阵，它俩相乘得到是一个N*1的向量，就是一个高高的向量， N行1列。xk本身也是一个高高的向量，它俩做减法完全没有问题。但是Hk矩阵的N*N就给我们带来了一个非常不喜欢的问题，这个矩阵太难运算的，当你有1000个维度的时候，就是1000×1000的矩阵，然后1万个维度的时候，10000×10000，维度就爆炸了。而牛顿法求驻点又是一个迭代算法,所以这个困难我们还要面临无限多次,导致了牛顿法求驻点在机器学习中无法使用.有没有什么解决办法呢。
6-BFGS算法
BFGS算法是通过迭代来逼近的算法.逼近的方式如下:
其中：。I是单位矩阵，
BFGS就是通过迭代来逼近Hk的逆矩阵矩阵,其中第一步的D矩阵是单位矩阵.所谓单位矩阵就是只有对角线元素为1，其余全为零的矩阵。根据之前的迭代公式：，除了第一步还是需要计算之前的逆矩阵H1，第一步迭代需要计算出x2=x1-H1*g1，s1=x2-x1，y1=g2-g1。因为有x2，所以g2是其梯度，也可求，有了s1,y1便可以求D2来近似代替H2的逆矩阵，然后一步步迭代，求出驻点。
我们要通过牛顿求驻点法和BFGS算法来求得一个函数的根,两个算法都需要迭代,慢慢逼近函数根,经过k次迭代以后,所得到的解就是机器学习中目标函数导函数的根.这种两个算法共同迭代的计算方式,我们称之为On The Fly。
回顾一下梯度下降的表达式，在BFGS算法迭代的第一步x2=x1-D1*g1,单位矩阵与梯度g相乘,就等于梯度g，形式上同梯度下降的表达式是相同的。相当于学习率等于1的梯度下降，所以BFGS算法可以理解为从梯度下降逐步转换为牛顿法求函数解的一个算法.
虽然我们使用了BFGS算法来利用单位矩阵逐步逼近H矩阵,但是根据Dk+1的公式，每次计算的时候都要存储上一代的Dk矩阵,Dk矩阵有多大呢.假设我们的数据集有十万个维度(不算特别大),那么每次迭代所要存储D矩阵的结果是74.5GB.
我们无法保存如此巨大的矩阵内容,如何解决呢? L-BFGS 算法上阵。
7-L-BFGS算法
我们每一次对D矩阵的迭代,都是通过迭代计算sk和yk得到的.既然存不下D矩阵,那么我们存储下所有的sk和yk,想要得到D10就用单位矩阵同存储下的s1和y1到s10和y10计算就可以了.这样一个时间换空间的办法，有效节省了内存空间。
但是,仅仅是这样还是不够的,因为当迭代次数非常大的时候,我们的内存同样存不下.这个时候只能丢掉一些存不下的数据.假设我们设置的存储向量数为100,当s和y迭代超过100时,就会扔掉第一个s和y,存储s2到s101和y2到y101,每多一次迭代就对应的扔掉最前边的s和y。假如最后收敛次数是1000的话，只保留s901，y901从而计算出D901 ，然后再保留s902，y902计算出D902，依次根据s1000,y1000，从而算出D1000，按理说需要D900才能计算出D901 ，此时我们粗暴的将D901置为单位矩阵I，这样虽然损失了精度,但确可以保证使用有限的内存将函数的解通过BFGS算法求得到。
虽然L-BFGS算法是线性收敛，但是每次迭代的开销非常小，因此L-BFGS算法执行速度还是很快的，而且由于每一步迭代都能保证近似矩阵的正定，因此算法的鲁棒性还是很强的。