Regularising high-variance models with synthetic data

Flexible learners such as regression trees and nearest neighbours have low bias but high variance on small samples: small changes in the training data produce very different models. This vignette shows that training on real plus synthetic data from a fitted synthesizer works as a regulariser for such models: it trades a little bias for a large reduction in variance. It also shows the limits of the idea, so that it can be used with care.

library(SynthesizerPlus)
library(rpart)
library(ggplot2)
data(Boston, package = "MASS")

rmse <- function(pred, obs) sqrt(mean((pred - obs)^2))
col_real <- "#2a78d6"
col_aug <- "#eb6834"

The idea

Empirical risk minimisation fits a model to the \(n\) training points only. Vicinal risk minimisation (Chapelle et al., 2001) replaces each training point by a distribution around it and fits the model to that smoothed distribution. A synthesizer provides exactly such a smoothed distribution, and it respects the data: marginal distributions, bounds, discreteness, categories and the dependence between variables. Training on many synthetic records therefore penalises models that react to the idiosyncrasies of individual training points, which is what regularisation does.

Two consequences follow, and the experiments below confirm both:

Throughout, the setting is a small training set of 60 records from the Boston housing data; the target is the median house value medv. The synthesizer is always fitted on the training records only, and models are evaluated on held-out real records.

# four learners with very different capacity
deep_tree <- function(d, te) {
  predict(rpart(medv ~ ., d, control = rpart.control(cp = 0, minsplit = 2,
                                                     minbucket = 1, xval = 0)), te)
}
pruned_tree <- function(d, te) {     # complexity chosen by 10-fold cross-validation
  f <- rpart(medv ~ ., d, control = rpart.control(cp = 0, minsplit = 2,
                                                  minbucket = 1, xval = 10))
  cp <- f$cptable[which.min(f$cptable[, "xerror"]), "CP"]
  predict(prune(f, cp = cp), te)
}
knn <- function(d, te, k) {
  X <- scale(as.matrix(d[, names(d) != "medv"]))
  Xt <- scale(as.matrix(te[, names(d) != "medv"]),
              attr(X, "scaled:center"), attr(X, "scaled:scale"))
  D <- outer(rowSums(Xt^2), rowSums(X^2), "+") - 2 * Xt %*% t(X)
  apply(D, 1, function(dd) mean(d$medv[order(dd)[seq_len(k)]]))
}
learners <- list(
  "fully grown tree" = deep_tree,
  "CV-pruned tree" = pruned_tree,
  "1-nearest neighbour" = function(d, te) knn(d, te, 1),
  "5-nearest neighbours" = function(d, te) knn(d, te, 5)
)

How much synthetic data?

For several random training sets we add 0 to 2000 synthetic records and record the test error of each learner.

sizes <- c(0, 100, 300, 1000, 2000)
curve <- do.call(rbind, lapply(1:8, function(r) {
  set.seed(r)
  idx <- sample(nrow(Boston), 60)
  train <- Boston[idx, ]
  test <- Boston[-idx, ]
  syn <- generate(fit_synthesizer(train), n = max(sizes), seed = r)
  do.call(rbind, lapply(sizes, function(s) {
    d <- rbind(train, syn[seq_len(s), ])
    data.frame(rep = r, synthetic = s, learner = names(learners),
               rmse = vapply(learners, function(f) rmse(f(d, test), test$medv), numeric(1)))
  }))
}))
curve$learner <- factor(curve$learner, levels = names(learners))
avg <- aggregate(rmse ~ synthetic + learner, curve, mean)
avg$rmse <- round(avg$rmse, 2)
reshape(avg, idvar = "learner", timevar = "synthetic", direction = "wide")
#>                 learner rmse.0 rmse.100 rmse.300 rmse.1000 rmse.2000
#> 1      fully grown tree   6.22     6.35     6.42      6.58      6.76
#> 6        CV-pruned tree   6.24     6.17     6.02      5.68      5.56
#> 11  1-nearest neighbour   6.51     6.72     6.78      7.04      7.05
#> 16 5-nearest neighbours   6.22     6.26     6.13      6.12      6.09
ggplot(avg, aes(synthetic, rmse)) +
  geom_line(colour = col_aug, linewidth = 0.8) +
  geom_point(colour = col_aug, size = 2.2) +
  facet_wrap(~learner, nrow = 1) +
  labs(x = "synthetic records added to 60 real ones", y = "test RMSE",
       title = "Augmentation helps learners that control their complexity") +
  theme_minimal(base_size = 10) +
  theme(panel.grid.minor = element_blank())

The two learners that can fit every single point — the fully grown tree and 1-nearest neighbour — get worse: they reproduce the synthetic records, including the noise that the synthesizer adds. The two learners that limit their complexity — the tree pruned by cross-validation and 5-nearest neighbours — get better as synthetic data are added, with most of the gain reached after one or two thousand records.

Augmentation is not a substitute for the usual complexity control; it works together with it. With more (synthetic) data, cross-validation can afford a larger tree whose predictions are still stable.

Where the gain comes from: bias and variance

To separate bias from variance, we fix a test set and draw many different training sets of 60 records from the remaining data. For every test point we record the mean and the spread of the predictions over the training sets.

set.seed(42)
test_idx <- sample(nrow(Boston), 200)
test <- Boston[test_idx, ]
pool <- Boston[-test_idx, ]

R <- 20
pred_real <- pred_aug <- matrix(NA_real_, R, nrow(test))
for (r in seq_len(R)) {
  train <- pool[sample(nrow(pool), 60), ]
  pred_real[r, ] <- pruned_tree(train, test)
  pred_aug[r, ] <- pruned_tree(augment_data(train, n = 2000, seed = r), test)
}
decompose <- function(P) {
  c(bias2 = mean((colMeans(P) - test$medv)^2),
    variance = mean(apply(P, 2, var)),
    mse = mean(sweep(P, 2, test$medv)^2))
}
bv <- rbind(real = decompose(pred_real), augmented = decompose(pred_aug))
round(bv, 1)
#>           bias2 variance  mse
#> real       23.4     21.8 44.1
#> augmented  32.4      7.6 39.6
bv_long <- data.frame(
  training = factor(rep(c("real only", "real + synthetic"), 3),
                    levels = c("real only", "real + synthetic")),
  component = factor(rep(c("squared bias", "variance", "mean squared error"), each = 2),
                     levels = c("squared bias", "variance", "mean squared error")),
  value = as.vector(bv)
)
ggplot(bv_long, aes(component, value, fill = training)) +
  geom_col(position = position_dodge2(padding = 0.15), width = 0.7) +
  scale_fill_manual(values = c(col_real, col_aug), name = NULL) +
  labs(x = NULL, y = NULL, title = "CV-pruned tree trained on 60 records") +
  theme_minimal(base_size = 10) +
  theme(legend.position = "top", panel.grid.minor = element_blank(),
        panel.grid.major.x = element_blank())

The variance of the predictions drops by roughly two thirds, while the squared bias increases noticeably: the synthesizer’s smoothed distribution is not the true one. Because the variance reduction is larger than the bias increase, the total error falls: a classical regularisation trade-off.

Seeing the variance

The effect is easy to see with a single predictor. Below, a CV-pruned tree predicts medv from lstat (share of lower-status population), trained on 15 different samples of 60 records. Each line is one fitted model.

grid <- data.frame(lstat = seq(2, 35, length.out = 200))
fits <- do.call(rbind, lapply(1:15, function(r) {
  set.seed(100 + r)
  train <- pool[sample(nrow(pool), 60), c("lstat", "medv")]
  aug <- augment_data(train, n = 2000, seed = r)
  rbind(data.frame(grid, medv = pruned_tree(train, grid), run = r, training = "real only"),
        data.frame(grid, medv = pruned_tree(aug, grid), run = r, training = "real + synthetic"))
}))
fits$training <- factor(fits$training, levels = c("real only", "real + synthetic"))
ggplot(fits, aes(lstat, medv, group = run, colour = training)) +
  geom_point(data = Boston, aes(lstat, medv), inherit.aes = FALSE,
             colour = "grey80", size = 0.6) +
  geom_step(alpha = 0.6, linewidth = 0.5) +
  facet_wrap(~training) +
  scale_colour_manual(values = c(col_real, col_aug), guide = "none") +
  labs(title = "15 trees, each trained on a different sample of 60 records") +
  theme_minimal(base_size = 10) +
  theme(panel.grid.minor = element_blank())

Trained on the real records alone, the trees disagree strongly, and many are coarse step functions. With synthetic records, the trees are larger and smoother and they agree much more with each other.

Compared with bagging

Bagging (bootstrap aggregation) is the classical way to reduce the variance of trees. How does augmentation compare, and can the two be combined?

bagged_trees <- function(d, te, B = 25) {
  rowMeans(sapply(seq_len(B), function(b) {
    db <- d[sample(nrow(d), replace = TRUE), ]
    predict(rpart(medv ~ ., db, control = rpart.control(cp = 0, minsplit = 5, xval = 0)), te)
  }))
}
cmp <- t(sapply(1:8, function(r) {
  set.seed(r)
  idx <- sample(nrow(Boston), 60)
  train <- Boston[idx, ]
  test <- Boston[-idx, ]
  aug <- augment_data(train, n = 2000, seed = r)
  c(`pruned tree` = rmse(pruned_tree(train, test), test$medv),
    `pruned tree + synthetic` = rmse(pruned_tree(aug, test), test$medv),
    `bagged trees` = rmse(bagged_trees(train, test), test$medv),
    `bagged trees + synthetic` = rmse(bagged_trees(aug, test), test$medv))
}))
round(colMeans(cmp), 2)
#>              pruned tree  pruned tree + synthetic             bagged trees 
#>                     6.12                     5.65                     5.15 
#> bagged trees + synthetic 
#>                     5.11

On this problem bagging alone is more effective than augmenting a single pruned tree, and combining both gives a further small improvement. Synthetic data are therefore best seen as one more regulariser in the toolbox — one that is particularly useful when the model class cannot easily be ensembled, or when the training pipeline expects a single, larger data set (e.g. a neural network).

Tuning the amount of synthetic data

The amount of synthetic data (and the synthesizer settings, such as by, copula or knots) are tuning parameters. Choose them by cross-validation inside the training data, refitting the synthesizer in every fold so that the validation fold is never seen by the synthesizer:

cv_augmentation <- function(train, sizes, learner, folds = 5, seed = 1) {
  set.seed(seed)
  fold <- sample(rep_len(seq_len(folds), nrow(train)))
  err <- sapply(sizes, function(s) {
    mean(sapply(seq_len(folds), function(k) {
      fit_part <- train[fold != k, ]
      valid <- train[fold == k, ]
      d <- if (s == 0) fit_part else augment_data(fit_part, n = s, seed = k)
      rmse(learner(d, valid), valid$medv)
    }))
  })
  data.frame(synthetic = sizes, cv_rmse = err)
}
set.seed(7)
train <- Boston[sample(nrow(Boston), 60), ]
cv_augmentation(train, sizes = c(0, 500, 2000), learner = pruned_tree)
#>   synthetic  cv_rmse
#> 1         0 6.880105
#> 2       500 5.452948
#> 3      2000 5.262143

Summary

References

Chapelle, O., Weston, J., Bottou, L. and Vapnik, V. (2001). Vicinal risk minimization. In Advances in Neural Information Processing Systems 13, 416–422.

Breiman, L. (1996). Bagging predictors. Machine Learning, 24(2), 123–140.