In tdhock/aum: Area Under Minimum of False Positives and Negatives

Accuracy comparison

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

Load data

data(package="mlbench")
data(Sonar, package="mlbench")
data(DNA, package="mlbench")
data(spam, package="kernlab")
data.list <- list(
  spam=list(
    input.mat=as.matrix(spam[1:57]),
    output.vec=ifelse(spam$type=="spam", 1, -1)),
  Sonar=list(
    input.mat=as.matrix(Sonar[,1:60]),
    output.vec=ifelse(Sonar$Class=="R", 1, -1)),
  DNA=list(
    input.mat=ifelse(as.matrix(DNA[,1:180])==0, 0, 1),
    output.vec=ifelse(DNA$Class=="n", -1, 1)))

Define functions for computing loss and gradient

N <- 10
set.seed(1)
rand.pred.vec <- rnorm(N)
subtrain.output.vec <- rep(c(-1, 1), l=N)
subtrain.diff.count.dt <- aum::aum_diffs_binary(subtrain.output.vec, denominator="count")
subtrain.diff.rate.dt <- aum::aum_diffs_binary(subtrain.output.vec, denominator="rate")
library(data.table)
PairsDT <- function(output.vec){
  is.positive <- output.vec == 1
  data.table(expand.grid(
    positive=which(is.positive),
    negative=which(!is.positive)))
}
subtrain.pairs.dt <- PairsDT(subtrain.output.vec)
margin <- 1
## Note: for efficiency subtrain labels are assumed to be pre-computed
## in the enclosing environment, once before the optimization starts.
equal.class.weights <- function(output.vec){
  otab <- table(output.vec)
  as.numeric(1/otab[paste(output.vec)])
}
subtrain.obs.weights <- equal.class.weights(subtrain.output.vec)
Logistic <- function(pred.vec, output.vec, obs.weights){
  list(
    gradient=-obs.weights*output.vec/(1+exp(output.vec*pred.vec)),
    loss=sum(obs.weights*log(1+exp(-output.vec*pred.vec))))
}
AUM <- function(pred.vec, diff.dt){
  L <- aum::aum(diff.dt, pred.vec)
  d <- L$derivative_mat
  non.diff <- abs(d[,1] - d[,2]) > 1e-6
  if(any(non.diff)){
    cat(sprintf("%d non-diff points\n", sum(non.diff)))
    print(d[non.diff, ])
  }
  with(L, list(gradient=derivative_mat[,1], loss=aum))
}  
loss.list <- list(
  logistic=function(pred.vec, output.vec=subtrain.output.vec, ...){
    Logistic(pred.vec, output.vec, 1/length(pred.vec))
  },
  logistic.weighted=
    function(pred.vec, output.vec=subtrain.output.vec,
             obs.weights=subtrain.obs.weights, ...){
      Logistic(pred.vec, output.vec, obs.weights)
    },
  aum.count=function(pred.vec, diff.count.dt=subtrain.diff.count.dt, ...){
    AUM(pred.vec, diff.count.dt)
  },
  aum.rate=function(pred.vec, diff.rate.dt=subtrain.diff.rate.dt, ...){
    AUM(pred.vec, diff.rate.dt)
  },
  squared.hinge.all.pairs=function(pred.vec, pairs.dt=subtrain.pairs.dt, ...){
    pairs.dt[, diff := pred.vec[positive]-pred.vec[negative]-margin]
    pairs.dt[, diff.clipped := ifelse(diff<0, diff, 0)]
    pairs.tall <- data.table::melt(
      pairs.dt,
      measure.vars=c("positive", "negative"),
      value.name="pred.i",
      variable.name="label")
    ## d/dx (x - y - m)^2 = x - y - m
    ## d/dy (x - y - m)^2 = -(x - y - m)
    pairs.tall[, grad.sign := ifelse(label=="positive", 1, -1)]
    N.pairs <- nrow(pairs.dt)
    grad.dt <- pairs.tall[, .(
      gradient=sum(grad.sign*diff.clipped)
    ), keyby=pred.i]
    list(
      gradient=grad.dt$gradient/N.pairs,
      loss=sum(pairs.dt$diff.clipped^2)/N.pairs)
  }
)
result.list <- list()
for(loss.name in names(loss.list)){
  fun <- loss.list[[loss.name]]
  result.list[[loss.name]] <- fun(rand.pred.vec)
}
str(result.list)
sapply(result.list, "[[", "gradient")

Optimization

out.loss.list <- list()
data.name <- "Sonar"
input.output.list <- data.list[[data.name]]
input.mat <- input.output.list[["input.mat"]]
full.input.mat <- scale(input.mat)
full.output.vec <- input.output.list[["output.vec"]]
stopifnot(full.output.vec %in% c(-1, 1))
set.seed(1)
n.folds <- 2
unique.folds <- 1:n.folds
fold.vec <- sample(rep(unique.folds, l=length(full.output.vec)))
validation.fold <- 1
is.set.list <- list(
  validation=fold.vec == validation.fold,
  subtrain=fold.vec != validation.fold)
set.data.list <- list()
for(set.name in names(is.set.list)){
  is.set <- is.set.list[[set.name]]
  output.vec <- full.output.vec[is.set]
  set.data.list[[set.name]] <- list(
    output.vec=output.vec,
    obs.weights=equal.class.weights(output.vec),
    input.mat=full.input.mat[is.set,],
    diff.rate.dt=aum::aum_diffs_binary(output.vec, denominator="rate"),
    diff.count.dt=aum::aum_diffs_binary(output.vec, denominator="count"),
    pairs.dt=PairsDT(output.vec))
}
X.mat <- set.data.list$subtrain$input.mat
for(loss.name in names(loss.list)){
  loss.grad.fun <- loss.list[[loss.name]]
  step.size <- 0.1
  cat(sprintf("loss=%s\n", loss.name))
  set.seed(1)
  weight.vec <- last.w <- rnorm(ncol(X.mat))
  done <- FALSE
  iteration <- 0
  prev.set.loss.vec <- rep(1e10, 2)
  while(!done){
    iteration <- iteration+1
    loss.for.weight <- function(w, set.data=set.data.list$subtrain){
      pred <- set.data$input.mat %*% w
      set.data$pred.vec <- pred
      out <- do.call(loss.grad.fun, set.data)
      out$pred <- pred
      out
    }
    loss.before.step <- loss.for.weight(weight.vec)
    direction <- -t(X.mat) %*% loss.before.step[["gradient"]]
    loss.for.step <- function(step.size){
      new.weight <- weight.vec + step.size * direction
      out <- loss.for.weight(new.weight)
      out$new.weight <- new.weight
      out$step.size <- step.size
      out
    }
    loss.after.step <- loss.for.step(step.size)
    weight.vec <- loss.after.step[["new.weight"]]
    diff.w <- sum(abs(weight.vec-last.w))
    last.w <- weight.vec
    set.loss.vec <- numeric()
    for(set.name in names(set.data.list)){
      set.data <- set.data.list[[set.name]]
      set.loss <- loss.for.weight(weight.vec, set.data)
      set.loss.vec[[set.name]] <- set.loss[["loss"]]
      roc.df <- WeightedROC::WeightedROC(
        set.loss[["pred"]],
        set.data[["output.vec"]])
      auc <- WeightedROC::WeightedAUC(roc.df)
      out.dt <- data.table(
        loss.name,
        iteration,
        set.name,
        auc,
        loss.value=set.loss$loss)
      for(aum.type in c("count", "rate")){
        diff.name <- paste0("diff.", aum.type, ".dt")
        aum.list <- aum::aum(set.data[[diff.name]], set.loss[["pred"]])
        out.col <- paste0("aum.", aum.type)
        out.dt[[out.col]] <- aum.list[["aum"]]
      }
      out.loss.list[[paste(
        loss.name,
        iteration,
        set.name
      )]] <- out.dt
      if(2000 < iteration || diff.w < 1e-6){
        done <- TRUE
      }
    }#for(set.name
  }#while(!done
}#for(loss.name

out.loss <- do.call(rbind, out.loss.list)
(max.valid.auc <- out.loss[
  set.name=="validation",
  .SD[which.max(auc), .(iteration, auc, set.name)],
  by=.(loss.name)])
library(ggplot2)
ggplot(,aes(
  iteration, auc, color=set.name))+
  geom_line(
    data=out.loss)+
  geom_point(
    data=max.valid.auc)+
  facet_grid(
    . ~ loss.name,
    labeller="label_both",
    scales='free',
    space='fixed')