# projection on l1 ball with weigths pib1w <- function(x,w,z=1){ if (sum(abs(x))>z & z != Inf){ v <- abs(x* w) u <- order(-v) sx <- cumsum(abs(x)[u]) sw <- cumsum(1/w[u]) rho <- max(which(v[u] > (sx-z) / sw)) theta <- (sx[rho] -z) / sw[rho] x<-sign(x)*pmax(abs(x) - theta/w, 0)} return(x) } # Stochastic Mirror Descent with projection on the l1 ball SMDproj <- function(a, b, init, iters = length(b), cost, instgrad, lambda,z =Inf) { ind<-sample(1:length(b),iters,replace=TRUE) a<-a[ind,] b<-b[ind] param <- data.frame(matrix(nrow = iters, ncol = length(init) + 1)) colnames(param) <- c(colnames(a), "Loss") x <- c(init) m <- x param[1, ] <- c(m, cost( m, a, b,lambda)) theta <- c(init) for (i in 2:iters) { eta <- 1/(lambda*i) # 1/i play with the learning rate theta <- theta - eta * instgrad(x, a[i,], b[i], lambda) x <- pib1(theta,z) m <- ((i-1)*m + x)/i param[i, ] <- c(m ,cost( m, a, b,lambda) ) } cat("Final cost: ", sprintf("%10.07f", param[nrow(param), ncol(param)]), "\n", sep = "") cat("Parameters:", as.numeric(param[nrow(param), 1:length(init) ]), sep = " ") param <- cbind(Iteration = 1:nrow(param), param) return(param) } # AdaGrad with projection Adaproj <- function(a, b, init, iters = length(b), cost, instgrad, lambda,z=Inf) { ind<-sample(1:length(b),iters,replace=TRUE) a<-a[ind,] b<-b[ind] param <- data.frame(matrix(nrow = iters, ncol = length(init) + 1)) colnames(param) <- c(colnames(a), "Loss") s <- rep(1/(4*length(init)),length(init)) eta <- 1/(lambda * sqrt(iters)) x <- c(init) m <- x param[1, ] <- c(m, cost(m, a, b,lambda)) for (i in 2:iters) { s <- s + instgrad(x, a[i,], b[i], lambda)^2 y <- x - eta * 1/sqrt(s) * instgrad(x, a[i,], b[i], lambda) x <- pib1w(y,sqrt(s),z) m <- ((i-1)*m+x)/i param[i,] <- c(m, cost(m, a, b,lambda)) } cat("Final cost: ", sprintf("%10.07f", param[nrow(param), ncol(param)]), "\n", sep = "") cat("Parameters:", as.numeric(param[nrow(param), 1:length(init) ]), sep = " ") param <- cbind(Iteration = 1:nrow(param), param) return(param) } # Generate the train sample atrain <- data.frame(intercept= rep(1,n = 100), a = rnorm(n = 100) * 5, b = rnorm(n = 100) * 3 + 1, c = rnorm(n = 100) * 2 + 2, d= rnorm(n = 100) * 5) # Add column names and intercept colnames(atrain) <- c("intercept", "a", "b", "c" ,"d") atrain <- as.matrix(atrain) # Setting up a perfect relationship (linear problem) with +1 or -1 for classes xstar <- c(-20,2,3,4,0) btrain <- 2 * as.integer(atrain%*%xstar > 0) - 1 #Comparison SMD proj and Ada on the same l1 ball and same seed start_time <- Sys.time() set.seed(100) paramAdaproj <- Adaproj(a = atrain, b = btrain, init = rep(0, dim(atrain)[2]), iters = 1000, # Play with the number of iterations cost = hingereg, instgrad = instgradreg, lambda = 1, # Play with the regularization parameter z=1) end_time <- Sys.time() end_time - start_time start_time <- Sys.time() set.seed(100) # Same seed! paramSMDproj <- SMDproj(a = atrain, b = btrain, init = rep(0, 5), iters = 1000, # Play with the number of iterations cost = hingereg, instgrad = instgradreg, lambda = 1, # Play with the regularization parameter z=1) end_time <- Sys.time() end_time - start_time # Different trajectories plot.ts(paramAdaproj[, 2:(ncol(paramAdaproj) - 1)]) plot.ts(paramSMDproj[, 2:(ncol(paramSMDproj) - 1)]) mean(paramAdaproj$Loss-paramSGDproj$Loss) # test # Explanatory variables and intercept in the test set atest <- data.frame(intercept= rep(1,n = 1000), a = rnorm(n = 1000) * 5, b = rnorm(n = 1000) * 3 + 1, c = rnorm(n = 1000) * 2 + 2, d = rnorm(n = 1000) * 5) atest <- as.matrix(atest) # test labels with +1 or -1 for classes btest <- 2 * as.integer(atest %*% xstar > 0) - 1 # How to analyse the convergence using a test set rateAdaproj <- colMeans(btest*( atest %*% t(as.matrix(paramAdaproj[, 2:(ncol(paramAdaproj) - 1)])))>0) rateSMDproj <- colMeans(btest*( atest %*% t(as.matrix(paramSMDproj[, 2:(ncol(paramSMDproj) - 1)])))>0) mean(rateAdaproj-rateSMDproj) plot(log(paramAdaproj$Iteration),log(1-rateAdaproj),type ="l",col="blue") lines(log(paramSMDproj$Iteration),log(1-rateSMDproj))