UA MATH567 高维统计专题3 含L1-norm的凸优化2 Proximal Gradient Descent

Proximal Gradient Descent的公式推导
Proximal Operator
- Indicator function
- LASSO
- Nuclear Norm

对于平滑的凸函数，写一个梯度下降的求最值的算法非常简单（上一讲介绍梯度下降但还在施工中）；但非平滑的函数求不出梯度，这时可以用Proximal Gradient Descent代替梯度下降求它的最值。

Proximal Gradient Descent的公式推导

假设F(x)F(x)F(x)是一个非平滑的凸函数，它可以被分解为一个平滑函数f(x)f(x)f(x)和一个形式相对比较简单的函数g(x)g(x)g(x)，
F(x)=f(x)+g(x)F(x)=f(x)+g(x)F(x)=f(x)+g(x)

假设∇f(x)\nabla f(x)∇f(x)是L-Lipschitz函数，则
F(x)≤F^(x,xk)=f(xk)+⟨∇f(xk),x−xk⟩+L2∥x−xk∥22⏟upperboundoff(x)+g(x)F(x) \le \hat F(x,x_k) \\=\underbrace{f(x_k)+\langle \nabla f(x_k),x-x_k \rangle +\frac{L}{2} \left\| x-x_k \right\|^2_2}_{upper\ bound\ of\ f(x)}+g(x)F(x)≤F^(x,xk)=upper bound of f(x)f(xk)+⟨∇f(xk),x−xk⟩+2L∥x−xk∥22+g(x)

Proximal gradient descent的思路用F^(x,xk)\hat F(x,x_k)F^(x,xk)的最小值点作为xk+1x_{k+1}xk+1，
xk+1=arg min⁡xF^(x,xk)=arg min⁡x⟨∇f(xk),x−xk⟩+L2∥x−xk∥22+g(x)=arg min⁡xL2∥x−xk+∇f(xk)L∥22+g(x)=arg min⁡xL2∥x−wk∥22+g(x)x_{k+1}=\argmin_{x} \hat F(x,x_k) \\ =\argmin_{x} \langle \nabla f(x_k),x-x_k \rangle +\frac{L}{2} \left\| x-x_k \right\|^2_2+g(x) \\ =\argmin_{x}\frac{L}{2} \left\| x-x_k+\frac{\nabla f(x_k)}{L} \right\|^2_2+g(x) \\ = \argmin_{x}\frac{L}{2} \left\| x-w_k\right\|^2_2+g(x)xk+1=xargminF^(x,xk)=xargmin⟨∇f(xk),x−xk⟩+2L∥x−xk∥22+g(x)=xargmin2L∥∥∥∥x−xk+L∇f(xk)∥∥∥∥22+g(x)=xargmin2L∥x−wk∥22+g(x)

其中
wk=xk−∇f(xk)Lw_k=x_k-\frac{\nabla f(x_k)}{L}wk=xk−L∇f(xk)

如果g(x)g(x)g(x)可以忽略，那上式的解就是xk+1=wkx_{k+1}=w_kxk+1=wk，这样操作起来就很容易；但即使g(x)g(x)g(x)不能忽略，它的形式也比较简单，所以这个最优化应该也不难解。

定义 Proximal Operator：称凸函数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用两个公式表示出来：
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)

Proximal Operator

现在我们有了Proximal Gradient Descent的公式，假设关于目标函数F(x)F(x)F(x)的分解F(x)=f(x)+g(x)F(x)=f(x)+g(x)F(x)=f(x)+g(x)已经找出来了，那么要写出最小化F(x)F(x)F(x)的算法，我们只需要确定ggg的Proximal Operator就可以了，下面介绍几个常见的例子：

Indicator function

g(x)=ID(x)={0,x∈D+∞,x∉Dg(x)=I_D(x)= \begin{cases} 0, x \in D \\ +\infty , x \notin D\end{cases}g(x)=ID(x)={0,x∈D+∞,x∈/D

其中DDD是凸闭集；
proxg(w)=arg min⁡xg(x)+12∥x−w∥22=arg min⁡xID(x)+12∥x−w∥22=arg min⁡x∈D∥x−w∥22=PD(w)prox_g(w)=\argmin_x \ g(x)+\frac{1}{2}\left\| x-w\right\|^2_2 \\ = \argmin_x \ I_D(x)+\frac{1}{2}\left\| x-w\right\|^2_2 \\ = \argmin_{x \in D}\left\| x-w\right\|^2_2 = P_D(w)proxg(w)=xargmin g(x)+21∥x−w∥22=xargmin ID(x)+21∥x−w∥22=x∈Dargmin∥x−w∥22=PD(w)

其中PD(x)P_D(x)PD(x)是www在DDD中的投影。

LASSO

g(x)=λ∥x∥1,x∈Rdg(x)=\lambda \left\|x \right\|_1,x \in \mathbb{R}^dg(x)=λ∥x∥1,x∈Rd

proxg(w)=arg min⁡xg(x)+12∥x−w∥22=arg min⁡xλ∥x∥1+12∥x−w∥22=arg min⁡x∑i=1dλ∣xj∣+12(xj−wj)2=arg min⁡x1,⋯,xdλ∣xj∣+12(xj−wj)2=(sign(wj)max⁡{∣wj∣−λ,0})d×1=(soft(wj,λ))d×1prox_g(w)=\argmin_x \ g(x)+\frac{1}{2}\left\| x-w\right\|^2_2 \\ = \argmin_x \lambda \left\|x \right\|_1+\frac{1}{2}\left\| x-w\right\|^2_2 \\ = \argmin_x \sum_{i=1}^d \lambda |x_j|+\frac{1}{2}(x_j-w_j)^2 \\ = \argmin_{x_1,\cdots,x_d} \lambda |x_j|+\frac{1}{2}(x_j-w_j)^2 \\ = (sign(w_j)\max\{|w_j|-\lambda,0\})_{d \times 1} = (soft(w_j,\lambda))_{d \times 1}proxg(w)=xargmin g(x)+21∥x−w∥22=xargminλ∥x∥1+21∥x−w∥22=xargmini=1∑dλ∣xj∣+21(xj−wj)2=x1,⋯,xdargminλ∣xj∣+21(xj−wj)2=(sign(wj)max{∣wj∣−λ,0})d×1=(soft(wj,λ))d×1

也就是entry-wise soft-threshold of wjw_jwj，用proximal gradient descent可以写LASSO的算法（但实际上肯定是R package效率更高，比如glmnet之类的包），下面是用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)
}PGLasso <- function(y,X,lambda,beta0,L,maxiter,eps){## L >= max eigen value if X'X## choose an appropriate value for L as input numiter = 0dis = 1while((numiter < maxiter) & (dis >= eps)){w = beta0 - (1/L)*t(X)%*%(X%*%beta0 - y)beta = soft.threshold(w,lambda/L)dis = norm(beta - beta0,type = "2")beta0 = betanumiter = numiter + 1}return(beta0)
}

这里简单介绍一下PGLasso的思路：
min⁡xF(β)=12∥y−Xβ∥22⏟f(β)+λ∥β∥1⏟g(β)\min_x F(\beta)= \underbrace{\frac{1}{2}\left\| y -X\beta \right\|_2^2}_{f(\beta)}+\underbrace{\lambda \left\| \beta\right\|_1}_{g(\beta)}xminF(β)=f(β)21∥y−Xβ∥22+g(β)λ∥β∥1

因此∇f(β)\nabla f(\beta)∇f(β)的Lipschitz范数为λ1(XTX)\lambda_1(X^TX)λ1(XTX)，我们可以把它作为LLL用在算法中，也可以取一个比λ1(XTX)\lambda_1(X^TX)λ1(XTX)稍微大一点点的常数作为LLL（但收敛会稍微变慢一点），代入Proximal gradient descent的公式：
wk=βk−1L∇f(βk)=βk−XTL(Xβk−y)βk+1=proxg/L(wk)=soft(wk,λ/L)w_k = \beta_k-\frac{1}{L}\nabla f(\beta_k)=\beta_k-\frac{X^T}{L}(X\beta_k-y) \\ \beta_{k+1} = prox_{g/L}(w_k)=soft(w_k,\lambda/L)wk=βk−L1∇f(βk)=βk−LXT(Xβk−y)βk+1=proxg/L(wk)=soft(wk,λ/L)

这就是PGLasso这个函数的思路。

用模拟数据试一下这个算法，输出solution path。

DGP <- function(n,d,s,sigma,sigma_epsilon){noise <- mvrnorm(1,rep(0,n),sigma_epsilon*diag(1,n,n))beta <- rbinom(d,1,s/d)*runif(d,-1,1)U <- sigma*diag(1,d,d)X = mvrnorm(n,rep(0,d),U)y = X%*%beta+noiseSimulation.data <- list(y = y,X = X,beta = beta)return(Simulation.data)
}library(MASS)
set.seed(2021)
Data <- DGP(100,10,3,1,0.2)
lambda = seq(1,50,1)
beta.path = matrix(rep(0,100*50),100,50)
L = max(eigen(t(Data$X)%*%(Data$X))$values)
start1 <- Sys.time()
for(i in 1:50){beta.path[,i] = PGLasso(Data$y,Data$X,lambda[i],rep(1,10),L,100,1e-6)
}
end1 <- Sys.time()
print(start1-end1)
colors = c("pink","green","yellow","orange","purple","red","blue","wheat","brown","grey")
for (kk in 1:10) {plot(beta.path[kk,]~lambda, type = "l", col = colors[kk],xlab = "lambda", ylab = "coordinates of beta",xlim = c(-1,51),ylim = c(-0.6,0.15))par(new = T)
}

Nuclear Norm

g(x)=λ∥x∥∗,x∈Rd1×d2g(x)=\lambda \left\|x\right\|_*,x \in \mathbb{R}^{d_1 \times d_2}g(x)=λ∥x∥∗,x∈Rd1×d2

则
proxg(w)=Usoft(Σ,λ)VTprox_g(w)=Usoft(\Sigma,\lambda)V^Tproxg(w)=Usoft(Σ,λ)VT

其中U,Σ,VU,\Sigma,VU,Σ,V是w∈Rd1×d2w\in \mathbb{R}^{d_1 \times d_2}w∈Rd1×d2的SVD，w=UΣVTw=U\Sigma V^Tw=UΣVT；简单推导如下：
proxg(w)=arg min⁡xg(x)+12∥x−w∥F2=arg min⁡xλ∥x∥∗+12∥x−w∥F2=arg min⁡xλ∥UTxV∥∗+12∥UTxV−Σ∥F2prox_g(w)=\argmin_x \ g(x)+\frac{1}{2}\left\| x-w\right\|^2_F \\ = \argmin_x \ \lambda \left\|x\right\|_*+\frac{1}{2}\left\| x-w\right\|^2_F \\ = \argmin_x \ \lambda \left\|U^TxV\right\|_*+\frac{1}{2}\left\| U^TxV-\Sigma\right\|^2_F proxg(w)=xargmin g(x)+21∥x−w∥F2=xargmin λ∥x∥∗+21∥x−w∥F2=xargmin λ∥∥UTxV∥∥∗+21∥∥UTxV−Σ∥∥F2

其中
12∥UTxV−Σ∥F2=12(∥UTxV∥F2+∥Σ∥F2−2⟨UTxV,Σ⟩)\frac{1}{2}\left\| U^TxV-\Sigma\right\|^2_F = \frac{1}{2}(\left\| U^TxV\right\|^2_F+\left\| \Sigma\right\|^2_F-2\langle U^TxV,\Sigma \rangle)21∥∥UTxV−Σ∥∥F2=21(∥∥UTxV∥∥F2+∥Σ∥F2−2⟨UTxV,Σ⟩)

根据von Neumann(1927)不等式，
⟨UTxV,Σ⟩≤⟨σ(x),σ(Σ)⟩\langle U^TxV,\Sigma \rangle \le \langle \sigma(x),\sigma(\Sigma) \rangle⟨UTxV,Σ⟩≤⟨σ(x),σ(Σ)⟩

当且仅当UTxV=diag(σ(X))U^TxV=diag(\sigma(X))UTxV=diag(σ(X))时等式成立，则
12∥UTxV−Σ∥F2≥12∥UTxV∥F2+12∥Σ∥F2−⟨σ(x),σ(Σ)⟩=12∥x∥F2+12∥Σ∥F2−⟨σ(x),σ(Σ)⟩=12∥σ(x)∥22+12∥σ(Σ)∥22−⟨σ(x),σ(Σ)⟩=12∥σ(x)−σ(Σ)∥22\frac{1}{2}\left\| U^TxV-\Sigma\right\|^2_F \ge \frac{1}{2}\left\| U^TxV\right\|^2_F+\frac{1}{2}\left\| \Sigma\right\|^2_F-\langle \sigma(x),\sigma(\Sigma) \rangle \\ =\frac{1}{2}\left\| x\right\|^2_F+\frac{1}{2}\left\| \Sigma\right\|^2_F-\langle \sigma(x),\sigma(\Sigma) \rangle \\ = \frac{1}{2} \left\| \sigma(x) \right\|_2^2+\frac{1}{2} \left\| \sigma(\Sigma) \right\|_2^2-\langle \sigma(x),\sigma(\Sigma) \rangle = \frac{1}{2}\left\| \sigma(x) - \sigma(\Sigma) \right\|_2^221∥∥UTxV−Σ∥∥F2≥21∥∥UTxV∥∥F2+21∥Σ∥F2−⟨σ(x),σ(Σ)⟩=21∥x∥F2+21∥Σ∥F2−⟨σ(x),σ(Σ)⟩=21∥σ(x)∥22+21∥σ(Σ)∥22−⟨σ(x),σ(Σ)⟩=21∥σ(x)−σ(Σ)∥22

所以
arg min⁡xλ∥UTxV∥∗+12∥UTxV−Σ∥F2=arg min⁡xλ∥σ(x)∥1+12∥σ(x)−σ(Σ)∥22=soft(σ(Σ),λ)\argmin_x \ \lambda \left\|U^TxV\right\|_*+\frac{1}{2}\left\| U^TxV-\Sigma\right\|^2_F \\ = \argmin_x \lambda \left\| \sigma(x) \right\|_1 + \frac{1}{2}\left\| \sigma(x) - \sigma(\Sigma) \right\|_2^2 \\ = soft(\sigma(\Sigma),\lambda)xargmin λ∥∥UTxV∥∥∗+21∥∥UTxV−Σ∥∥F2=xargminλ∥σ(x)∥1+21∥σ(x)−σ(Σ)∥22=soft(σ(Σ),λ)