Some Gradient Method example

Gradient Descent Using Backtracking Line Search

I use backtracking line search gradient descent method to optimize logistic regression.

We consider \((X_i,Y_i)\), and i=1,2,…,n, and \(X_i\in R^d\), and \(Y_i\in \{-1,1\}\). We assume that \(X_i\) is from standard normal distribution. We assume \(\beta^*=(1/\sqrt(d),1/\sqrt(d),...,1/\sqrt(d)\), and \(Y_i\) is from Bernoulli distribution.

\[P(Y_i=1)=1-P(Y_i=-1)=\frac{1}{1+exp(-X_i^T\beta^*)}\]

We use gradient descent and Backtracking Line Search to find the maximum likelihood estimator of this model.

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,100)
   
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_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)
}

Then, to test the function, I first construct a data set where n=2000, number of features is 5

n1=2000
d1=5
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,100,alpha=0.25,beta=0.9,theta_s)
plot(out$error_output[-100])

out$theta_output

##           [,1]
## [1,] 0.5373592
## [2,] 0.4608630
## [3,] 0.4735111
## [4,] 0.4378819
## [5,] 0.3203328

Frank Wolfe Method and Projected Gradient Descent Method

We consider \((X_i,Y_i)\), and i=1,2,…,n, and \(X_i\in R^d\), and \(Y_i \in R\). We assume that \((X_i,Y_i)\) satisfies the following linear model.

\[Y_i=X_i^T\beta+\sigma_i\] Where \(\sigma_i\) is drawn from standard normal distribution. We can use optimaization method to solve above problem.

\[min_\beta\sum_{i=1}^n(Y_i-X_i^T\beta)^2\] \[s.t. \space ||\beta||_1<=r\]

The first part is Frank wolfe method to solve optimization problem.

Frank_Wolfe<-function(x,y,maxiter,r)
{
  n<-nrow(x)
  d<-ncol(x)
  theta=matrix(rep(0,d),d,1) 
  iter<-0
  theta_trace=matrix(rep(0,d*maxiter),maxiter,d)
  
   
while(iter<maxiter){
    iter<-iter+1
    y_k=matrix(rep(0,d),d,1) 
    grad=2*t(-x)%*%(y-x%*%theta) #d*n*n*1=d*1
    #A=matrix(sign_func(grad),nrow=1,byrow=T)
    #lp.sol<-lp(direction="min",objective.in=grad,const.mat=A,const.dir="<=",const.rhs=r,compute.sens=1)
    #y_k=lp.sol$solution
    #y_k[which.min(grad)]=r
    if (abs(min(grad))>=abs(max(grad))){
      y_k[which.min(grad)]=1
    }else if(abs(min(grad))<abs(max(grad))){
      y_k[which.max(grad)]=-1
    }
    step=2/(iter+2)
    theta_new=(1-step)*theta+step*y_k
    theta_trace[iter,]=theta_new
    theta=theta_new
    
}
      return(theta_trace)
}

The second part is projected gradient descent method.

L1_ball_pro function is to compute the projection of L1 norm.

L1_ball_pro=function(h,r){
  temp1=sort(abs(h),decreasing=TRUE)
  temp2=cumsum(temp1)-r
  temp3=temp1-temp2/c(1:length(temp1))
  p=length(temp3[temp3>0])
  t=(sum(temp1[1:p])-r)/p
  h_abs=abs(h)-t
  h_abs[(h_abs)<0]=0
  theta_new=sign(h)*h_abs
  return(theta_new)
}

Projected_desc<-function(x,y,maxiter,r)
{
  n<-nrow(x)
  d<-ncol(x)
  theta=matrix(rep(0,d),d,1) 
  iter<-0
  theta_trace=matrix(rep(0,d*maxiter),maxiter,d)
  
   
while(iter<maxiter){
    iter<-iter+1
    grad=2*t(x)%*%(x%*%theta-y) #d*n*n*1=d*1
    grad=grad/n
    step=1/100000
    h=theta-step*grad
    theta_trace[iter,]=theta
    if (sum(abs(h))>r){
      theta_new=L1_ball_pro(h,r)
    }else if(sum(abs(h))<=r){
      theta_new=h
    }
    theta_trace[iter,]=theta_new
    theta=theta_new
}
    return(theta_trace)
}

By using fat data set in faraway library, I use two features as predictor

The first part is Frank Wolfe method.

y=as.matrix(fat$siri)
x=as.matrix(fat[,4:5])
t1<-Sys.time()
out_1=Frank_Wolfe(x,y,5000,1)
t2<-Sys.time()
time_1=t2-t1
tail(out_1)

##               [,1]       [,2]
## [4995,] 0.07511176 0.09004490
## [4996,] 0.07508170 0.09040903
## [4997,] 0.07505166 0.09077294
## [4998,] 0.07502164 0.09033663
## [4999,] 0.07499164 0.09070042
## [5000,] 0.07496166 0.09026431

Then, I use projected gradient descent method

y=as.matrix(fat$siri)
x=as.matrix(fat[,4:5])
t1<-Sys.time()
out_2=Projected_desc(x,y,5000,1)
t2<-Sys.time()
tail(out_2)

##               [,1]       [,2]
## [4995,] 0.07960534 0.08935974
## [4996,] 0.07960534 0.08935974
## [4997,] 0.07960534 0.08935974
## [4998,] 0.07960534 0.08935974
## [4999,] 0.07960534 0.08935974
## [5000,] 0.07960534 0.08935974

time_2=t2-t1

I use lm function output as standard coefficients.

lm(y~0+x)

## 
## Call:
## lm(formula = y ~ 0 + x)
## 
## Coefficients:
##    xage  xweight  
## 0.07961  0.08936

Then, plot these two output by using MSE.

output_1=sqrt(rowSums((out_1-c(0.07961,0.08936))^2))
output_2=sqrt(rowSums((out_2-c(0.07961,0.08936))^2))
plot(output_1,col=1,pch=20)
lines(output_2,col=2)
legend("topright", legend=c("Frank Wolfe method", "projected gradient descent"),
       col=c(1,2), lty=1:1, cex=0.8)

We can notice that, projected gradient descent method seems to have a faster.

And we can check running time, projected method has a smaller running time.

time_1

## Time difference of 0.151875 secs

time_2

## Time difference of 0.0999341 secs