最优化算法 - NNLS(非负最小二乘)

优质
小牛编辑
135浏览
2023-12-01

  spark中的非负正则化最小二乘法并不是wiki中介绍的NNLS的实现,而是做了相应的优化。它使用改进投影梯度法结合共轭梯度法来求解非负最小二乘。
在介绍spark的源码之前,我们要先了解什么是最小二乘法以及共轭梯度法。

1 最小二乘法

1.1 最小二乘问题

  在某些最优化问题中,目标函数由若干个函数的平方和构成,它的一般形式如下所示:

1.1

  其中x=(x1,x2,…,xn),一般假设m>=n。把极小化这类函数的问题称为最小二乘问题。

1.2

  当$f{i}(x)$为x的线性函数时,称(1.2)为线性最小二乘问题,当$f{i}(x)$为x的非线性函数时,称(1.2)为非线性最小二乘问题。

1.2 线性最小二乘问题

  在公式(1.1)中,假设

1.3

  其中,p是n维列向量,bi是实数,这样我们可以用矩阵的形式表示(1.1)式。令

1.3

  A是m * n矩阵,bm维列向量。则

1.4

  因为F(x)是凸的,所以对(1.4)求导可以得到全局极小值,令其导数为0,我们可以得到这个极小值。

1.5

  假设A为满秩,$A^{T}A$为n阶对称正定矩阵,我们可以求得x的值为以下的形式:

1.6

1.3 非线性最小二乘问题

  假设在(1.1)中,$f{i}(x)$为非线性函数,且F(x)有连续偏导数。由于$f{i}(x)$为非线性函数,所以(1.2)中的非线性最小二乘无法套用(1.6)中的公式求得。
解这类问题的基本思想是,通过解一系列线性最小二乘问题求非线性最小二乘问题的解。设$x^{(k)}$是解的第k次近似。在$x^{(k)}$时,将函数$f_{i}(x)$线性化,从而将非线性最小二乘转换为线性最小二乘问题,
用(1.6)中的公式求解极小点$x^{(k+1)}$ ,把它作为非线性最小二乘问题解的第k+1次近似。然后再从$x^{(k+1)}$出发,继续迭代。下面将来推导迭代公式。令

1.7

  上式右端是$f_{i}(x)$在$x^{(k)}$处展开的一阶泰勒级数多项式。令

1.8

  用∅(x)近似F(x),从而用∅(x)的极小点作为目标函数F(x)的极小点的估计。现在求解线性最小二乘问题

1.9

  把(1.9)写成

1.10

  在公式(1.10)中,

1.10

  将$A_{k}$和b带入公式(1.5)中,可以得到,

1.11

  如果$A{k}$为列满秩,且$A{k}^{T}A_{k}$是对称正定矩阵,那么由(1.11)可以得到x的极小值。

1.12

  可以推导出$2A{k}^{T}f{k}$是目标函数F(x)在$x^{(k)}$处的梯度,$2A{k}^{T}A{k}$是函数∅(x)的海森矩阵。所以(1.12)又可以写为如下形式。

1.13

  公式(1.13)称为Gauss-Newton公式。向量

1.14

  称为点$x^{(k)}$处的Gauss-Newton方向。为保证每次迭代能使目标函数值下降(至少不能上升),在求出$d^{(k)}$后,不直接使用$x^{(k)}+d^{(k)}$作为k+1次近似,而是从$x^{(k)}$出发,沿$d^{(k)}$方向进行一维搜索。

1.15

  求出步长$\lambda ^{(k)}$后,令

1.16

  最小二乘的计算步骤如下:

  • (1) 给定初始点$x^{(1)}$,允许误差ε>0,k=1

  • (2) 计算函数值$f{i}(x)$,得到向量$f^{(k)}$,再计算一阶偏导,得到m*n矩阵$A{(k)}$

  • (3) 解方程组(1.14)求得Gauss-Newton方向$d^{(k)}$

  • (4)从$x^{(k)}$出发,沿着$d^{(k)}$作一维搜索,求出步长$\lambda ^{(k)}$ ,并令$x^{(k+1)}=x^{(k)}-\lambda d^{(k)}$

  • (5)若$||x^{(k+1)}-x^{(k)}||<=\epsilon$停止迭代,求出x,否则,k=k+1,返回步骤(2)

  在某些情况下,矩阵$A^{T}A$是奇异的,这种情况下,我们无法求出它的逆矩阵,因此我们需要对其进行修改。用到的基本技巧是将一个正定对角矩阵添加到$A^{T}A$上,改变原来矩阵的特征值结构,使其变成条件较好的对称正定矩阵。
典型的算法是Marquardt

1.17

  其中,In阶单位矩阵,alpha是一个正实数。当alpha为0时,$d^{(k)}$就是Gauss-Newton方向,当alpha充分大时,这时$d^{(k)}$接近F(x)在$x^{(k)}$处的最速下降方向。算法的具体过程见参考文献【1】。

2 共轭梯度法

2.1 共轭方向

  在讲解共轭梯度法之前,我们需要先知道什么是共轭方向,下面的定义给出了答案。

  定义2.1An*n对称正定矩阵,若两个方向$d^{(1)}$ 和$d^{(2)}$满足

2.1

则称这两个方向关于A共轭。若$d^{(1)},d^{(2)},…,d^{(k)}$是k个方向,它们两两关于A共轭,则称这组方向是关于A共轭的。即

2.2

在上述定义中,如果A是单位矩阵,那么两个方向关于A共轭等价于两个方向正交。如果A是一般的对称正定矩阵,$d^{(i)}$与$d^{(j)}$共轭,就是$d^{(i)}$与$Ad^{(j)}$正交。共轭方向有一些重要的性质。

  定理2.1An阶对称正定矩阵,$d^{(1)},d^{(2)},…,d^{(k)}$是kA的共轭的非零向量,则这个向量组线性无关。

  定理2.2 (扩张子空间定理) 设有函数

2.2.a1

  其中,An阶对称正定矩阵,$d^{(1)},d^{(2)},…,d^{(k)}$是kA的共轭的非零向量,以任意的$x^{(1)}$为初始点,
沿$d^{(1)},d^{(2)},…,d^{(k)}$进行一维搜索,得到$x^{(2)},x^{(3)},…,x^{(k+1)}$,则$d^{(k+1)}$是线性流型
$x^{(1)}+H_{k}$上的唯一极小点,特别的,当k=n时,$x^{(n+1)}$是函数f(x)的唯一极小点。其中,

2.2.a2

是$d^{(1)},d^{(2)},…,d^{(k)}$生成的子空间。

  这两个定理在文献【1】中有详细证明。

2.2 共轭梯度法

  共轭梯度法的基本思想是把共轭性与最速下降方法相结合,利用已知点处的梯度构造一组共轭方向,并沿这组方向进行搜索,求出目标函数的极小点。这里我们仅仅讨论二次凸函数的共轭梯度法。

  考虑问题

2.3

  其中A是对称正定矩阵,c是常数。

  具体求解方式如下:

  首先给定任何一个初始点$x^{(1)}$ ,计算目标函数f(x)在这点的梯度$g{(1)}$ ,若$||g{(1)}||=0$ ,则停止计算;否则令

2.4

  沿方向$d^{(1)}$搜索,得到点$x^{(2)}$ 。计算$x^{(2)}$处的梯度,若$||g{(2)}||!=0$,则利用$g{(2)}$和$d^{(1)}$构造第二个搜索方向$d^{(2)}$,再沿$d^{(2)}$搜索。

  一般的,若已知$x^{(k)}$和搜索方向$d^{(k)}$,则从$x^{(k))}$出发,沿方向$d^{(k)}$搜索,得到

2.5

  其中步长lambda满足

2.5.a1

  此时可以求得lambda的显式表达。令

2.5.a2

  通过求导可以求上面公式的极小值,即

2.6

  根据二次函数梯度表达式,(2.6)式可以推出如下方程

2.7

  由(2.7)式可以得到

2.8

  计算f(x)在$x^{(k+1)}$处的梯度,若$||g{k+1}||=0$ ,则停止计算,
否则用$g
{(k+1)}$和$d^{(k)}$构造下一个搜索方向$d^{(k+1)}$ ,并使$d^{(k)}$与$d^{(k+1)}$共轭。按照这种设想,令

2.9

  在公式(2.9)两端同时乘以$d^{(k)T}A$,并令

2.10

  可以求得

2.11

  再从$x^{(k+1)}$出发,沿$d^{(k+1)}$方向搜索。综上分析 ,在第1个搜索方向取负梯度的前提下,重复使用公式(2.5)、(2.8)、(2.9)、(2.11),我们就能够构造一组搜索方向。当然,前提是这组方向是关于A共轭的。
定理2.3说明了这组方向是共轭的。

  定理2.3 对于正定二次函数(2.3),具有精确一维搜索的的共轭梯度法在m<=n次一维搜索后终止,并且对于所有i(1<=i<=m),下列关系成立:

2.11.a1

  还可以证明,对于正定二次函数,运用共轭梯度法时,不做矩阵运算也可以计算变量beta_k

  定理2.4 对于正定二次函数,共轭梯度法中因子beta_k具有下列表达式

2.12

  对于二次凸函数,共轭梯度法的计算步骤如下:

NNLS.2.1

3 最小二乘法在spark中的具体实现

  Spark ml中解决最小二乘可以选择两种方式,一种是非负正则化最小二乘,一种是乔里斯基分解(Cholesky)。

  乔里斯基分解分解是把一个对称正定的矩阵表示成一个上三角矩阵U的转置和其本身的乘积的分解。在ml代码中,直接调用netlib-java封装的dppsv方法实现。

  1. lapack.dppsv(“u”, k, 1, ne.ata, ne.atb, k, info)

  可以深入dppsv代码(Fortran代码)了解更深的细节。我们分析的重点是非负正则化最小二乘的实现,因为在某些情况下,方程组的解为负数是没有意义的。虽然方程组可以得到精确解,但却不能取负值解。在这种情况下,其非负最小二乘解比方程的精确解更有意义。非负最小二乘问题要求解的问题如下公式

3.1

其中ata是半正定矩阵。

  在ml代码中,org.apache.spark.mllib.optimization.NNLS对象实现了非负最小二乘算法。该算法结合了投影梯度算法和共轭梯度算法来求解。

  首先介绍NNLS对象中的Workspace类。

  1. class Workspace(val n: Int) {
  2. val scratch = new Array[Double](n)
  3. val grad = new Array[Double](n) //投影梯度
  4. val x = new Array[Double](n)
  5. val dir = new Array[Double](n) //搜索方向
  6. val lastDir = new Array[Double](n)
  7. val res = new Array[Double](n) //梯度
  8. }

  在Workspace中,res表示梯度,grad表示梯度的投影,dir表示迭代过程中的搜索方向(共轭梯度中的搜索方向$d^{(k)}$),scratch代表公式(2.8)中的
$d^{(k)T}A$。

  NNLS对象中,sort方法用来解最小二乘,它通过迭代求解极小值。我们将分步骤剖析该方法。

  • (1)确定迭代次数。
  1. val iterMax = math.max(400, 20 * n)
  • (2)求梯度。

  在每次迭代内部,第一步会求梯度res,代码如下

  1. //y := alpha*A*x + beta*y 即 y:=1.0 * ata * x + 0.0 * res
  2. blas.dgemv("N", n, n, 1.0, ata, n, x, 1, 0.0, res, 1)
  3. // ata*x - atb
  4. blas.daxpy(n, -1.0, atb, 1, res, 1)

  dgemv方法的作用是得到y := alpha*A*x + beta*y,在本代码中表示res=ata*xdaxpy方法的作用是得到y:=step*x +y,在本代码中表示res=ata*x-atb,即梯度。

  • (3)求梯度的投影矩阵

  求梯度矩阵的投影矩阵的依据如下公式。

3.2

  详细代码如下所示:

  1. //转换为投影矩阵
  2. i = 0
  3. while (i < n) {
  4. if (grad(i) > 0.0 && x(i) == 0.0) {
  5. grad(i) = 0.0
  6. }
  7. i = i + 1
  8. }
  • (4)求搜索方向。

  在第一次迭代中,搜索方向即为梯度方向。如下面代码所示。

  1. //在第一次迭代中,搜索方向dir即为梯度方向
  2. blas.dcopy(n, grad, 1, dir, 1)

  在第k次迭代中,搜索方向由梯度方向和前一步的搜索方向共同确定,计算依赖的公式是(2.9)。具体代码有两行

  1. val alpha = ngrad / lastNorm
  2. //alpha*lastDir + dir,此时dir为梯度方向
  3. blas.daxpy(n, alpha, lastDir, 1, dir, 1)

  此处的alpha就是根据公式(2.12)计算的。

  • (5)计算步长。

  知道了搜索方向,我们就可以依据公式(2.8)来计算步长。

  1. def steplen(dir: Array[Double], res: Array[Double]): Double = {
  2. //top = g * d
  3. val top = blas.ddot(n, dir, 1, res, 1)
  4. // y := alpha*A*x + beta*y.
  5. // scratch = d * ata
  6. blas.dgemv("N", n, n, 1.0, ata, n, dir, 1, 0.0, scratch, 1)
  7. //公式(2.8),添加1e-20是为了避免分母为0
  8. //g*d/d*ata*d
  9. top / (blas.ddot(n, scratch, 1, dir, 1) + 1e-20)
  10. }
  • (6)调整步长并修改迭代值。

  因为解是非负的,所以步长需要做一定的处理,如果步长与搜索方向的乘积大于x的值,那么重置步长。重置逻辑如下:

  1. i = 0
  2. while (i < n) {
  3. if (step * dir(i) > x(i)) {
  4. //如果步长过大,那么二者的商替代
  5. step = x(i) / dir(i)
  6. }
  7. i = i + 1
  8. }

  最后,修改x的值,完成该次迭代。

  1. i = 0
  2. while (i < n) {
  3. // x(i)趋向为0
  4. if (step * dir(i) > x(i) * (1 - 1e-14)) {
  5. x(i) = 0
  6. lastWall = iterno
  7. } else {
  8. x(i) -= step * dir(i)
  9. }
  10. i = i + 1
  11. }

4 参考文献

【1】陈宝林,最优化理论和算法

【2】Boris T. Polyak,The conjugate gradient method in extreme problems

【3】https://en.wikipedia.org/wiki/Non-negative_least_squares