project-4-mine.R

data(zip.train, package="ElemStatLearn")
table(zip.train[, 1])
table(zip.train[, 2])
all.y.vec <- zip.train[, 1]
is.01 <- all.y.vec %in% c(0,1)
table(is.01)
y.vec <- all.y.vec[is.01]
X.mat <- zip.train[is.01, -1]
str(y.vec)
str(X.mat)

set.seed(1)
n.folds <- 5
fold.vec <- sample(rep(1:n.folds, l=nrow(X.mat)))
validation.fold <- 1
is.train <- fold.vec != validation.fold
is.validation <- fold.vec == validation.fold

X.train <- X.mat[is.train,]
y.train <- y.vec[is.train]

X.sc <- scale(X.train)
y.tilde <- ifelse(y.train==1, 1, -1)
table(y.tilde, y.train)
X.filtered <- X.sc[, attr(X.sc, "scaled:scale") != 0]

## NB: first element is bias/intercept.
X.int <- cbind(1, X.filtered)
w.vec <- rep(0, l=ncol(X.int))
sigmoid <- function(z){
  1/(1+exp(-z))
}
step.size <- 0.1
lambda <- 0.5
ppart <- function(x){
  ifelse(x<0, 0, x)
}
soft <- function(x, l){
  sign(x) * ppart(abs(x)-l)
}
prox <- function(x, l){
  c(x[1], soft(x[-1], l))
}
wstep <- function(size){
  prox(w.vec+size*d.vec, lambda*size)
}
log.loss <- function(pred){
  log(1+exp(-y.tilde * pred))
}
cstep <- function(size){
  mean(log.loss(X.int %*% wstep(size)))
}
subdiff.crit <- function(w,d){
  ifelse(
    w==0,
    ppart(abs(d)-lambda),
    abs(d-sign(w)*lambda))
}

pred.vec <- X.int %*% w.vec
prob.vec <- sigmoid(-pred.vec * y.tilde)
grad.vec <- -t(X.int) %*% (y.tilde * prob.vec) / nrow(X.int)
d.vec <- -grad.vec
curve(sapply(x, cstep), 0, 10)
step.fac <- 2
while(cstep(step.size/step.fac) < cstep(step.size)){
  step.size <- step.size/step.fac
}
while(cstep(step.size*step.fac) < cstep(step.size)){
  step.size <- step.size*step.fac
}
points(step.size, cstep(step.size))
crit.vec <- c(abs(d.vec[1]), subdiff.crit(w.vec[-1], d.vec[-1]))
cat(sprintf("crit=%f cost=%f step=%f\n", max(crit.vec), cstep(0), step.size))
w.vec <- wstep(step.size)