UA MATH567 高维统计专题3 含L1-norm的凸优化4 Nesterov方法与Accelerate Proximal Gradient

UA MATH567 高维统计专题3 含L1-norm的凸优化4 一阶方法的加速

Nesterov方法
Accelerate Proximal Gradient (APG)

梯度下降与Proximal gradient descent都是一阶方法，也就是借助目标函数的某种一阶导进行迭代找最小值的方法，假设目标函数f(x)f(x)f(x)是一个平滑的凸函数，则对于min⁡xf(x)\min_x f(x)minxf(x)，一阶方法可以用下面的递推式表示
xk+1=ϕk(x0,x1,⋯,xk,∇f(x0),∇f(x1),⋯,∇f(xk))x_{k+1}=\phi_k(x_0,x_1,\cdots,x_k,\nabla f(x_0),\nabla f(x_1),\cdots,\nabla f(x_k))xk+1=ϕk(x0,x1,⋯,xk,∇f(x0),∇f(x1),⋯,∇f(xk))

但是一阶方法只利用了梯度信息，也就是单调性的信息（kkk次迭代后误差为O(1/k)O(1/k)O(1/k)，而误差的下界为O(1/k2)O(1/k^2)O(1/k2)），所以一阶方法还有可以提高的空间。

Nesterov方法

Nesterov (1983)认为，一阶方法会出现oscillation，也就是如下图所示的，下降路径存在很严重的“震荡”，他觉得在当次下降中，比如在x2x_2x2到x3x_3x3的下降过程中，下降方向不应该只受到x2x_2x2处梯度或者次梯度的影响，而是应该综合x2x_2x2处的梯度、次梯度以及x1x_1x1到x2x_2x2下降的“惯性”，所以Nesterov引入了一种Momentum用来修正下降方向。

基于Nesterov的思想，递推式可以写成
xk+1=xk−α∇f(xk)⏟GradientDescent+β(xk−xk−1)⏟Momentumx_{k+1}=\underbrace{x_k-\alpha\nabla f(x_k)}_{Gradient\ Descent}+\underbrace{\beta(x_k-x_{k-1})}_{Momentum}xk+1=Gradient Descentxk−α∇f(xk)+Momentumβ(xk−xk−1)

这个式子只是简单描述Nesterov的思想，具体算法实现有很多种形式，比较常用的一种是：

假设∇f\nabla f∇f是LLL-Lipschitz函数，取α=1/L\alpha=1/Lα=1/L, t1=1t_1=1t1=1，按下面的递推式迭代
tk+1=1+1+4tk22βk+1=tk−1tk+1pk+1=xk+βk+1(xk−xk−1)xk+1=pk+1−α∇f(pk+1)t_{k+1}=\frac{1+\sqrt{1+4t_k^2}}{2} \\ \beta_{k+1}=\frac{t_k-1}{t_{k+1}} \\ p_{k+1}=x_k+\beta_{k+1}(x_k-x_{k-1}) \\ x_{k+1}=p_{k+1}-\alpha\nabla f(p_{k+1})tk+1=21+1+4tk2βk+1=tk+1tk−1pk+1=xk+βk+1(xk−xk−1)xk+1=pk+1−α∇f(pk+1)

Accelerate Proximal Gradient (APG)

第二讲介绍了Proximal Gradient Descent，在分解F=f+gF=f+gF=f+g以及∇f\nabla f∇f是LLL-Lipschitz函数的条件下，Proximal Gradient Descent可以用两个公式表示出来：
wk=xk−∇f(xk)Lxk+1=proxg/L(wk)w_k=x_k-\frac{\nabla f(x_k)}{L} \\ x_{k+1}=prox_{g/L}(w_k)wk=xk−L∇f(xk)xk+1=proxg/L(wk)

其中proxgprox_gproxg是凸函数ggg的proximal operator为
proxg(w)=arg min⁡xg(x)+12∥x−w∥22prox_g(w)=\argmin_x \ g(x)+\frac{1}{2}\left\| x-w\right\|^2_2proxg(w)=xargmin g(x)+21∥x−w∥22

既然Proximal Gradient Descent的第一步也是一阶方法，我们也就可以用Nesterov方法做加速：

t1=1,tk+1=1+1+4tk22βk+1=tk−1tk+1pk+1=xk+βk+1(xk−xk−1)wk+1=pk+1−1L∇f(pk+1)xk+1=proxg/L(wk)t_1= 1,t_{k+1}=\frac{1+\sqrt{1+4t_k^2}}{2} \\ \beta_{k+1}=\frac{t_k-1}{t_{k+1}} \\ p_{k+1}=x_k+\beta_{k+1}(x_k-x_{k-1}) \\ w_{k+1}=p_{k+1}-\frac{1}{L}\nabla f(p_{k+1}) \\ x_{k+1}=prox_{g/L}(w_k)t1=1,tk+1=21+1+4tk2βk+1=tk+1tk−1pk+1=xk+βk+1(xk−xk−1)wk+1=pk+1−L1∇f(pk+1)xk+1=proxg/L(wk)

这个算法被称为Accelerate Proximal Gradient (APG)，这里附一个用APG做Basis Pursuit问题
min⁡x∥x∥1s.t.y=Ax\min_x \ \left\| x\right\|_1 \\ s.t.\ \ y=Axxmin ∥x∥1s.t. y=Ax

的R代码示例（请读者自行调试）

soft.threshold <- function(x,tau){len0 = length(x)y = rep(0,len0)for(index0 in 1:len0){if(x[index0]>=tau){y[index0] = x[index0]-tau}if(x[index0]<=-tau){y[index0] = x[index0]+tau}}return(y)
}APG <- function(y,A,x0,lambda,mu,maxiter = 100,eps = 1e-3){x1 = x0t1 = 1L = abs(max(eigen(t(A)%*%A)$values))numiter = 0diff.x = 1while(numiter<maxiter & diff.x>=eps){t = 0.5*(1+sqrt(1+4*t1^2))beta = (t1 - 1)/tp = x1 + beta*(x1 - x0)w = p - (t(A)%*%lambda + mu*t(A)%*%(A%*%p-y))/Lx = soft.threshold(w,1/L)diff.x = sqrt(sum((x-x1)^2))numiter = numiter + 1x0 = x1x1 = x}return(x1)
}