Group Lasso

In this article, I compare proximal gradient descent method with accelerated proximal gradient descent method and accelerated gradient descent method smoothed problem. The code is show in part (c) and some experiments are conducted in part (d)

(a) Smoothing

We need to Design a \(1/\mu\)-smooth approximation of \(h(\beta)\)

\[h(\beta)=\sum_{j=1}^J\sqrt{d_j}||\beta_{(j)}||_2\]

We know that the subgradient of \(h(\beta)\) is

\[\partial h_j(\beta)=\sqrt{d_j}\frac{\beta_{(j)}}{||\beta_{(j)}||_2}\leq\sqrt{d_j }I_j\] The subgradient is bounded, so we know that \(h(\beta)\) is L-Lipschitz

We then use Moreau-Yosida regularization to get the \(1/\mu\)-smooth approximation of \(h(\beta)\)

\[h_{\mu}(\beta)=min_y\big(h(y)+\frac{1}{2\mu}||\beta-y||^2\big)=min_y\big(\sum_{j=1}^J\sqrt{d_j}||y_{(j)}||_2+\frac{1}{2\mu}||\beta-y||^2\big)\]

We use proximal operator to solve this, as we have proved it in the homework4, the proximal operator for \(h_{\mu}(\beta)\) is

\[y_{(j)}=Prox_{\mu h_j}(\beta_{(j)})=\big(1-\frac{\mu\sqrt{d_j}}{||\beta_{(j)}||_2}\big)_+\beta_{(j)}\] Where \((x)_+=max(x,0)\).

So, the Moreau envelop is

\[h_{\mu}(\beta)=\sum_{j=1}^J\sqrt{d_j}||\big(1-\ frac{\mu\sqrt{d_j}}{||\beta_{(j)}||_2}\big)_+\beta_{(j)}||_2+\frac{1}{2\mu}||\beta-y||^2\]

Where \(y=[y_{(1)},....y_{(j)}]\)

Because we know that \(h(\beta)\) is L-Lipschitz, \(h_{\mu}(\beta)\) is the \(1/\mu\)-smooth approximation of \(h(\beta)\)

About the choice of \(\mu\), we can fix \(d\), \(n\), \(d\) and \(\lambda\), then, search the best \(\mu\) from 0 to K, step is \(\rho\), after we test the model and get the MSE of the test data set, we can obtain the best \(\mu\) by choose the \(\mu\) of the lowest MSE.

(b) Algorithm steps

(i) FISTA

As in previous work, by defining \(\widetilde{Y_i}=\frac{Y_i+1}{2}\) and \(h(X_i^T\beta)=\frac{1}{1+exp(-X_i^T\beta)}\), the likelihood function is

\[f(\beta)=\sum_{i=1}^n\widetilde{Y_i}X_i^T\beta+log(h(-X_i^T\beta))\]

And the gradient is

\[\triangledown f(\beta)=\sum_{i=1}^n[\widetilde{Y_i}-h(X_i^T\beta)]X_i\]

For the accelerated proximal gradient descent method, I use FISTA method, below is steps. For the objective function

\[min_{\beta} \space f(\beta)+\lambda h(\beta)\]

  1. Choose initial point \(\beta^{(0)}\) and \(s^{(0)}=1\), let \(k=0\)

  2. Compute \(t_k\) and \(\triangledown(f(\beta)\)

  3. Then update \(x\) as \(x^{(k+1)}=Prox_{t_k\lambda h}(\beta^{(k)}-t_k\triangledown f(\beta^{(k)}))\)

  4. \(s^{(k+1)}=\frac{1+\sqrt{1+4s^{(k)^2}}}{2}\)

  5. Update \(\beta^{(k+1)}=x^{(k+1)}+\frac{s^{(k)}-1}{s^{(k+1)}}(x^{(k+1)}-x^{(k)})\)

  6. \(k=k+1\), return to 2

  7. Output: \(\beta^{(k+1)}\)

(ii) Nesterov’s accelerated gradient descent method

The objective function is

\[F(\beta)=f(\beta)+\lambda h_{\mu}(\beta)=\sum_{i=1}^n\widetilde{Y_i}X_i^T\beta+log(h(-X_i^T\beta))+\lambda(\sum_{j=1}^J\sqrt{d_j}||\big(1-\frac{\mu\sqrt{d_j}}{||\beta_{(j)}||_2}\big)_+\beta_{(j)}||_2+\frac{1}{2\mu}||\beta-y||^2)\]

\[\triangledown F(\beta)=\sum_{i=1}^n[\widetilde{Y_i}-h(X_i^T\beta)]X_i+\lambda\triangledown h_{\mu}(\beta)\]

If \(1-\frac{\mu\sqrt{d_j}}{||\beta_{(j)}||_2}>0\)

We know that \(||y_{(j)}||=||\beta_{(j)}||-\mu \sqrt{d_j}\), and

\[\triangledown_{\beta_{(j)}}||y_{(j)}||=\frac{\beta_{(j)}}{||\beta_{(j)}||}\]

\[\triangledown_{\beta_{(j)}} h_{\mu}(\beta)=\sqrt{d_j}\triangledown_{\beta_{(j)}}||y_{(j)}||+\triangledown_{\beta_{(j)}}\frac{1}{2\mu}||\frac{\mu\sqrt{d_j}}{||\beta_{(j)}||_2}\beta_{(j)}||^2\\=\sqrt{d_j}\frac{\beta_{(j)}}{||\beta_{(j)}||}+0\]

If \(1-\frac{\mu\sqrt{d_j}}{||\beta_{(j)}||_2}<0\)

\[\triangledown_{\beta_{(j)}} h_{\mu}(\beta)=\frac{\beta_{(j)}}{\mu}\]

For Nesterov’s accelerated gradient descent method smoothed problem, we get the algorithm below.

  1. Choose initial point \(\beta^{(0)}\) and \(s^{(0)}=1\), let \(k=0\)

  2. Compute \(t_k\) and \(\triangledown(F(\beta))\)

  3. Then update \(x\) as \(x^{(k+1)}=\beta^{(k)}-t_k\triangledown F(\beta^{(k)})\)

  4. Update \(\beta^{(k+1)}=x^{(k+1)}+\frac{k}{k+3}(x^{(k+1)}-x^{(k)})\)

  5. \(k=k+1\), return to 2

  6. Output: \(\beta^{(k+1)}\)

(c) Code

Accelerated proximal gradient descent method

Group_lasso=function(beta,sstep,lam,group_d){
  temp=beta
  
  for (i in c(1:(length(group_d)-1))){
    if (sqrt(sum(beta[(group_d[i]+1):sum(group_d[1:(i+1)])]^2))>lam*sqrt(group_d[i+1]))
      {
    temp[(group_d[i]+1):sum(group_d[1:(i+1)])]=beta[(group_d[i]+1):sum(group_d[1:(i+1)])]-sstep*lam*sqrt(group_d[i])*(beta[(group_d[i]+1):sum(group_d[1:(i+1)])]/sqrt(sum(beta[(group_d[i]+1):sum(group_d[1:(i+1)])]^2)))
    }
    else{
    temp[(group_d[i]+1):sum(group_d[1:(i+1)])]=0
    }

  }
  return(temp)
}


GradientDescent_Accelerated<-function(x,y,maxiter,alpha,beta,theta_s,lam)
{
  sigmoid <- function(z) { 1 / (1 + exp(-z))}
  n<-nrow(x)
  d<-ncol(x)
  theta=matrix(rep(0,d),d,1) 
  iter<-1
  error_sum=rep(1,maxiter)
  s=rep(0,maxiter)
  group_d=matrix(c(0,5,0,4,8,0,3,2,5,0,3))
   
while(iter<maxiter){
    iter<-iter+1
    s[iter]=(1+sqrt(1+4*s[iter-1]^2))/2
    h<-sigmoid(x%*%theta)  
    grad=(1/n)*t(x)%*%(h-y) 
    des= grad    
    sstep=1
    new_theta<-theta-sstep*des #5*1
    new_theta=Group_lasso(new_theta,sstep,lam,group_d)
    new_theta_temp=new_theta+(s[iter-1]-1)/s[iter]*(new_theta-theta)
    new_theta=new_theta_temp
    new_h<-sigmoid(x%*%new_theta) # 2000*1
    costfunction<-(1/n)*(-t(y)%*%log(h))-(t(1-y)%*%log(1-h))
    new_costfunction<-(1/n)*(-t(y)%*%log(new_h))-(t(1-y)%*%log(1-new_h))
    while(new_costfunction>costfunction-alpha*sstep*t(grad)%*%grad){  # backtrack line search
          sstep<-sstep*beta
          new_theta<-theta-sstep*des
          new_h<-sigmoid(x%*%new_theta)
          new_costfunction<-(1/n)*(-t(y)%*%log(new_h))-(t(1-y)%*%log(1-new_h))
          }
    theta<-new_theta 
    error_sum[iter-1]=log(sqrt(sum((theta-theta_s)^2)))
   
}
      out=list(theta_output=theta,error_output=error_sum)
      return(out)
}

Nesterov’s accelerated gradient descent method smoothed problem

Group_lasso_grad=function(beta,sstep,mu,lam,group_d){
  temp=beta
  
  for (i in c(1:(length(group_d)-1))){
    if (sqrt(sum(beta[(group_d[i]+1):sum(group_d[1:(i+1)])]^2))>lam*sqrt(group_d[i+1]))
      {
    temp[(group_d[i]+1):sum(group_d[1:(i+1)])]=sqrt(group_d[i])*(beta[(group_d[i]+1):sum(group_d[1:(i+1)])]/sqrt(sum(beta[(group_d[i]+1):sum(group_d[1:(i+1)])]^2)))
    }
    else{
    temp[(group_d[i]+1):sum(group_d[1:(i+1)])]=beta[(group_d[i]+1):sum(group_d[1:(i+1)])]/mu
    }

  }
  return(temp)
}


Nesterov<-function(x,y,maxiter,alpha,beta,mu,lam,theta_s)
{
  sigmoid <- function(z) { 1 / (1 + exp(-z))}
  n<-nrow(x)
  d<-ncol(x)
  theta=matrix(rep(0,d),d,1) 
  iter<-1
  k=0
  error_sum=rep(1,maxiter)
  s=rep(0,maxiter)
  group_d=matrix(c(0,5,0,4,8,0,3,2,5,0,3))
   
while(iter<maxiter){
    iter<-iter+1
    s[iter]=(1+sqrt(1+4*s[iter-1]^2))/2
    h<-sigmoid(x%*%theta)  
    grad=(1/n)*t(x)%*%(h-y)+lam*Group_lasso_grad(theta,sstep,mu,lam,group_d)
    des= grad    
    sstep=1
    new_theta<-theta-sstep*des #5*1
    new_theta_temp=new_theta+(k/(k+3))*(new_theta-theta)
    new_theta=new_theta_temp
    new_h<-sigmoid(x%*%new_theta) 
    costfunction<-(1/n)*(-t(y)%*%log(h))-(t(1-y)%*%log(1-h))
    new_costfunction<-(1/n)*(-t(y)%*%log(new_h))-(t(1-y)%*%log(1-new_h))
    while(new_costfunction>costfunction-alpha*sstep*t(grad)%*%grad){  # backtrack line search
          sstep<-sstep*beta
          new_theta<-theta-sstep*des
          new_h<-sigmoid(x%*%new_theta)
          new_costfunction<-(1/n)*(-t(y)%*%log(new_h))-(t(1-y)%*%log(1-new_h))
          }
    theta<-new_theta
    k=k+1
    error_sum[iter-1]=log(sqrt(sum((theta-theta_s)^2)))
   
}
      out=list(theta_output=theta,error_output=error_sum)
      return(out)
}

(d) Results

The group is (5,0,4,8,0,3,2,5,0,3)

The experiments for Accelerated proximal gradient descent method

For 2000 observations and each observation has 30 features, we set lambda to (0.1, 0.3, 0.5, 0.7), we have the following results.

n1=2000
d1=30
theta_s<-matrix(rep(1/sqrt(d1),d1),d1,1)
mean<-rep(0,d1)
sigma<-diag(d1)
x = mvrnorm(n1, mean, sigma)
y=1/(1+exp(-x%*%theta_s))
judge_5=function(y){
  if(y>0.5){
    x=1
  }
  else{
    x=0
  }
  return(x)
}

y_f=lapply(y,judge_5 )
y_f=as.vector(unlist(y_f))
Y=matrix(0,n1)
Z=matrix(0,n1)
for (i in 1:n1){
  Z[i]=rbern(1,prob=y[i])
  Y[i]=2*Z[i]-1
}
out1_a=GradientDescent_Accelerated(x,Y,30,alpha=0.25,beta=0.9,lam=0.1,theta_s)
out2_a=GradientDescent_Accelerated(x,Y,30,alpha=0.25,beta=0.9,lam=0.3,theta_s)
out3_a=GradientDescent_Accelerated(x,Y,30,alpha=0.25,beta=0.9,lam=0.5,theta_s)
out4_a=GradientDescent_Accelerated(x,Y,30,alpha=0.25,beta=0.9,lam=0.7,theta_s)
plot(out1_a$error_output[-30],col=1,lty=1,type="l")
lines(out2_a$error_output[-30],col=2,lty=2,type="l")
lines(out3_a$error_output[-30],col=3,lty=3,type="l")
lines(out4_a$error_output[-30],col=4,lty=4,type="l")
legend("topright", legend=c("lambda=0.1", "lambda=0.3","lambda=0.5","lambda=0.7"),col=c(1,2,3,4),lty=c(1,2,3,4), cex=0.8)

For n=2000, 3000, 4000, and we set lambda=0.3 and d=30

n1=3000
d1=30
theta_s<-matrix(rep(1/sqrt(d1),d1),d1,1)
mean<-rep(0,d1)
sigma<-diag(d1)
x = mvrnorm(n1, mean, sigma)
y=1/(1+exp(-x%*%theta_s))
judge_5=function(y){
  if(y>0.5){
    x=1
  }
  else{
    x=0
  }
  return(x)
}

y_f=lapply(y,judge_5 )
y_f=as.vector(unlist(y_f))
Y=matrix(0,n1)
Z=matrix(0,n1)
for (i in 1:n1){
  Z[i]=rbern(1,prob=y[i])
  Y[i]=2*Z[i]-1
}
out5_a=GradientDescent_Accelerated(x,Y,30,alpha=0.25,beta=0.9,lam=0.3,theta_s)


n1=4000
d1=30
theta_s<-matrix(rep(1/sqrt(d1),d1),d1,1)
mean<-rep(0,d1)
sigma<-diag(d1)
x = mvrnorm(n1, mean, sigma)
y=1/(1+exp(-x%*%theta_s))
judge_5=function(y){
  if(y>0.5){
    x=1
  }
  else{
    x=0
  }
  return(x)
}

y_f=lapply(y,judge_5 )
y_f=as.vector(unlist(y_f))
Y=matrix(0,n1)
Z=matrix(0,n1)
for (i in 1:n1){
  Z[i]=rbern(1,prob=y[i])
  Y[i]=2*Z[i]-1
}
out6_a=GradientDescent_Accelerated(x,Y,30,alpha=0.25,beta=0.9,lam=0.3,theta_s)

plot(out2_a$error_output[-30],col=1,lty=1,type="l",ylim=c(-1.6,-0.5))
lines(out5_a$error_output[-30],col=2,lty=2,type="l")
lines(out6_a$error_output[-30],col=3,lty=3,type="l")

legend("topright", legend=c("n=2000", "n=3000","n=4000"),col=c(1,2,3),lty=c(1,2,3), cex=0.8)

The experiments for Nesterov’s accelerated gradient descent method smoothed problem

For 2000 observations and each observation has 30 features, we set lambda to (0.1, 0.3, 0.5, 0.7), we have the following results.

out1_n=Nesterov(x,Y,30,alpha=0.25,beta=0.9,mu=1,lam=0.1,theta_s)
out2_n=Nesterov(x,Y,30,alpha=0.25,beta=0.9,mu=1,lam=0.3,theta_s)
out3_n=Nesterov(x,Y,30,alpha=0.25,beta=0.9,mu=1,lam=0.5,theta_s)
out4_n=Nesterov(x,Y,30,alpha=0.25,beta=0.9,mu=1,lam=0.7,theta_s)
plot(out1_n$error_output[-30],col=1,lty=1,type="l")
lines(out2_n$error_output[-30],col=2,lty=2,type="l")
lines(out3_n$error_output[-30],col=3,lty=3,type="l")
lines(out4_n$error_output[-30],col=4,lty=4,type="l")
legend("topright", legend=c("lambda=0.1", "lambda=0.3","lambda=0.5","lambda=0.7"),col=c(1,2,3,4),lty=c(1,2,3,4), cex=0.8)

Then, we check \(\mu\), we set \(\mu=1,2,3,4\), and fix \(n\), \(d\), \(\lambda=0.3\)

out5_n=Nesterov(x,Y,30,alpha=0.25,beta=0.9,mu=0.5,lam=0.3,theta_s)
out6_n=Nesterov(x,Y,30,alpha=0.25,beta=0.9,mu=1,lam=0.3,theta_s)
out7_n=Nesterov(x,Y,30,alpha=0.25,beta=0.9,mu=2,lam=0.3,theta_s)
out8_n=Nesterov(x,Y,30,alpha=0.25,beta=0.9,mu=3,lam=0.3,theta_s)
out9_n=Nesterov(x,Y,30,alpha=0.25,beta=0.9,mu=4,lam=0.3,theta_s)
plot(out5_n$error_output[-30],col=1,lty=1,type="l")
lines(out6_n$error_output[-30],col=2,lty=2,type="l")
lines(out7_n$error_output[-30],col=3,lty=3,type="l")
lines(out8_n$error_output[-30],col=4,lty=4,type="l")
lines(out9_n$error_output[-30],col=5,lty=5,type="l")
legend("topright", legend=c("mu=0.5", "mu=1","mu=2","mu=3","mu=4"),col=c(1,2,3,4,5),lty=c(1,2,3,4,5), cex=0.8)

Proximal Gradient method

We use our results in previous work

Group_lasso=function(beta,sstep,lam,group_d){
  temp=beta
  
  for (i in c(1:(length(group_d)-1))){
    if (sqrt(sum(beta[(group_d[i]+1):sum(group_d[1:(i+1)])]^2))>lam*sqrt(group_d[i+1]))
      {
    temp[(group_d[i]+1):sum(group_d[1:(i+1)])]=beta[(group_d[i]+1):sum(group_d[1:(i+1)])]-sstep*lam*sqrt(group_d[i])*(beta[(group_d[i]+1):sum(group_d[1:(i+1)])]/sqrt(sum(beta[(group_d[i]+1):sum(group_d[1:(i+1)])]^2)))
    }
    else{
    temp[(group_d[i]+1):sum(group_d[1:(i+1)])]=0
    }

  }
  return(temp)
}


GradientDescent_test<-function(x,y,maxiter,alpha,beta,theta_s)
{
  sigmoid <- function(z) { 1 / (1 + exp(-z))}
  n<-nrow(x)
  d<-ncol(x)
  theta=matrix(rep(0,d),d,1) 
  iter<-1
  error_sum=rep(0,maxiter)
  group_d=matrix(c(0,5,0,4,8,0,3,2,5,0,3))
   
while(iter<maxiter){
    iter<-iter+1
    h<-sigmoid(x%*%theta)  #2000*1
    grad=(1/n)*t(x)%*%(h-y) #5*1
    des= grad    # 5*1
    sstep=1
    new_theta<-theta-sstep*des #5*1
    new_theta=Group_lasso(new_theta,sstep,lam=0.5,group_d)
    new_h<-sigmoid(x%*%new_theta) # 2000*1
    costfunction<-(1/n)*(-t(y)%*%log(h))-(t(1-y)%*%log(1-h))
    new_costfunction<-(1/n)*(-t(y)%*%log(new_h))-(t(1-y)%*%log(1-new_h))
    while(new_costfunction>costfunction-alpha*sstep*t(grad)%*%grad){
          sstep<-sstep*beta
          new_theta<-theta-sstep*des
          new_h<-sigmoid(x%*%new_theta)
          new_costfunction<-(1/n)*(-t(y)%*%log(new_h))-(t(1-y)%*%log(1-new_h))
          }
    theta<-new_theta 
    error_sum[iter-1]=log(sqrt(sum((theta-theta_s)^2)))
   
}
      out=list(theta_output=theta,error_output=error_sum)
      return(out)
}
n1=3000
d1=30
theta_s<-matrix(rep(1/sqrt(d1),d1),d1,1)
mean<-rep(0,d1)
sigma<-diag(d1)
x = mvrnorm(n1, mean, sigma)
y=1/(1+exp(-x%*%theta_s))
judge_5=function(y){
  if(y>0.5){
    x=1
  }
  else{
    x=0
  }
  return(x)
}

y_f=lapply(y,judge_5 )
y_f=as.vector(unlist(y_f))
Y=matrix(0,n1)
Z=matrix(0,n1)
for (i in 1:n1){
  Z[i]=rbern(1,prob=y[i])
  Y[i]=2*Z[i]-1
}
out=GradientDescent_test(x,Y,30,alpha=0.25,beta=0.9,theta_s)

Finally we use \(n=3000\), \(d=30\), \(group=(5,0,4,8,0,3,2,5,0,3)\) to compare these three methods

plot(out$error_output[-30],col=1,lty=1,type="l")
lines(out3_a$error_output[-30],col=2,lty=2,type="l")
lines(out6_n$error_output[-30],col=3,lty=3,type="l")

legend("topright", legend=c("Proximal gradient method", "Accelerated proximal gradient method","Nesterov's accelerated method"),col=c(1,2,3),lty=c(1,2,3), cex=0.8)

We can conclude that Accelerated proximal gradient method and Nesterov’s accelerated method converge faster than proximal gradient method, but the accuracy of proximal gradient method is better than other two methods.