The building blocks

The machine learning (ML) methods are the numerical building blocks the selectors use: scalers, ridge and elastic-net regression with cross-validation (the \(\lambda_2\) tuners), the singular value decomposition (SVD) and principal component analysis (PCA) behind the sparse PCA selector, and the hierarchical clustering behind the dependency-aware and the grouped selectors. They are exposed as classes and functions of their own, built on Eigen [1].

set.seed(8)
n <- 120
p <- 30
X <- matrix(rnorm(n * p), n, p)
beta <- c(2, -1.5, 0, 0, 1, rep(0, p - 5))
y <- as.numeric(X %*% beta + rnorm(n))

Scalers

MLZScoreScaler centers and divides by the sample standard deviation, MLLpNormScaler centers and divides by the L1 or L2 norm. Both work in place on a double matrix (or a memory-mapped one) and can undo their transformation.

z <- MLZScoreScaler$new()
z$fit(X)
round(z$get_centers()[1:3], 3)
#> [1] -0.103  0.009 -0.104
round(z$get_scales()[1:3], 3)
#> [1] 1.067 1.008 1.007
X_before <- X + 0
invisible(z$transform_inplace(X))
round(colMeans(X)[1:3], 12)
#> [1] 0 0 0
round(apply(X, 2, sd)[1:3], 12)
#> [1] 1 1 1
invisible(z$inverse_transform_inplace(X))
max(abs(X - X_before))
#> [1] 4.440892e-16

fit_transform_inplace() does both in one call, get_dropped_indices() names constant columns, and save() / load() persist the statistics.

l2 <- MLLpNormScaler$new(norm_type = 2, center = TRUE)
invisible(l2$fit_transform_inplace(X))
sqrt(colSums(X^2))[1:3]
#> [1] 1 1 1
invisible(l2$inverse_transform_inplace(X))

Ridge regression

MLRidge solves \(\min_{\boldsymbol\beta} \|\boldsymbol y - \boldsymbol X \boldsymbol\beta\|_2^2 + \lambda \|\boldsymbol\beta\|_2^2\) for one penalty [2], by the primal or the dual Cholesky path depending on the shape of X. MLRidgeCV chooses the penalty by \(K\)-fold cross-validation over a geometric grid and reports the minimum and the one-standard-error rule.

ridge <- MLRidge$new()
round(ridge$solve(X, y, lambda = 1)[1:5], 3)
#> [1]  2.143 -1.639 -0.159  0.077  0.920
cv <- MLRidgeCV$new()
cv$fit(X, y, n_folds = 5, n_lambda = 50)
c(lambda_min = cv$cv_min(), lambda_1se = cv$cv_1se())
#> lambda_min lambda_1se 
#> 0.08377766 0.25877608

plot(cv$get_lambdas(), cv$get_cv_errors(), log = "x", type = "l", col = "#336C68", lwd = 2,
     xlab = "lambda", ylab = "CV error")
abline(v = c(cv$cv_min(), cv$cv_1se()), lty = c(1, 2), col = "red")

Cross-validation error of the ridge penalty with the minimum and the one-standard-error choice.

MLRidgeCV$new(low_memory = TRUE) runs the dual (kernel) form on column blocks of X read in place, with memory \(O(n^2)\) instead of \(O(np)\), the route for \(p \gg n\).

The elastic net

MLElasticNet fits the elastic-net path [3] by coordinate descent over a glmnet-style grid [4]. alpha = 1 is the lasso, alpha = 0 the ridge. The coefficients are reported on the original scale.

en <- MLElasticNet$new()
en$fit(X, y, alpha = 0.5, n_lambda = 50)
dim(en$get_coef())
#> [1] 30 39
length(en$get_lambdas())
#> [1] 39
en$converged()
#> [1] TRUE
B <- en$get_coef()               # the path can stop before n_lambda points
round(B[1:5, ncol(B)], 3)
#> [1]  2.159 -1.650 -0.167  0.071  0.916
round(en$predict(X[1:3, ])[, ncol(B)], 3)
#> [1]  3.120  3.338 -1.793

MLElasticNetCV selects the penalty by cross-validation on a shared grid. The strong rules [5] speed up the sweeps.

encv <- MLElasticNetCV$new()
encv$fit(X, y, alpha = 0.5, n_folds = 5, n_lambda = 50)
c(lambda_min = encv$cv_min(), lambda_1se = encv$cv_1se())
#> lambda_min lambda_1se 
#>  0.2125322  0.4507704

Both classes accept low_memory = TRUE and then run on X itself, a memory-mapped matrix included, without an allocation of the size of X.

SVD and PCA

MLSVD$new()$compute(X, M) returns the top-\(M\) truncated SVD, by a direct, a Gram or a randomized path depending on the shape of X.

svd_res <- MLSVD$new()$compute(X, M = 3)
names(svd_res)
#> [1] "U" "S" "V"
dim(svd_res$U)
#> [1] 120   3
round(svd_res$S, 3)
#> [1] 16.113 14.803 14.607

MLPCA fits the top-\(M\) components. It centers and normalizes X in place at construction, the same rule as the T-solvers, and restore() undoes it.

X_before <- X + 0
pca <- MLPCA$new(X, M = 3, center = TRUE, normalize = TRUE)
pca$fit()
dim(pca$get_scores())
#> [1] 120   3
dim(pca$get_loadings())
#> [1] 30  3
round(pca$get_explained_variance(), 3)
#> [1] 0.016 0.015 0.014
round(pca$transform(X_before[1:2, ]), 3)
#>        [,1]   [,2]   [,3]
#> [1,]  0.056  0.143 -0.064
#> [2,] -0.140 -0.067  0.073
pca$restore()
max(abs(X - X_before))
#> [1] 4.440892e-16

Hierarchical clustering

ml_agglomerative_cluster() runs agglomerative clustering on the rows of a matrix with a linkage from ml_linkage_method (Ward, average, complete, single, WPGMA, median, centroid) and a metric from ml_distance_metric (Euclidean, correlation, Manhattan, and two locality-sensitive hashing (LSH) variants of the correlation distance for very wide data). Single linkage runs the SLINK algorithm [6], the other linkages the nearest-neighbour chain [7], [8] with the Lance-Williams update [9]. The result is a linkage matrix, and ml_cut_tree() cuts it into a given number of clusters.

set.seed(9)
centers <- rep(c(-3, 0, 3), each = 20)
Z <- centers + matrix(rnorm(60 * 4), 60, 4)
linkage <- ml_agglomerative_cluster(Z, method = ml_linkage_method$Ward,
                                    metric = ml_distance_metric$Euclidean)
dim(linkage)
#> [1] 59  4
labels <- ml_cut_tree(linkage, num_orig_objs = 60, num_clusters = 3)
table(labels, rep(1:3, each = 20))
#>       
#> labels  1  2  3
#>      1 20  0  0
#>      2  0  1 20
#>      3  0 19  0

To cluster variables instead of observations, cluster the transposed design with the correlation distance, which is what the grouped and the dependency-aware selectors do.

link_vars <- ml_agglomerative_cluster(t(X), method = ml_linkage_method$Single,
                                      metric = ml_distance_metric$Correlation)
head(ml_cut_tree(link_vars, num_orig_objs = p, num_clusters = 5))
#> [1] 1 1 1 2 1 1

use_mmap = TRUE keeps the intermediate distances in memory-mapped files, and a memory-mapped data matrix is clustered in place.