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)$和一個形式相對比較簡單的函數$g(x)$，
$$F(x)=f(x)+g(x)$$

假設$?f(x)$是L-Lipschitz函數，則
$$\begin{gathered} F(x)\le \hat{F}(x,x_k) \\ =\underset{upperboundoff(x)}{f(x_k)+??f(x_k),x?x_k?+\frac{L}{2}{\Vert x?x_k\Vert }_2^2}+g(x) \end{gathered}$$

Proximal gradient descent的思路用$\hat{F}(x,x_k)$的最小值點作為$x_{k+1}$，
$$\begin{gathered} x_{k+1}=\underset{x}{\mathrm{argmin}}\hat{F}(x,x_k) \\ =\underset{x}{\mathrm{argmin}}??f(x_k),x?x_k?+\frac{L}{2}{\Vert x?x_k\Vert }_2^2+g(x) \\ =\underset{x}{\mathrm{argmin}}\frac{L}{2}{\Vert x?x_k+\frac{?f(x_k)}{L}\Vert }_2^2+g(x) \\ =\underset{x}{\mathrm{argmin}}\frac{L}{2}{\Vert x?w_k\Vert }_2^2+g(x) \end{gathered}$$

其中
$$w_k=x_k?\frac{?f(x_k)}{L}$$

如果$g(x)$可以忽略，那上式的解就是$x_{k+1}=w_k$，這樣操作起來就很容易；但即使$g(x)$不能忽略，它的形式也比較簡單，所以這個最優化應該也不難解。

定義 Proximal Operator：稱凸函數$g$的proximal operator為
$$pro{x}_g(w)=\underset{x}{\mathrm{argmin}}g(x)+\frac{1}{2}{\Vert x?w\Vert }_2^2$$

評注 在這個定義的幫助下，我們可以把Proximal Gradient Descent用兩個公式表示出來：
$$\begin{gathered} w_k=x_k?\frac{?f(x_k)}{L} \\ {x}_{k+1}=pro{x}_{g/L}(w_k) \end{gathered}$$

Proximal Operator

現在我們有了Proximal Gradient Descent的公式，假設關于目標函數$F(x)$的分解$F(x)=f(x)+g(x)$已經找出來了，那么要寫出最小化$F(x)$的算法，我們只需要確定$g$的Proximal Operator就可以了，下面介紹幾個常見的例子：

Indicator function

$$g(x)=I_D(x)=\{\begin{array}{l}0,x\in D\\ +\infty ,x\in /D\end{array}$$

其中$D$是凸閉集；
$$\begin{gathered} pro{x}_g(w)=\underset{x}{\mathrm{argmin}}g(x)+\frac{1}{2}{\Vert x?w\Vert }_2^2 \\ =\underset{x}{\mathrm{argmin}}{I}_D(x)+\frac{1}{2}{\Vert x?w\Vert }_2^2 \\ =\underset{x\in D}{\mathrm{argmin}}{\Vert x?w\Vert }_2^2=P_D(w) \end{gathered}$$

其中$P_D(x)$是$w$在$D$中的投影。

LASSO

$$g(x)=\lambda {\Vert x\Vert }_1,x\in {\mathbb{R}}^d$$

$$\begin{gathered} pro{x}_g(w)=\underset{x}{\mathrm{argmin}}g(x)+\frac{1}{2}{\Vert x?w\Vert }_2^2 \\ =\underset{x}{\mathrm{argmin}}\lambda {\Vert x\Vert }_1+\frac{1}{2}{\Vert x?w\Vert }_2^2 \\ =\underset{x}{\mathrm{argmin}}\sum _{i=1}^d\lambda |x_j|+\frac{1}{2}(x_j?w_j{)}^2 \\ =\underset{x_1,?,x_d}{\mathrm{argmin}}\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} \end{gathered}$$

也就是entry-wise soft-threshold of $w_j$，用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的思路：
$$\underset{x}{\min }F(\beta )=\underset{f(\beta )}{\frac{1}{2}{\Vert y?X\beta \Vert }_2^2}+\underset{g(\beta )}{\lambda {\Vert \beta \Vert }_1}$$

其中$?f(\beta )=X^{T}X\beta ?X^{T}y$，我們需要計算它的Lipschitz范數
$$\begin{gathered} {\Vert ?f(\beta )??f({\beta }^{\prime })\Vert }_2={\Vert X^{T}X(\beta ?{\beta }^{\prime })\Vert }_2 \\ \le {\Vert X^{T}X\Vert }_2{\Vert \beta ?{\beta }^{\prime }\Vert }_2\le {\Vert X^{T}X\Vert }_1{\Vert \beta ?{\beta }^{\prime }\Vert }_2 \\ ={\lambda }_1(X^{T}X){\Vert \beta ?{\beta }^{\prime }\Vert }_2 \end{gathered}$$

因此$?f(\beta )$的Lipschitz范數為${\lambda }_1(X^{T}X)$，我們可以把它作為$L$用在算法中，也可以取一個比${\lambda }_1(X^{T}X)$稍微大一點點的常數作為$L$（但收斂會稍微變慢一點），代入Proximal gradient descent的公式：
$$\begin{gathered} w_k={\beta }_k?\frac{1}{L}?f({\beta }_k)={\beta }_k?\frac{X^T}{L}(X{\beta }_k?y) \\ {\beta }_{k+1}=pro{x}_{g/L}(w_k)=soft(w_k,\lambda /L) \end{gathered}$$

這就是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)=\lambda {\Vert x\Vert }_{?},x\in {\mathbb{R}}^{d_1\times d_2}$$

則
$$pro{x}_g(w)=Usoft(\Sigma ,\lambda )V^T$$

其中$U,\Sigma ,V$是$w\in {\mathbb{R}}^{d_1\times d_2}$的SVD，$w=U\Sigma V^T$；簡單推導如下：
$$\begin{gathered} pro{x}_g(w)=\underset{x}{\mathrm{argmin}}g(x)+\frac{1}{2}{\Vert x?w\Vert }_F^2 \\ =\underset{x}{\mathrm{argmin}}\lambda {\Vert x\Vert }_{?}+\frac{1}{2}{\Vert x?w\Vert }_F^2 \\ =\underset{x}{\mathrm{argmin}}\lambda {\Vert U^{T}xV\Vert }_{?}+\frac{1}{2}{\Vert U^{T}xV?\Sigma \Vert }_F^2 \end{gathered}$$

其中
$$\frac{1}{2}{\Vert U^{T}xV?\Sigma \Vert }_F^2=\frac{1}{2}({\Vert U^{T}xV\Vert }_F^2+{\Vert \Sigma \Vert }_F^2?2?U^{T}xV,\Sigma ?)$$

根據von Neumann(1927)不等式，
$$?U^{T}xV,\Sigma ?\le ?\sigma (x),\sigma (\Sigma )?$$

當且僅當$U^{T}xV=diag(\sigma (X))$時等式成立，則
$$\begin{gathered} \frac{1}{2}{\Vert U^{T}xV?\Sigma \Vert }_F^2\ge \frac{1}{2}{\Vert U^{T}xV\Vert }_F^2+\frac{1}{2}{\Vert \Sigma \Vert }_F^2??\sigma (x),\sigma (\Sigma )? \\ =\frac{1}{2}{\Vert x\Vert }_F^2+\frac{1}{2}{\Vert \Sigma \Vert }_F^2??\sigma (x),\sigma (\Sigma )? \\ =\frac{1}{2}{\Vert \sigma (x)\Vert }_2^2+\frac{1}{2}{\Vert \sigma (\Sigma )\Vert }_2^2??\sigma (x),\sigma (\Sigma )?=\frac{1}{2}{\Vert \sigma (x)?\sigma (\Sigma )\Vert }_2^2 \end{gathered}$$

所以
$$\begin{gathered} \underset{x}{\mathrm{argmin}}\lambda {\Vert U^{T}xV\Vert }_{?}+\frac{1}{2}{\Vert U^{T}xV?\Sigma \Vert }_F^2 \\ =\underset{x}{\mathrm{argmin}}\lambda {\Vert \sigma (x)\Vert }_1+\frac{1}{2}{\Vert \sigma (x)?\sigma (\Sigma )\Vert }_2^2 \\ =soft(\sigma (\Sigma ),\lambda ) \end{gathered}$$

總結

以上是生活随笔為你收集整理的UA MATH567 高维统计专题3 含L1-norm的凸优化2 Proximal Gradient Descent的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。