Calculating mse for multiple dimensionality reduction technique

Question

I'm trying to find a metric to compare multiple dimensionality reduction techniques similar to what was done in this blog post pca-vs-autoencoders-for-dimensionality-reduction...

Specifically this part of the comparison

# pCA reconstruction
pca.recon <- function(pca, x, k){
  mu <- matrix(rep(pca$center, nrow(pca$x)), nrow = nrow(pca$x), byrow = T)
  recon <- pca$x[,1:k] %*% t(pca$rotation[,1:k]) + mu
  mse <- mean((recon - x)^2)
  return(list(x = recon, mse = mse))
}

xhat <- rep(NA, 10)
for(k in 1:10){
  xhat[k] <- pca.recon(pca, x_train, k)$mse
}

ae.mse <- rep(NA, 5)
for(k in 1:5){
  modelk <- keras_model_sequential()
  modelk %>%
    layer_dense(units = 6, activation = "tanh", input_shape = ncol(x_train)) %>%
    layer_dense(units = k, activation = "tanh", name = "bottleneck") %>%
    layer_dense(units = 6, activation = "tanh") %>%
    layer_dense(units = ncol(x_train))

  modelk %>% compile(
    loss = "mean_squared_error", 
    optimizer = "adam"
  )

  modelk %>% fit(
    x = x_train, 
    y = x_train, 
    epochs = 5000,
    verbose = 0
  )

  ae.mse[k] <- unname(evaluate(modelk, x_train, x_train))
}

df <- data.frame(k = c(1:10, 1:5), mse = c(xhat, ae.mse), method = c(rep("pca", 10), rep("autoencoder", 5)))
ggplot(df, aes(x = k, y = mse, col = method)) + geom_line()

i will like to add other techniques to the mix such as TSNE from Rtsne package, UMAP from the umap package and IVIS from the ivis package (currently not on CRAN but can be installed like so ->

devtools::install_github("beringresearch/ivis/R-package")
library(ivis)
install_ivis()

The data input and processing for all the techniques are similar but it seems some of them already have mse determination baked into their functions (e.g. autoencoder). I'm wondering if anyone has experience with what i'm trying to do.

Answer 1

Different decomposition methods can be regarded as interchangeable gears in a statistical machine, that serves to be useful to you, the creator.

To pick the best gear, the metric you evaluate by should not necessarily concern the gears, but how well the machine performs overall with each gear inserted respectively.

Disregard the gear specs: You have a couple of gears which all come with their own validation specs from their factories (packages). Those numbers/summeries/specs are likely not what you want. Likely gears will not be provided with the same metrics, so it is hard to do a fair comparison. Also, those metrics will be all about the gears and not about your specific machine anyways. Do not do as the blog suggest and mix machine metrics with the gears in the pca.recon(). Let the gears be gears, and delay the metric evaluation onto the machine level.

Does the gear fit at all?: You need to check for your specific machine, that all your candidate gears will actually fit inside. Your gears for your composition/reconstruction machine, must be able to turn both ways. The t-sne is only designed to turn forward and do decomposition, so it is not possible to do a meaningful evaluation. Likewise for the UMAP. 

Maybe the whole reconstruction loss benchmarking thing was not the actual machine you wanted to use in the first place. Maybe just a side project to pick gears for another machine... 
If your machine is to make beautiful plots then good quantitative benchmarking is hard to come by. If your machine is some initial decomposition mixed with a simple classifier, then the t-sne gear will fit just fine and some prediction accuracy metric is likely useful to pick gears with

.

Interfacing various gears: The gears will actually not fit out-of-the-box into your machine because the size and form is not the same. Each gear needs individual adaptation. You may be tempted to re-fit the machine to the gear and that will do for a few gears. That would be to literal copy-paste your machine code, insert and adapt each gear. A more scalable way is to just interface the gears, such you can leave them in a bag next to the machine and let a robot insert one gear at the time and write you a report. That is a main selling point of frameworks such as sklearn, caret and keras. You could also code it your self. Here is a simple example

:

rm(list=ls())
#some data
X <- iris[,c("Sepal.Length","Sepal.Width","Petal.Length","Petal.Width")]

#my_gear, prcomp wrapped in an interface
#any gear must have the gear(X, N, ...) signature
pca_decompose <- function(X, N=2, ...) {
  
  #implement gear forward (decompose)
  pca <- prcomp(
    X, rank. = N,
    scale = FALSE #must be false, beacuse reconstructor below does not support re-scaling, because I'm lazy.
  )
  
  #implement gear backward (reconstruct)
  reconstruct <- function(Xnew = pca$x) {
    
    # a pca reconstructor implementation similar to function from the blog, pca already in closure
    # I think the blog mistankenly referred to pca$x instead of x sometimes
    pca.recon <- function(x, k){
      x_recon <- x[,1:k] %*% t(pca$rotation[,1:k])
      #slightly more effecient way to reapply center
      for(i in seq_along(pca$center)) x_recon[,i] <- x_recon[,i] + pca$center[i] 
      return(x_recon)
    }
    X_rc <- pca.recon(Xnew, k=N)
    return(X_rc)
  }
  
  #wrap up the interface
  self <- list(
    X_decomposed = pca$x,  # any decomposition must be named X_dc
    reconstruct = reconstruct
  )
  
  class(self) <- c("my_pca","my_universal_gear")
  return(self)
}

#define a machine with the relevant use case
my_machine <- function(gear, data, ...) {
 dc_obj <- gear(data, ...)
 data_rc <- dc_obj$reconstruct(dc_obj$X_decomposed)
}

#define the most useful metric
my_metric <- function(X,Y) {
  # this 'multivariate' mse, is not commonly used I think.
  # but whatever floats the boat
  mean((X-Y)^2) 
}

#define how to evaluate.
#try the gear in the mahine and meassure outcome with metric
my_evaluation <- function(gear, machine, data, metric, ...) {
  data <- as.matrix(data)
  output <- machine(gear,data, ...)
  my_metric(data,output)
}

#useful syntactic sugar
set_params <- function(gear, ...) {
  params = list(...)
  function(...) do.call(gear,c(list(...),params))
}

#evaluate a gear
my_evaluation(
  gear = pca_decompose,
  machine = my_machine,
  data = X,
  #gear params
  N=2
)

#the same as
my_evaluation(
  gear = set_params(pca_decompose,N=2), #nice to preset gear params
  machine = my_machine,
  data = X
)

#define all gears to evaluate
#the gearbag could also in another usecase be a grid search of optimal hyper-parameters
my_gearbag = list(
  pca_dc_N1 = set_params(pca_decompose,N=1),
  pca_dc_N2 = set_params(pca_decompose,N=2),
  pca_dc_N3 = set_params(pca_decompose,N=3),
  pca_dc_N4 = set_params(pca_decompose,N=4)
  #put also autoencoder or what ever in the gearbag
)

my_robot <- function(evaluation, machine, gearbag, data) {
  results <- sapply(
    X = gearbag, #this X is not the data put placeholder for what to iterate
    FUN = evaluation,
    machine = machine,
    data = X
  )
  
  report = list(
    README = "metric results for gears",
    results = results
  )
}


my_report <- my_robot(my_evaluation, my_machine, my_gearbag, data)

print(my_report)

printout

$README
[1] "metric results for gears"

$results
   pca_dc_N1    pca_dc_N2    pca_dc_N3    pca_dc_N4 
8.560431e-02 2.534107e-02 5.919048e-03 1.692109e-31