Subspace ensembles for scalable statistical analyses

## Loading required package: ggplot2

What is so challenging about "Big Data"?

What are the current approaches to large-scale statistical analyses?

Subspace ensembles - a unifying framework for large class of scalable analyses

Aggregation of "weak" estimators

Applications: Medical record linkage and network intrusion detection

My current work and open questions

A data set is "large" if it exceeds 20% of the available RAM for a single machine but can be analyzed on a single machine

• "R is not well-suited for working with data structures larger than 10-20% of a computer's RAM." -The R Installation and Administration Manual
• You can may be able to use standard tools but you need to be careful

A data set is "massive" if it requires more than one machine to process

• Need an approach that will "scale-out"
• Applications using "massive" data platforms is a subset of those who actually need it

s - single machine solution
c - cluster solution

1. Updatable/Downdatable

   • Allows new data to be incorporated incrementally
   • Does not reduce total complexity

2. Iterative approximation

   • Successively approximates the analysis of interest
   • Doesn't always provide the full analysis, only the essential ones

3. Subspace ensembles

   • Perform analyses on subspaces of the data
   • Predict by using an average of the individual estimators
   • Aggregate (ensemble) is often equivalent to the reference estimator

Consider a computational procedure with a complexity of \(O(n^c)\) for \(c \geq 1\) for a data matrix of size \(n \times p\).

Matloff (2012) showed the computational complexity of running the same procedure on \(n/q\) "chunks," where \(q \geq 1\) is

\[ O(n^c q^{-(c-1)}). \]

Create a lot of these "weak estimators" and aggregate them to create an estimator that is competitive with the reference estimator.

Same approach can be used for columns instead of rows.

Let \(\{P_r\}\) be a set of \(b\) be a set of orthogonal projection matrices in \(\mathcal{R}^{n \times p}\) such that:

1. Each element of each matrix in \(\{P_r\}\) is either zero or one.
2. The sum of the matrices is identity up to a permutation of the rows of \(I\). \[ \sum_{i=1}^b P_{r_i} = I\]
3. The sum of the projected \(X\) matrices is \(X\) up to a permutation of the rows of \(X\). \[ \sum_{i=1}^b P_{r_i} X = X\]

Then \(\{P_r\}\) partitions \(X\) by rows.

Ie. \(P_{r_i} X\) picks rows.

Let \(B\) be an oblique projection matrix in \(\mathcal{R}^{n \times n}\) such that:

1. Each element of each matrix in \(\{P_r\}\) is either zero or one.
2. There is at most a single value of 1 in each row of \(X\)

This can easily be generalized to \(m\) out of \(n\) bootstrapping.

Is there a valid sampling scheme where \(P\) is not a projection?

Basic approach

• Partition the rows of a data set into (approximately) evenly divided ``chunks''
• Train estimators on individual chunks
• Prediction is performed by aggregating estimators

Assumptions

• IID samples
• Reference estimator has asymptotic normal distribution with variance \(\sigma\)
• Ensemble estimates are normal and have average variance \(\sigma / c\) where \(c\) is the number of chunks.

Result: The ensemble estimator has the same asymptotic variance and mean as the ensemble estimator.

Basic approach

• Create a random partition of the rows with all subsets having roughly the same size
• Train bootstrap estimators on partitions
• Prediction is performed on aggregate estimators

Assumptions

• IID samples
• Distribution of estimator is in the set of Hadamard differentiable functions
• Subset size is asymptotically large

Result: The empirical CDF of the BLB estimator converges to the CDF of the bootstrap reference estimator. Both estimators converge to the CDF of the function at the same rate \(O_p(1/n)\).

Available from the UCI Machine Learning Repository

Element-wise comparison of medical records with personal data

5,749,132 record pairs

1. Phonetic equality of first name and family name, equality of date of birth.
2. Phonetic equality of first name, equality of day of birth.
3. Phonetic equality of first name, equality of month of birth.
4. Phonetic equality of first name, equality of year of birth.
5. Equality of complete date of birth.
6. Phonetic equality of family name, equality of sex.

require(e1017)
require(foreach)
source("eutils.r")

# Given a data set z, make a training and testing subset.
train <- sample(nrow(z), floor(nrow(z)/2))
test <- setdiff(1:nrow(z), train)

# Create the ensemble.  library(doMC) registerDoMC(cores=6)
e <- foreach(it = irpart(z[train, 1], z[train, -1], 100)) %dopar% {
    svm(x = it$x, y = as.factor(it$y))
}

# Get the individual predictions.
preds <- foreach(fit = e, .combine = cbind) %dopar% {
    as.logical(predict(fit, z[test, -1]))
}

# Aggregate the predictions
agPred <- foreach(1:nrow(preds), .combine = c) %dopar% {
    as.logical(classVoteWinner(preds[i, ]))
}

From a linear algebra viewpoint, the row-sampling techniques can just as easily be performed on columns.

However, from a statistical standpoint this is very different.

• Subsampling rows means each sample is complete
• Subsampling columns you get more samples but you don't get all of the data associated with each sample.

Sampling/partioning rows

• Chunk and average (Matloff, 2012)
• Bag of little bootstraps (Kleiner et al., 2012)

Partitioning columns

• Weighted random subspaces (Li and Zhao, 2009)

For a fixed estimator, the procedure and characteristics are defined by the subspaces the data are projected into

The combination of estimator and subspaces determine the computational savings and convergence properties.

A "column-chunking" approach to creating ensembles.

Basic approach

• Create a random partition of the columns with all subsets having roughly the same size.
• Train estimators on subspaces.
• Use cross validation to get an estimate of the error
• Create weights based on error
• Prediction is performed by the weighted aggregation of estimators

Has been found empirically to be better than equal weighting.

During Training:

1. Estimate the out-of-sample error
2. Define what you want to minimize to create a set of weights on the simplex for each estimator in the ensemble
   • Misclassification rate
   • Squared-error loss

During Estimation:

1. Create a prediction for each estimator in the ensemble
2. Weight the individual estimates
   • Weighted voting for classification
   • Multiply estimate by weight for regression

Let \(\{\widehat{\theta}\}\) be a set of \(b\) independent estimators with classification error rates \(e_1, e_2, ... e_b\). The weight assignment \[ w_i = \frac{1}{N} \log \left(\frac{1-e_i}{e_i}\right) \] with \[ N = \sum_{i=1}^b w_i \] achieves the minimal error rate.

Let \(\{\widehat{\theta}\}\) be a set of \(b\) estimators. The goal is to find \(\{w\}\) that minimize \[ \sum_{i=1}^b w_i \mathbb{P} \widehat{\theta}_i^2 + \sum_{1 \leq i < j \leq b} w_i w_j \mathbb{P} \widehat{\theta}_i \widehat{\theta}_j \] subject to \[ \sum_{i=1}^b w_i = 1 \] and \(w_i \geq 0\) \(\forall\) \(i \in \{1, ..., b\}\).

This can easily be reformulated to be a quadratic programming problem which minimizes \[ -a^{T} x + 1/2 x^T A x \] subject to \(C^T x \geq c\).

From a former Knowledge Discover and Data-mining competition

42 variables

Traffic was categorized as "good" and "bad" WRT intrusion

Task is to categorize traffic

Approach: apply random subspace ensembles

duration protocol_type service flag src_bytes dst_bytes land wrong_fragment urgent hot num_failed_logins logged_in num_compromised root_shell su_attempted num_root num_file_creations num_shells num_access_files num_outbound_cmds is_host_login is_guest_login count srv_count, serror_rate srv_serror_rate rerror_rate srv_rerror_rate same_srv_rate diff_srv_rate srv_diff_host_rate dst_host_count dst_host_srv_count dst_host_same_srv_rate dst_host_diff_srv_rate dst_host_same_src_port_rate dst_host_srv_diff_host_rate, dst_host_serror_rate dst_host_srv_serror_rate dst_host_rerror_rate dst_host_srv_rerror_rate type

Theoretical

• General approach for reducing computational complexity
• Makes use of existing results
• Work goes into proving properties of ensemble

Technological

• Turns a large class of problems from non-embarrassingly parallel to EP
• Makes use of existing implementations
• Work goes into integrating with parallel computing platform

1. Is there a general approach for assessing column sampling?
   • What about model stability?
   • Can we get general convergence results?
2. What about unsupervised learning?
3. Can we expect ensemble estimator consistency if we use both row and column subspaces?
4. Are there better bases we could be project onto?
5. How do we handle samples that are not IID?

Forthcoming

Same model abstraction as the caret package

Will use foreach for parallelism

Support is planned for data frames, matrices, bigmemory, and SciDB

Supervised and unsupervised learning will be supported

Email me for an alpha version

Email: michael.kane@yale.edu

Github: http://github.com/kaneplusplus

This talk: http://euler.stat.yale.edu/~mjk56/SubspaceEnsembles

Download the ensemble package: email me

Email me for the slidify code for this talk

Aglib (2013) The Aglib User's Guide. http://www.alglib.net/#book

Algorithm AS274 (1992) Applied Statitsics. 41(2).

F. R. Bach and M. I. Jordan (2005) Predictive low-rank decomposition for kernel methods. Proceedings of the \(22^{nd}\) International Conference on Machine Learning.

J. Baglama and L. Reichel (2005) Augmented Implicitly Restarted Lanczos Bidiagonalization Methods. SIAM Journal of Scientific Computing.

O. Chapelle (2006) Training a Support Vector Machine in the Primal. Neural Computation. 19(5) 1155-1178.

X. Li and H. Zhao (2009) Weighted random subspace methods for high dimensional data classification. Statistics and its Interface. 2 153-159.

N. Matloff (2012) Efficient R Parallel Loops on Long-Latency Platforms. Talk given at the \(43^{rd}\) Symposium on the Interface of Computing Science and Statistics.

N. Matloff (2013) Software Alchemy: Turning Complex Statistical Computations into Embarrassingly Parallel Ones. Preprint.