LogisticSlope <- function(X, y, lambda, options=list())
{
# -------------------------------------------------------------
# Start times
# -------------------------------------------------------------
t0 <- proc.time()[3]
# -------------------------------------------------------------
# Define function for retrieving option fields with defaults
# -------------------------------------------------------------
getDefaultField <- function(options,name,default)
{ if (!is.null(options[[name]]))
{ return(options[[name]]) }
else
{ return(default) }
}
# -------------------------------------------------------------
# Parse parameters
# -------------------------------------------------------------
iterations <- getDefaultField(options,"iterations", 10000)
verbosity <- getDefaultField(options,"verbosity" , 1)
optimIter <- getDefaultField(options,"optimIter" , 1)
tolInfeas <- getDefaultField(options,"tolInfeas" , 1e-6)
tolRelGap <- getDefaultField(options,"tolRelGap" , 1e-6)
wInit <- getDefaultField(options,"wInit" , vector())
# -------------------------------------------------------------
# Ensure that lambda is non-increasing
# -------------------------------------------------------------
n <- length(lambda)
matrix(lambda,c(1,n))
if ((n > 1) && any(lambda[2:n] > lambda[1:n-1]))
{ stop("Lambda must be non-increasing.")
}
# -------------------------------------------------------------
# Initialize
# -------------------------------------------------------------
# Get problem dimension
n <- ncol(X)
# Get initial lower bound on the Lipschitz constant
rndKind <- RNGkind()
rndState <- tryCatch({.Random.seed},error=function(m) {} )
set.seed(0,kind="Mersenne-Twister",normal.kind="Inversion")
w <- matrix(rnorm(n),c(n,1))
w <- w / sqrt(sum(w^2))
w <- t(X) %*% (X %*% w)
L <- sqrt(sum(w^2))/4
if (!is.null(rndState))
set.seed(rndState[-1],kind=rndKind[1],normal.kind=rndKind[2])
# Set constants
STATUS_RUNNING <- 0
STATUS_OPTIMAL <- 1
STATUS_ITERATIONS <- 2
STATUS_MSG <- c('Optimal','Iteration limit reached')
# Initialize parameters and iterates
if (length(wInit) == 0) wInit <- matrix(0,n,1)
t <- 1
eta <- 2
lambda <- matrix(lambda,nrow=length(lambda))
y <- matrix(y,nrow=length(y))
Y <- diag(as.vector(y), length(y))
YX <- Y %*% X
w <- wInit
w0 <- 0
v <- w
iter <- 0
status <- STATUS_RUNNING
# Deal with Lasso case
modeLasso <- (length(lambda) == 1)
if (modeLasso)
proxFunction <- function(x,lambda) { return(sign(x) * pmax(abs(x) - lambda,0)) }
else
proxFunction <- function(v1,v2) { return(prox_sorted_L1(v1,v2)) }
if (verbosity > 0)
{ printf <- function(...) invisible(cat(sprintf(...)))
printf('%5s %9s %9s %9s\n','Iter','Gap','Infeas.','Rel. gap')
}
# -------------------------------------------------------------
# Main loop
# -------------------------------------------------------------
while (TRUE)
{
# Compute the gradient at f(v)
r <- 1 / (1 + exp(YX %*% v))
g <- -as.double(crossprod(YX, r))
log1mr <- log(1-r)
f <- -sum(log1mr)
# Increment iteration count
iter <- iter + 1
# Check optimality conditions
if ((iter %% optimIter) == 0)
{
# Compute 'dual', check infeasibility and gap
if (modeLasso) {
infeas <- max(abs(g)-lambda,0)
objPrimal <- f + lambda * norm(v,'1')
} else {
gs <- sort(abs(g), decreasing=TRUE)
vs <- sort(abs(v), decreasing=TRUE)
infeas <- max(max(cumsum(gs-lambda)),0)
# Compute primal and dual objective
objPrimal <- f + as.double(crossprod(lambda,vs))
}
objDual <- as.double(crossprod(r-1, log1mr)) - as.double(crossprod(r, log(r)))
# Format string
if (verbosity > 0)
str <- sprintf(' %9.2e %9.2e %9.2e',objPrimal - objDual, infeas/lambda[[1]], abs(objPrimal - objDual) / max(1,objPrimal))
# Check primal-dual gap
if ((abs(objPrimal - objDual)/max(1,objPrimal) < tolRelGap) &&
(infeas < tolInfeas * lambda[[1]]))
status <- STATUS_OPTIMAL
}
else
{ str <- ''
}
if (verbosity > 0)
{ if ((verbosity == 2) ||
((verbosity == 1) && ((iter %% optimIter) == 0)))
{ printf('%5d %s\n', iter, str)
}
}
# Stopping criteria
if ((status == 0) && (iter >= iterations))
status <- STATUS_ITERATIONS
if (status != 0)
{ if (verbosity > 0)
printf('Exiting with status %d -- %s\n', status, STATUS_MSG[[status]])
break
}
# Keep copies of previous values
wPrev <- w
fPrev <- f
tPrev <- t
# Lipschitz search
while (TRUE)
{ # Compute prox mapping
w <- proxFunction(v - (1/L)*g, lambda/L)
d <- w - v
r <- 1/(1 + exp(YX %*% w))
f <- -sum(log(1-r))
q <- fPrev + as.double(crossprod(d,g)) + (L/2)*as.double(crossprod(d))
if (q >= f*(1-1e-12))
break
else
L <- L * eta
} # Lipschitz search
# Update
t <- (1 + sqrt(1 + 4*t^2)) / 2
v <- w + ((tPrev - 1) / t) * (w - wPrev)
} # While (TRUE)
# Information structure
info <- c()
info$runtime <- proc.time()[3] - t0
info$objPrimal <- objPrimal
info$objDual <- objDual
info$infeas <- infeas
info$status <- status
info$w <- v
info$L <- L
return(info)
} # Function Adlas
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.