| Alink漫谈(十一) ：线性回归 之 L | 您所在的位置：网站首页 › lbfgs算法原理 › Alink漫谈(十一) ：线性回归 之 L |
|
Alink漫谈(十一) ：线性回归 之 L-BFGS优化
目录Alink漫谈(十一) ：线性回归 之 L-BFGS优化0x00 摘要0x01 回顾1.1 优化基本思路1.2 各类优化方法0x02 基本概念2.1 泰勒展开如何通俗推理？2.2 牛顿法2.2.1 泰勒一阶展开2.2.2 泰勒二阶展开2.2.3 高维空间2.2.4 牛顿法基本流程2.2.5 问题点及解决2.3 拟牛顿法2.4 L-BFGS算法0x03 优化模型 -- L-BFGS算法3.1 如何分布式实施3.2 CalcGradient3.3 AllReduce3.4 CalDirection3.4.1 预先分配3.4.2 计算方向3.5 CalcLosses3.6 UpdateModel3.7 OutputModel0x04 准备模型元数据0x05 建立模型0x06 使用模型预测0x07 本系列其他文章0xFF 参考
0x00 摘要
Alink 是阿里巴巴基于实时计算引擎 Flink 研发的新一代机器学习算法平台，是业界首个同时支持批式算法、流式算法的机器学习平台。本文介绍了线性回归的L-BFGS优化在Alink是如何实现的，希望可以作为大家看线性回归代码的Roadmap。 因为Alink的公开资料太少，所以以下均为自行揣测，肯定会有疏漏错误，希望大家指出，我会随时更新。 本系列目前已有十一篇，欢迎大家指点 Alink漫谈(一) : 从KMeans算法实现不同看Alink设计思想 Alink漫谈(二) : 从源码看机器学习平台Alink设计和架构 [Alink漫谈之三] AllReduce通信模型 Alink漫谈(四) : 模型的来龙去脉 Alink漫谈(五) : 迭代计算和Superstep Alink漫谈(六) : TF-IDF算法的实现 Alink漫谈(七) : 如何划分训练数据集和测试数据集 Alink漫谈(八) : 二分类评估 AUC、K-S、PRC、Precision、Recall、LiftChart 如何实现 Alink漫谈(九) ：特征工程 之 特征哈希/标准化缩放 Alink漫谈(十) ：线性回归实现 之 数据预处理 0x01 回顾到目前为止，已经处理完毕输入，接下来就是优化。优化的主要目标是找到一个方向，参数朝这个方向移动之后使得损失函数的值能够减小，这个方向往往由一阶偏导或者二阶偏导各种组合求得。 所以我们再次复习下基本思路。 1.1 优化基本思路对于线性回归模型 f(x) = w'x+e，我们构造一个Cost函数（损失函数）J(θ)，并且通过找到 J(θ) 函数的最小值，就可以确定线性模型的系数 w 了。 最终的优化函数是：min(L(Y, f(x)) + J(x)) ，即最优化经验风险和结构风险，而这个函数就被称为目标函数。 我们要做的是依据我们的训练集，选取最优的θ，在我们的训练集中让f(x)尽可能接近真实的值。我们定义了一个函数来描述 “f(x)和真实的值y之间的差距“，这个函数称为目标函数，表达式如下： \[J(θ)≈\frac{1}{2} \sum_{i≈1}^m(f_θ(x^{(i)} — y^{(i)})^2\\ \]这里的目标函数就是著名的最小二乘函数。 我们要选择最优的θ，使得f(x)最近进真实值。这个问题就转化为求解最优的θ，使目标函数 J(θ) 取最小值。 1.2 各类优化方法寻找合适的W令目标函数f(W) 最小，是一个无约束最优化问题，解决这个问题的通用做法是随机给定一个初始的W0，通过迭代，在每次迭代中计算目标函数的下降方向并更新W，直到目标函数稳定在最小的点。 不同的优化算法的区别就在于目标函数下降方向Dt的计算。下降方向是通过对目标函数在当前的W下求一阶倒数（梯度，Gradient）和求二阶导数（海森矩阵，Hessian Matrix）得到。常见的算法有梯度下降法、牛顿法、拟牛顿法。 梯度下降法直接采用目标函数在当前W的梯度的反方向作为下降方向。 牛顿法是在当前W下，利用二次泰勒展开近似目标函数，然后利用该近似函数来求解目标函数的下降方向。其中Bt为目标函数f(W)在Wt处的海森矩阵。这个搜索方向也称作牛顿方向。 拟牛顿法只要求每一步迭代中计算目标函数的梯度，通过拟合的方式找到一个近似的海森矩阵用于计算牛顿方向。 L-BFGS(Limited-memory BFGS)则是解决了BFGS中每次迭代后都需要保存N*N阶海森逆矩阵的问题，只需要保存每次迭代的两组向量和一组标量即可。Alink中，UnaryLossObjFunc是目标函数，SquareLossFunc 是损失函数，使用L-BFGS算法对模型进行优化。 即 优化方向由拟牛顿法L-BFGS搞定(具体如何弄就是看f(x)的泰勒二阶展开)，损失函数最小值是平方损失函数来计算。 0x02 基本概念因为L-BFGS算法是拟牛顿法的一种，所以我们先从牛顿法的本质泰勒展开开始介绍。 2.1 泰勒展开泰勒展开是希望基于某区间一点x_0展开，用一组简单的幂函数来近似一个复杂的函数f(x)在该区间的局部。泰勒展开的应用场景例如：我们很难直接求得sin(1)的值，但我们可以通过sin的泰勒级数求得sin(1)的近似值，且展开项越多，精度越高。计算机一般都是把 sin(x)进行泰勒展开进行计算的。 泰勒当年为什么要发明这条公式？ 因为当时数学界对简单函数的研究和应用已经趋于成熟，而复杂函数，比如：f(x) = sin(x)ln(1+x) 这种一看就头疼的函数，还有那种根本就找不到表达式的曲线。除了代入一个x可以得到它的y，就啥事都很难干了。所以泰勒同学就迎难而上！决定让这些式子统统现出原形，统统变简单。 要让一个复杂函数变简单，能不能把它转换成别的表达式？泰勒公式一句话描述：就是用多项式函数去逼近光滑函数。即，根据“以直代曲、化整为零”的数学思想，产生了泰勒公式。 泰勒公式通过把【任意函数表达式】转换（重写）为【多项式】形式，是一种极其强大的函数近似工具。为什么说它强大呢？ 多项式非常【友好】，三易，易计算，易求导，易积分 几何感觉和计算感觉都很直观，如抛物线和几次方就是底数自己乘自己乘几次 如何通俗推理？泰勒公式干的事情就是：使用多项式表达式估计（近似） 在 附近的值。 当我们想要仿造一个东西的时候，无形之中都会按照如下思路，即先保证大体上相似，再保证局部相似，再保证细节相似，再保证更细微的地方相似……不断地细化下去，无穷次细化以后，仿造的东西将无限接近真品。真假难辨。 物理学家得出结论：把生活中关于“仿造”的经验运用到运动学问题中，如果想仿造一段曲线，那么首先应该保证曲线的起始点一样，其次保证起始点处位移随时间的变化率一样（速度相同），再次应该保证前两者相等的同时关于时间的二阶变化率一样（加速度相同）……如果随时间每一阶变化率（每一阶导数）都一样，那这俩曲线肯定是完全等价的。 所以如果泰勒想一个办法让自己避免接触 sin(x)这类函数，即把这类函数替换掉。 就可以根据这类函数的图像，仿造一个图像，与原来的图像相类似，这种行为在数学上叫近似。不扯这个名词。讲讲如何仿造图像。 仿造的第一步，就是让仿造的曲线也过这个点。 完成了仿造的第一步，很粗糙，甚至完全看不出来这俩有什么相似的地方，那就继续细节化。开始考虑曲线的变化趋势，即导数，保证在此处的导数相等。 经历了第二步，现在起始点相同了，整体变化趋势相近了，可能看起来有那么点意思了。想进一步精确化，应该考虑凹凸性。高中学过：表征图像的凹凸性的参数为“导数的导数”。所以，下一步就让二者的导数的导数相等。 起始点相同，增减性相同，凹凸性相同后，仿造的函数更像了。如果再继续细化下去，应该会无限接近。所以泰勒认为“仿造一段曲线，要先保证起点相同，再保证在此处导数相同，继续保证在此处的导数的导数相同……” 泰勒展开式就是把一个三角函数或者指数函数或者其他比较难缠的函数用多项式替换掉。 也就是说，有一个原函数f(x) ，我再造一个图像与原函数图像相似的多项式函数 g(x) ，为了保证相似，我只需要保证这俩函数在某一点的初始值相等，1阶导数相等，2阶导数相等，……n阶导数相等。 2.2 牛顿法牛顿法的基本思路是，在现有极小点估计值的附近对 f(x) 做二阶泰勒展开，进而找到极小点的下一个估计值。其核心思想是利用迭代点 x_k 处的一阶导数(梯度)和二阶导数(Hessen矩阵)对目标函数进行二次函数近似，然后把二次模型的极小点作为新的迭代点，并不断重复这一过程，直至求得满足精度的近似极小值。 梯度下降算法是将函数在 xn 位置进行一次函数近似，也就是一条直线。计算梯度，从而决定下一步优化的方向是梯度的反方向。而牛顿法是将函数在 xn 位置进行二阶函数近似，也就是二次曲线。计算梯度和二阶导数，从而决定下一步的优化方向。 我们要优化的都是多元函数，x往往不是一个实数，而是一个向量。所以f(x) 的一阶导数也是一个向量，再求导的二阶导数是一个矩阵，就是Hessian矩阵。 2.2.1 泰勒一阶展开牛顿法求根可以按照泰勒一阶展开。例如对于方程 f(x) = 0，我们在任意一点 x0 处，进行一阶泰勒展开： \[f(x) = f(x_0) + (x - x_0)f^{ '}(x_0) \]令 f(x) = 0，带入上式，即可得到： \[x = x_0 - \frac{f(x_0)}{f^{'}(x_0)} \]注意，这里使用了近似，得到的 x 并不是方程的根，只是近似解。但是，x 比 x0 更接近于方程的根。 然后，利用迭代方法求解，以 x 为 x0，求解下一个距离方程的根更近的位置。迭代公式可以写成： \[x_{n+1} = x_n - \frac{f(x_n)}{f^{'}(x_n)} \]2.2.2 泰勒二阶展开机器学习、深度学习中，损失函数的优化问题一般是基于一阶导数梯度下降的。现在，从另一个角度来看，想要让损失函数最小化，这其实是一个最值问题，对应函数的一阶导数 f'(x) = 0。也就是说，如果我们找到了能让 f'(x) = 0 的点 x，损失函数取得最小值，也就实现了模型优化目标。 现在的目标是计算 f’(x) = 0 对应的 x，即 f’(x) = 0 的根。转化为求根问题，就可以利用上一节的牛顿法了。 与上一节有所不同，首先，对 f(x) 在 x0 处进行二阶泰勒展开： \[f(x) = f(x_0) + (x - x_0)f^{'}(x_0) + \frac{1}{2}(x-x_0)^2f^{''}(x_0) \]上式成立的条件是 f(x) 近似等于 f(x0)。令 f(x) = f(x0)，并对 (x - x0) 求导，可得： \[x = x_0 - \frac{f^{'}(x_0)}{f^{''}(x_0)} \]同样，虽然 x 并不是最优解点，但是 x 比 x0 更接近 f'(x) = 0 的根。这样，就能得到最优化的迭代公式： \[x_{n+1} = x_n - \frac{f^{'}(x_n)}{f^{''}(x_n)} \]通过迭代公式，就能不断地找到 f'(x) = 0 的近似根，从而也就实现了损失函数最小化的优化目标。 2.2.3 高维空间在机器学习的优化问题中，我们要优化的都是多元函数，所以x往往不是一个实数，而是一个向量。所以将牛顿求根法利用到机器学习中时,x是一个向量，f(x) 的一阶导数也是一个向量，再求导的二阶导数是一个矩阵，就是Hessian矩阵。 在高维空间，我们用梯度替代一阶导数，用Hessian矩阵替代二阶导数，牛顿法的迭代公式不变： \[x_{k+1} = x_k - [Hf(x_k)^{-1}].J_f(x_k) \]其中 J 定义为 雅克比矩阵，对应一阶偏导数。 推导如下 ： 我们假设f(x)是二阶可微实函数，把f(x)在xk处Taylor展开并取二阶近似为 \[\begin{aligned}&f(x)≈f(x^k)+∇f(x^k)T(x−x^k)+\frac{1}{2}(x−x^k)^T∇^2f(x^k)(x−x^k)\\&x^k为当前的极小值估计值\\&∇f(x^k)是f(x)在x^k处的一阶导数\\&∇^2f(x) 是f(x)在x^k处的Hessen矩阵。\\\end{aligned} \]我们的目标是求f(x)的最小值，而导数为0的点极有可能为极值点，故在此对f(x)求导，并令其导数为0，即∇f(x)=0，可得 \[∇f(x)=∇f(x^k)+∇^2f(x^k)(x−x^k)=0 \]设∇2f(x)可逆，由(2)可以得到牛顿法的迭代公式 \[\begin{aligned}&x^{k+1}=x^k−∇^2f(x^k)^{−1}∇f(x^k) \\&d= −∇^2f(x^k)^{−1}∇f(x^k)被称为牛顿方向 \\\end{aligned} \]当原函数存在一阶二阶连续可导时，可以采用牛顿法进行一维搜索，收敛速度快，具有局部二阶收敛速度。 2.2.4 牛顿法基本流程总结(模仿)一下使用牛顿法求根的步骤: a.已知函数的情况下,随机产生点. b.由已知的点按照的公式进行k次迭代. c.如果与上一次的迭代结果相同或者相差结果小于一个阈值时,本次结果就是函数的根. 伪代码如下 def newton(feature, label, iterMax, sigma, delta): '''牛顿法 input: feature(mat):特征 label(mat):标签 iterMax(int):最大迭代次数 sigma(float), delta(float):牛顿法中的参数 output: w(mat):回归系数 ''' n = np.shape(feature)[1] w = np.mat(np.zeros((n, 1))) it = 0 while it m ? k - m : 0; int l = k = 0; i--) { int j = (i + delta) % m; double dot = sK[j].dot(yK[j]); if (Math.abs(dot) > 0.0) { double rhoJ = 1.0 / dot; alpha[i] = rhoJ * (sK[j].dot(dir.f0)); // 计算alpha dir.f0.plusScaleEqual(yK[j], -alpha[i]); // 重新修正d } } for (int i = 0; i < l; i++) { int j = (i + delta) % m; double dot = sK[j].dot(yK[j]); if (Math.abs(dot) > 0.0) { double rhoJ = 1.0 / dot; double betaI = rhoJ * (yK[j].dot(dir.f0)); // 乘以rho dir.f0.plusScaleEqual(sK[j], (alpha[i] - betaI));// 重新修正d } } } //最后是存储在 OptimVariable.dir 3.5 CalcLosses根据更新的 dir 和 当前系数 计算损失函数误差值，这个损失是为后续的线性搜索准备的。目的是如果损失函数误差值达到允许的范围，那么停止迭代，否则重复迭代。 CalcLosses基本逻辑如下： 1)得到本次步长 Double beta = dir.f1[1] / numSearchStep; 后续UpdateModel 中会对下一次计算的步长（learning rate）进行更新，比如 dir.f1[1] *= 1.0 / (numSearchStep * numSearchStep); 或者 dir.f1[1] = Math.min(dir.f1[1], numSearchStep); 2)调用目标函数的 calcSearchValues 来计算当前系数对应的损失； 3)calcSearchValues 遍历输入labelVectors，对于每个 labelVector 按照线性搜索的步骤进行计算损失。vec[i] += weight * this.unaryLossFunc.loss(etaCoef - i * etaDelta, labelVector.f1); 循环内部如下： 3.1)用x-vec和coef计算出来的 Y ，etaCoef = getEta(labelVector, coefVector); 3.2)以x-vec和dirVec计算出来的 deltaY，etaDelta = getEta(labelVector, dirVec) * beta; 3.3)按照线性搜索的步骤进行计算损失。vec[i] += weight * this.unaryLossFunc.loss(etaCoef - i * etaDelta, labelVector.f1); 联系损失函数可知，etaCoef - i * etaDelta, labelVector.f1 是 训练数据预测值 与 实际类别 的偏差； 4)为后续聚合 lossAllReduce 准备数据；代码如下： public class CalcLosses extends ComputeFunction { @Override public void calc(ComContext context) { Iterable labledVectors = context.getObj(OptimVariable.trainData); Tuple2 dir = context.getObj(OptimVariable.dir); Tuple2 coef = context.getObj(OptimVariable.currentCoef); if (objFunc == null) { objFunc = ((List)context.getObj(OptimVariable.objFunc)).get(0); } /** * calculate losses of current coefficient. * if optimizer is owlqn, constraint search will used, else without constraint. */ Double beta = dir.f1[1] / numSearchStep; double[] vec = method.equals(OptimMethod.OWLQN) ? objFunc.constraintCalcSearchValues(labledVectors, coef.f0, dir.f0, beta, numSearchStep) : objFunc.calcSearchValues(labledVectors, coef.f0, dir.f0, beta, numSearchStep); // 变量为 dir = {Tuple2@9988} "(-9.599 -9.599 -3.5515274823133662 -3.5911942493220335,[4.0, 0.1])" coef = {Tuple2@9989} "(0.001 0.0 0.0 0.0,1.7976931348623157E308)" beta = {Double@10014} 0.025 vec = {double[5]@10015} 0 = 0.0 1 = 0.0 2 = 0.0 3 = 0.0 4 = 0.0 // prepare buffer vec for allReduce. double[] buffer = context.getObj(OptimVariable.lossAllReduce); if (buffer == null) { buffer = vec.clone(); context.putObj(OptimVariable.lossAllReduce, buffer); } else { System.arraycopy(vec, 0, buffer, 0, vec.length); } } }其中搜索是一元目标函数提供的，其又调用了损失函数。 public class UnaryLossObjFunc extends OptimObjFunc { /** * Calculate loss values for line search in optimization. * * @param labelVectors train data. * @param coefVector coefficient of current time. * @param dirVec descend direction of optimization problem. * @param beta step length of line search. * @param numStep num of line search step. * @return double[] losses. */ @Override public double[] calcSearchValues(Iterable labelVectors, DenseVector coefVector, DenseVector dirVec, double beta, int numStep) { double[] vec = new double[numStep + 1]; // labelVector是三元组Tuple3 labelVectors = {ArrayList@10007} size = 4 0 = {Tuple3@10027} "(1.0,16.8,1.0 1.0 1.4657097546055162 1.4770978917519928)" 1 = {Tuple3@10034} "(1.0,6.7,1.0 1.0 -0.338240712601273 -0.7385489458759964)" 2 = {Tuple3@10035} "(1.0,6.9,1.0 1.0 -0.7892283294029703 -0.3692744729379982)" 3 = {Tuple3@10036} "(1.0,8.0,1.0 1.0 -0.338240712601273 -0.3692744729379982)" coefVector = {DenseVector@10008} "0.001 0.0 0.0 0.0" dirVec = {DenseVector@10009} "-9.599 -9.599 -3.5515274823133662 -3.5911942493220335" beta = 0.025 numStep = 4 vec = {double[5]@10026} 0 = 0.0 1 = 0.0 2 = 0.0 3 = 0.0 4 = 0.0 for (Tuple3 labelVector : labelVectors) { double weight = labelVector.f0; //用x-vec和coef计算出来的 Y double etaCoef = getEta(labelVector, coefVector); //以x-vec和dirVec计算出来的 deltaY double etaDelta = getEta(labelVector, dirVec) * beta; weight = 1.0 etaCoef = 0.001 etaDelta = -0.7427013482280431 for (int i = 0; i < numStep + 1; ++i) { //labelVector.f1就是label y //联系下面损失函数可知，etaCoef - i * etaDelta, labelVector.f1 是 训练数据预测值 与 实际类别 的偏差 vec[i] += weight * this.unaryLossFunc.loss(etaCoef - i * etaDelta, labelVector.f1); } } return vec; } private double getEta(Tuple3 labelVector, DenseVector coefVector) { //labelVector.f2 = {DenseVector@9972} "1.0 1.0 1.4657097546055162 1.4770978917519928" return MatVecOp.dot(labelVector.f2, coefVector); } } vec = {double[5]@10026} 0 = 219.33160199999998 1 = 198.85962729259512 2 = 179.40202828917856 3 = 160.95880498975038 4 = 143.52995739431051回顾损失函数如下 /** * Squared loss function. * https://en.wikipedia.org/wiki/Loss_functions_for_classification#Square_loss */ public class SquareLossFunc implements UnaryLossFunc { public SquareLossFunc() { } @Override public double loss(double eta, double y) { return 0.5 * (eta - y) * (eta - y); } @Override public double derivative(double eta, double y) { return eta - y; } @Override public double secondDerivative(double eta, double y) { return 1; } } 3.6 UpdateModel本模块做两件事 基于dir和step length来更新coefficient，即依据方向和步长计算。 如果简化理解，参数更新公式为 ：下一刻参数 = 上一时刻参数 - 学习率 * (损失函数对这个参数的导数)。 判断循环的收敛。因为变量太多，所以有时候就忘记谁是谁了，所以再次标示。 OptimVariable.dir是CalcGradient计算出来的梯度做修正之后的结果 OptimVariable.lossAllReduce 这个会变化，此时是上一阶段计算的损失代码逻辑大致如下： 1)得出"最新损失"的最小值位置pos 2)得出学习率 beta = dir.f1[1] / numSearchStep; 3)根据"最新损失pos"和 上一个最小值 last loss value判断来进行分别处理： 3.1)如果所有损失都比 last loss value 大，则 3.1.1)缩减学习率 multiply 1.0 / (numSearchStep*numSearchStep) 3.1.2)把 eta 设置为 0；这个就是步长了 3.1.3)把当前 loss 设置为 last loss value 3.2)如果losses[numSearchStep] 比 last loss value 小，则 3.2.1)增大学习率 multiply numSearchStep 3.2.2)设置eta是smallest value pos，eta = beta * pos; 这个eta就是步长了 3.2.3)把当前 loss 设置为当前loss最小值 losses[numSearchStep] 3.3)否则 3.3.1)学习率不更改 3.3.2)设置eta是smallest value pos，eta = beta * pos; 这个eta就是步长了 3.3.3)把当前 loss 设置为当前loss最小值 losses[pos] 4)修正Hessian矩阵 5)用方向向量和步长来来更新系数向量 curCoef.f0.plusScaleEqual(dir.f0, -eta); 6)如果当前loss比 min loss 小，则 用 current loss 更新 min loss 6.1)minCoef.f1 = curCoef.f1; 更新最小loss 6.2)minCoef.f0.setEqual(curCoef.f0); 更新最小loss所对应的Coef，即线性模型最后需要的系数在这里求步长，我没有发现 Wolf-Powell 准则的使用，Alink做了某种优化。如果有朋友能看出来Wolf-Powell如何使用，还请留言指点 ，谢谢。 代码如下： public class UpdateModel extends ComputeFunction { @Override public void calc(ComContext context) { double[] losses = context.getObj(OptimVariable.lossAllReduce); Tuple2 dir = context.getObj(OptimVariable.dir); Tuple2 pseGrad = context.getObj(OptimVariable.pseGrad); Tuple2 curCoef = context.getObj(OptimVariable.currentCoef); Tuple2 minCoef = context.getObj(OptimVariable.minCoef); double lossChangeRatio = 1.0; double[] lossCurve = context.getObj(OptimVariable.lossCurve); for (int j = 0; j < losses.length; ++j) { losses[j] /= dir.f1[0]; //dir.f1[0]是权重 } int pos = -1; //get the min value of losses, and remember the position. for (int j = 0; j < losses.length; ++j) { if (losses[j] < losses[0]) { losses[0] = losses[j]; pos = j; } } // adaptive learningRate strategy double beta = dir.f1[1] / numSearchStep; double eta; // 变量如下 losses = {double[5]@10001} 0 = 35.88248934857763 1 = 49.71490682314878 2 = 44.85050707229464 3 = 40.239701247437594 4 = 35.88248934857763 dir = {Tuple2@10002} "(-9.599 -9.599 -3.5515274823133662 -3.5911942493220335,[4.0, 0.1])" pseGrad = null curCoef = {Tuple2@10003} "(0.001 0.0 0.0 0.0,1.7976931348623157E308)" minCoef = {Tuple2@10004} "(0.001 0.0 0.0 0.0,1.7976931348623157E308)" lossChangeRatio = 1.0 lossCurve = {double[100]@10005} pos = 4 beta = 0.025 curCoef.f1 = {Double@10006} 1.7976931348623157E308 if (pos == -1) { /** * if all losses larger than last loss value. we'll do the below things: * 1. reduce learning rate by multiply 1.0 / (numSearchStep*numSearchStep). * 2. set eta with zero. * 3. set current loss equals last loss value. */ eta = 0; dir.f1[1] *= 1.0 / (numSearchStep * numSearchStep); // 学习率 curCoef.f1 = losses[0]; // 最小loss } else if (pos == numSearchStep) { /** * if losses[numSearchStep] smaller than last loss value. we'll do the below things: * 1. enlarge learning rate by multiply numSearchStep. * 2. set eta with the smallest value pos. * 3. set current loss equals smallest loss value. */ eta = beta * pos; dir.f1[1] *= numSearchStep; dir.f1[1] = Math.min(dir.f1[1], numSearchStep); // 学习率 lossChangeRatio = Math.abs((curCoef.f1 - losses[pos]) / curCoef.f1); curCoef.f1 = losses[numSearchStep]; // 最小loss // 当前数值 numSearchStep = 4 是NUM_SEARCH_STEP缺省值，是线性搜索的参数，Line search parameter, which define the search value num of one step. lossChangeRatio = 1.0 pos = 4 eta = 0.1 curCoef.f1 = {Double@10049} 35.88248934857763 dir.f1[1] = 0.4 } else { /** * else : * 1. learning rate not changed. * 2. set eta with the smallest value pos. * 3. set current loss equals smallest loss value. */ eta = beta * pos; lossChangeRatio = Math.abs((curCoef.f1 - losses[pos]) / curCoef.f1); curCoef.f1 = losses[pos]; // 最小loss } /* update loss value in loss curve at this step */ lossCurve[context.getStepNo() - 1] = curCoef.f1; lossCurve = {double[100]@9998} 0 = 35.88248934857763 1 = Infinity if (method.equals(OptimMethod.OWLQN)) { ..... } else if (method.equals(OptimMethod.LBFGS)) { Tuple2 sKyK = context.getObj(OptimVariable.sKyK); int size = dir.f0.size(); int k = context.getStepNo() - 1; DenseVector[] sK = sKyK.f0; for (int s = 0; s < size; ++s) { // 修正矩阵 sK[k % OptimVariable.numCorrections].set(s, dir.f0.get(s) * (-eta)); } curCoef.f0.plusScaleEqual(dir.f0, -eta); // 这里是用方向向量和步长来更新系数向量 sKyK = {Tuple2@10043} "([0.9599000000000001 0.9599000000000001 0.35515274823133663 0.3591194249322034, 0.0 0.0 0.0 0.0, 0.0 0.0 0.0 0.0, 0.0 0.0 0.0 0.0, 0.0 0.0 0.0 0.0, 0.0 0.0 0.0 0.0, 0.0 0.0 0.0 0.0, 0.0 0.0 0.0 0.0, 0.0 0.0 0.0 0.0, 0.0 0.0 0.0 0.0],[0.0 0.0 0.0 0.0, 0.0 0.0 0.0 0.0, 0.0 0.0 0.0 0.0, 0.0 0.0 0.0 0.0, 0.0 0.0 0.0 0.0, 0.0 0.0 0.0 0.0, 0.0 0.0 0.0 0.0, 0.0 0.0 0.0 0.0, 0.0 0.0 0.0 0.0, 0.0 0.0 0.0 0.0])" f0 = {DenseVector[10]@10044} 0 = {DenseVector@10074} "0.9599000000000001 0.9599000000000001 0.35515274823133663 0.3591194249322034" } /** * if current loss is smaller than min loss, then update the min loss and min coefficient by current. */ if (curCoef.f1 < minCoef.f1) { minCoef.f1 = curCoef.f1; // 最小loss minCoef.f0.setEqual(curCoef.f0); // 最小loss所对应的Coef，即线性模型最后需要的系数 } curCoef = {Tuple2@9996} "(0.9609000000000001 0.9599000000000001 0.35515274823133663 0.3591194249322034,35.88248934857763)" f0 = {DenseVector@10059} "0.9609000000000001 0.9599000000000001 0.35515274823133663 0.3591194249322034" f1 = {Double@10048} 35.88248934857763 minCoef = {Tuple2@9997} "(0.9609000000000001 0.9599000000000001 0.35515274823133663 0.3591194249322034,35.88248934857763)" f0 = {DenseVector@10059} "0.9609000000000001 0.9599000000000001 0.35515274823133663 0.3591194249322034" f1 = {Double@10048} 35.88248934857763 // judge the convergence of iteration. filter(dir, curCoef, minCoef, context, lossChangeRatio); } }filter 判断是否收敛，里面的打印log很清晰的说明了函数逻辑。 /** * judge the convergence of iteration. */ public void filter(Tuple2 grad, Tuple2 c, Tuple2 m, ComContext context, double lossChangeRatio) { double epsilon = params.get(HasEpsilonDv0000001.EPSILON); int maxIter = params.get(HasMaxIterDefaultAs100.MAX_ITER); double gradNorm = ((Tuple2)context.getObj(gradName)).f0.normL2(); if (c.f1 < epsilon || gradNorm < epsilon) { printLog(" method converged at step : ", c.f1, m.f1, grad.f1[1], gradNorm, context, lossChangeRatio); grad.f1[0] = -1.0; } else if (context.getStepNo() > maxIter - 1) { printLog(" method stop at max step : ", c.f1, m.f1, grad.f1[1], gradNorm, context, lossChangeRatio); grad.f1[0] = -1.0; } else if (grad.f1[1] < EPS) { printLog(" learning rate is too small, method stops at step : ", c.f1, m.f1, grad.f1[1], gradNorm, context, lossChangeRatio); grad.f1[0] = -1.0; } else if (lossChangeRatio < epsilon && gradNorm < Math.sqrt(epsilon)) { printLog(" loss change ratio is too small, method stops at step : ", c.f1, m.f1, grad.f1[1], gradNorm, context, lossChangeRatio); grad.f1[0] = -1.0; } else { printLog(" method continue at step : ", c.f1, m.f1, grad.f1[1], gradNorm, context, lossChangeRatio); } } 3.7 OutputModel.closeWith(new OutputModel()) 这部分是每次迭代结束，临时输出模型，把数据转换成Flink通用的Row类型，Transfer the state to model rows。 public class OutputModel extends CompleteResultFunction { @Override public List calc(ComContext context) { // get the coefficient of min loss. Tuple2 minCoef = context.getObj(OptimVariable.minCoef); double[] lossCurve = context.getObj(OptimVariable.lossCurve); int effectiveSize = lossCurve.length; for (int i = 0; i < lossCurve.length; ++i) { if (Double.isInfinite(lossCurve[i])) { effectiveSize = i; break; } } double[] effectiveCurve = new double[effectiveSize]; System.arraycopy(lossCurve, 0, effectiveCurve, 0, effectiveSize); Params params = new Params(); params.set(ModelParamName.COEF, minCoef.f0);// 重点在这里，minCoef是我们真正需要的 params.set(ModelParamName.LOSS_CURVE, effectiveCurve); List model = new ArrayList (1); model.add(Row.of(params.toJson())); return model; } } 0x04 准备模型元数据这里设置了并行度为1。 // Prepare the meta info of linear model. DataSet meta = labelInfo.f0 .mapPartition(new CreateMeta(modelName, linearModelType, isRegProc, params)) .setParallelism(1);具体代码 public static class CreateMeta implements MapPartitionFunction { @Override public void mapPartition(Iterable rows, Collector metas) throws Exception { Object[] labels = null; if (!this.isRegProc) { labels = orderLabels(rows); } Params meta = new Params(); meta.set(ModelParamName.MODEL_NAME, this.modelName); meta.set(ModelParamName.LINEAR_MODEL_TYPE, this.modelType); meta.set(ModelParamName.LABEL_VALUES, labels); meta.set(ModelParamName.HAS_INTERCEPT_ITEM, this.hasInterceptItem); meta.set(ModelParamName.VECTOR_COL_NAME, vectorColName); meta.set(LinearTrainParams.LABEL_COL, labelName); metas.collect(meta); } } // 变量为 meta = {Params@9667} "Params {hasInterceptItem=true, vectorColName=null, modelName="Linear Regression", labelValues=null, labelCol="label", linearModelType="LinearReg"}" 0x05 建立模型当迭代循环结束之后，Alink就根据Coef数据来建立模型。 /** * build the linear model rows, the format to be output. */ public static class BuildModelFromCoefs extends AbstractRichFunction implements @Override public void mapPartition(Iterable iterable, Collector collector) throws Exception { for (Tuple2 coefVector : iterable) { LinearModelData modelData = buildLinearModelData(meta, featureNames, labelType, meanVar, hasIntercept, standardization, coefVector); new LinearModelDataConverter(this.labelType).save(modelData, collector); } } }得到模型数据为，里面coef就是 f(x)=w^Tx+b 中的 w, b。是最终用来计算的。 modelData = {LinearModelData@10584} featureNames = {String[3]@9787} featureTypes = null vectorColName = null coefVector = {DenseVector@10485} "-3.938937407856857 4.799499941426075 0.8929571907809862 1.078169576770847" coefVectors = null vectorSize = 3 modelName = "Linear Regression" labelName = null labelValues = null linearModelType = {LinearModelType@4674} "LinearReg" hasInterceptItem = true lossCurve = {double[12]@10593} 0 = 35.88248934857763 1 = 12.807366842002144 2 = 0.5228366663917704 3 = 0.031112070740366038 4 = 0.01098914933042993 5 = 0.009765757443537283 6 = 0.008750523231785415 7 = 0.004210085397869248 8 = 0.0039042232755530704 9 = 0.0038821509860327537 10 = 0.003882042680010676 11 = 0.0038820422536391033 labelType = {FractionalTypeInfo@9790} "Double" 0x06 使用模型预测预测时候，使用的是LinearModelMapper，其内部部分变量打印出来如下，能够看出来模型数据。 this = {LinearModelMapper@10704} vectorColIndex = -1 model = {LinearModelData@10597} featureNames = {String[3]@10632} featureTypes = null vectorColName = null coefVector = {DenseVector@10612} "-3.938937407856857 4.799499941426075 0.8929571907809862 1.078169576770847" coefVectors = null vectorSize = 0 modelName = "Linear Regression" labelName = null labelValues = {Object[0]@10613} linearModelType = {LinearModelType@10611} "LinearReg" hasInterceptItem = true lossCurve = null labelType = null具体预测是在 LinearModelMapper.predict 中完成，具体如下： 对应原始输入 Row.of("$3$0:1.0 1:7.0 2:9.0", "1.0 7.0 9.0", 1.0, 7.0, 9.0, 16.8), ， 通过 FeatureLabelUtil.getFeatureVector 处理之后， 得到的四元组是 "1.0 1.0 7.0 9.0"，其中第一个 1.0 是通过 aVector.set(0, 1.0); 专门设定的固定值。比如模型是 f(x) = ax + b，这个固定值 1.0 就是 b 的初始化值，随着优化过程会得到 b。所以这里还是需要有一个 1.0 来进行预测。 模型系数是："-3.938937407856857 4.799499941426075 0.8929571907809862 1.078169576770847" 。 四元组 和 模型系数 点积的结果就是 dotValue = 16.814789059973744。这样就能看出来模型系数如何使用的了。 public Object predict(Vector vector) throws Exception { double dotValue = MatVecOp.dot(vector, model.coefVector); switch (model.linearModelType) { case LR: case SVM: case Perceptron: return dotValue >= 0 ? model.labelValues[0] : model.labelValues[1]; case LinearReg: case SVR: return dotValue; } } vector = {DenseVector@10610} "1.0 1.0 7.0 9.0" data = {double[4]@10619} 0 = 1.0 1 = 1.0 2 = 7.0 3 = 9.0 model.coefVector = {DenseVector@10612} "-3.938937407856857 4.799499941426075 0.8929571907809862 1.078169576770847" data = {double[4]@10620} 0 = -3.938937407856857 1 = 4.799499941426075 2 = 0.8929571907809862 3 = 1.078169576770847 0x07 本系列其他文章Alink漫谈(一) : 从KMeans算法实现不同看Alink设计思想 Alink漫谈(二) : 从源码看机器学习平台Alink设计和架构 [Alink漫谈之三] AllReduce通信模型 Alink漫谈(四) : 模型的来龙去脉 Alink漫谈(五) : 迭代计算和Superstep Alink漫谈(六) : TF-IDF算法的实现 Alink漫谈(七) : 如何划分训练数据集和测试数据集 Alink漫谈(八) : 二分类评估 AUC、K-S、PRC、Precision、Recall、LiftChart 如何实现 Alink漫谈(九) ：特征工程 之 特征哈希/标准化缩放 Alink漫谈(十) ：线性回归实现 之 数据预处理 0xFF 参考http://www.mamicode.com/info-detail-2508527.html) 机器学习算法实现解析——liblbfgs之L-BFGS算法 深入机器学习系列之BFGS & L-BFGS 拟牛顿法公式推导以及python代码实现（一） 浅显易懂——泰勒展开式 泰勒展开 — Taylor Expansion 从牛顿插值法到泰勒公式 怎样更好地理解并记忆泰勒展开式？ 优化算法——牛顿法(Newton Method) 优化算法——拟牛顿法之L-BFGS算法 一文读懂L-BFGS算法 《分布式机器学习算法、理论与实践_刘铁岩》 LogisticRegressionWithLBFGS--逻辑回归 逻辑回归优化方法-L-BFGS 机器学习与运筹优化（三）从牛顿法到L-BFGS 机器学习中牛顿法凸优化的通俗解释 深入机器学习系列17-BFGS & L-BFGS 优化方法基础系列-精确的一维搜索技术 优化方法基础系列-非精确的一维搜索技术 [原创]用“人话”解释不精确线搜索中的Armijo-Goldstein准则及Wolfe-Powell准则 https://www.zhihu.com/question/36425542 一维搜索算法介绍及其实现 数值优化|笔记整理（1）——引入，线搜索：步长选取条件 https://zhuanlan.zhihu.com/p/32821110 线搜索（一）：步长的选取 最优化问题——一维搜索(二) CRF L-BFGS Line Search原理及代码分析 步长与学习率 机器学习优化算法L-BFGS及其分布式实现 L-BFGS算法详解（逻辑回归的默认优化算法） 线性回归的原理及实践（牛顿法） 线性回归、梯度下降（Linear Regression、Gradient Descent） L-BFGS算法 并行逻辑回归 |
| CopyRight 2018-2019 实验室设备网 版权所有 |