stochastic-gradient-maxLik.R
In maxLik: Maximum Likelihood Estimation and Related Tools

### R code from vignette source 'stochastic-gradient-maxLik.Rnw'

###################################################
### code chunk number 1: foo
###################################################
options(keep.source = TRUE, width = 60,
        try.outFile=stdout()  # make try to produce error messages
        )
foo <- packageDescription("maxLik")


###################################################
### code chunk number 2: stochastic-gradient-maxLik.Rnw:163-164 (eval = FALSE)
###################################################
## maxSGA(fn, grad, start, nObs, control)


###################################################
### code chunk number 3: stochastic-gradient-maxLik.Rnw:375-379
###################################################
i <- which(names(MASS::Boston) == "medv")
X <- as.matrix(MASS::Boston[,-i])
X <- cbind("const"=1, X)  # add constant
y <- MASS::Boston[,i]


###################################################
### code chunk number 4: stochastic-gradient-maxLik.Rnw:384-385
###################################################
colMeans(X)


###################################################
### code chunk number 5: stochastic-gradient-maxLik.Rnw:392-393
###################################################
eigenvals <- eigen(crossprod(X))$values


###################################################
### code chunk number 6: stochastic-gradient-maxLik.Rnw:404-407
###################################################
betaX <- solve(crossprod(X)) %*% crossprod(X, y)
betaX <- drop(betaX)  # matrix to vector
betaX


###################################################
### code chunk number 7: stochastic-gradient-maxLik.Rnw:431-436
###################################################
gradloss <- function(theta, index)  {
   e <- y[index] - X[index,,drop=FALSE] %*% theta
   g <- t(e) %*% X[index,,drop=FALSE]
   2*g/length(index)
}


###################################################
### code chunk number 8: gradonly
###################################################
library(maxLik)
set.seed(3)
start <- setNames(rnorm(ncol(X), sd=0.1), colnames(X))
                           # add names for better reference
res <- try(maxSGA(grad=gradloss,
           start=start,
           nObs=nrow(X),
           control=list(iterlim=1000)
           )
    )


###################################################
### code chunk number 9: stochastic-gradient-maxLik.Rnw:476-483
###################################################
res <- maxSGA(grad=gradloss,
              start=start,
              nObs=nrow(X),
              control=list(iterlim=1000,
                           SG_clip=1e4)  # limit ||g|| <= 100
              )
summary(res)


###################################################
### code chunk number 10: stochastic-gradient-maxLik.Rnw:495-499
###################################################
loss <- function(theta, index) {
   e <- y[index] - X[index,] %*% theta
   -crossprod(e)/length(index)
}


###################################################
### code chunk number 11: stochastic-gradient-maxLik.Rnw:514-538
###################################################
res <- maxSGA(loss, gradloss,
              start=start,
              nObs=nrow(X),
              control=list(iterlim=1000,
                           # will misbehave with larger numbers
                           SG_clip=1e4,
                           SG_learningRate=0.001,
                           storeParameters=TRUE,
                           storeValues=TRUE
                           )  
              )
par <- storedParameters(res)
val <- storedValues(res)
par(mfrow=c(1,2))
plot(par[,1], par[,2], type="b", pch=".",
     xlab=names(start)[1], ylab=names(start)[2], main="Parameters")
## add some arrows to see which way the parameters move
iB <- c(40, nrow(par)/2, nrow(par))
iA <- iB - 10
arrows(par[iA,1], par[iA,2], par[iB,1], par[iB,2], length=0.1)
##
plot(seq(length=length(val))-1, -val, type="l",
     xlab="epoch", ylab="MSE", main="Loss",
     log="y")


###################################################
### code chunk number 12: stochastic-gradient-maxLik.Rnw:565-570
###################################################
i <- sample(nrow(X), 0.8*nrow(X))  # training indices, 80% of data
Xt <- X[i,]  # training data
yt <- y[i]
Xv <- X[-i,]  # validation data
yv <- y[-i]


###################################################
### code chunk number 13: stochastic-gradient-maxLik.Rnw:575-584
###################################################
gradloss <- function(theta, index)  {
   e <- yt[index] - Xt[index,,drop=FALSE] %*% theta
   g <- -2*t(e) %*% Xt[index,,drop=FALSE]
   -g/length(index)
}
loss <- function(theta, index) {
   e <- yv - Xv %*% theta
   -crossprod(e)/length(yv)
}


###################################################
### code chunk number 14: batch1
###################################################
res <- maxSGA(loss, gradloss,
              start=start,
              nObs=nrow(Xt),  # note: only training data now
              control=list(iterlim=100,
                           SG_batchSize=1,
                           SG_learningRate=1e-5,
                           SG_clip=1e4,
                           storeParameters=TRUE,
                           storeValues=TRUE
                           )  
              )
par <- storedParameters(res)
val <- storedValues(res)
par(mfrow=c(1,2))
plot(par[,1], par[,2], type="b", pch=".",
     xlab=names(start)[1], ylab=names(start)[2], main="Parameters")
iB <- c(40, nrow(par)/2, nrow(par))
iA <- iB - 1
arrows(par[iA,1], par[iA,2], par[iB,1], par[iB,2], length=0.1)
plot(seq(length=length(val))-1, -val, type="l",
     xlab="epoch", ylab="MSE", main="Loss",
     log="y")


###################################################
### code chunk number 15: momentum
###################################################
res <- maxSGA(loss, gradloss,
              start=start,
              nObs=nrow(Xt),
              control=list(iterlim=100,
                           SG_batchSize=1,
                           SG_learningRate=1e-6,
                           SG_clip=1e4,
                           SGA_momentum = 0.99,
                           storeParameters=TRUE,
                           storeValues=TRUE
                           )  
              )
par <- storedParameters(res)
val <- storedValues(res)
par(mfrow=c(1,2))
plot(par[,1], par[,2], type="b", pch=".",
     xlab=names(start)[1], ylab=names(start)[2], main="Parameters")
iB <- c(40, nrow(par)/2, nrow(par))
iA <- iB - 1
arrows(par[iA,1], par[iA,2], par[iB,1], par[iB,2], length=0.1)
plot(seq(length=length(val))-1, -val, type="l",
     xlab="epoch", ylab="MSE", main="Loss",
     log="y")


###################################################
### code chunk number 16: Adam
###################################################
res <- maxAdam(loss, gradloss,
              start=start,
              nObs=nrow(Xt),
              control=list(iterlim=100,
                           SG_batchSize=1,
                           SG_learningRate=1e-6,
                           SG_clip=1e4,
                           storeParameters=TRUE,
                           storeValues=TRUE
                           )  
              )
par <- storedParameters(res)
val <- storedValues(res)
par(mfrow=c(1,2))
plot(par[,1], par[,2], type="b", pch=".",
     xlab=names(start)[1], ylab=names(start)[2], main="Parameters")
iB <- c(40, nrow(par)/2, nrow(par))
iA <- iB - 1
arrows(par[iA,1], par[iA,2], par[iB,1], par[iB,2], length=0.1)
plot(seq(length=length(val))-1, -val, type="l",
     xlab="epoch", ylab="MSE", main="Loss",
     log="y")


###################################################
### code chunk number 17: SANN
###################################################
val <- NULL
# loop over batch sizes
for(B in c(1,10,100)) {
   res <- maxAdam(loss, gradloss,
                  start=start,
                  nObs=nrow(Xt),
                  control=list(iterlim=200,
                               SG_batchSize=1,
                               SG_learningRate=1e-6,
                               SG_clip=1e4,
                               SG_patience=5,
                           # worse value allowed only 5 times
                               storeValues=TRUE
                               )  
                  )
   cat("Batch size", B, ",", nIter(res),
       "epochs, function value", maxValue(res), "\n")
   val <- c(val, na.omit(storedValues(res)))
   start <- coef(res)
}
plot(seq(length=length(val))-1, -val, type="l",
     xlab="epoch", ylab="MSE", main="Loss",
     log="y")
summary(res)