--- title: "Proximity matrices when n is large" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Proximity matrices when n is large} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r, include = FALSE} knitr::opts_chunk$set(collapse = TRUE, comment = "#>") run <- requireNamespace("randomForest", quietly = TRUE) knitr::opts_chunk$set(eval = run, purl = run) ``` A proximity matrix is $n \times n$. At $n = 10^4$ that is 800 MB in double precision; at $n = 10^5$ it is 80 GB. The ensemble happily fits data at both sizes, so the matrix, not the model, is what makes the method unusable. This is not a hypothetical: it is the wall the `e2tree` work hit on the Fannie Mae and HMDA data. Three ways round it are implemented, and they give up different things. ```{r setup} library(Proximum) library(randomForest) set.seed(1) n <- 600 X <- data.frame(matrix(rnorm(n * 6), n, 6)) y <- factor(ifelse(X$X1 + X$X2 + rnorm(n) > 0, "a", "b")) training <- cbind(X, y = y) forest <- randomForest(y ~ ., data = training, ntree = 500) px <- as_proximity(forest, newdata = training) format(object.size(px), units = "auto") ``` ## Sparsity: keep every value above a threshold Thresholding and storing the result as a sparse matrix is lossless for every pair above the threshold and stores nothing for the pairs below it. ```{r sparsify} sp <- sparsify(px, threshold = 0.05) sp ``` ```{r sparsify-summary} summary(sp) ``` The saving is real. What it is not is durable: ```{r sparsify-cost} c( dense = format(object.size(px), units = "auto"), sparse = format(object.size(sp), units = "auto") ) ``` The result is a `proximity_sparse`, not a `proximity`, and no statistic in the package will take it: ```{r sparsify-refusal, error = TRUE} mantel_test(sp, px, n_perm = 99) ``` That refusal is the point rather than an omission. Every statistic here runs on the induced dissimilarity or on the doubly centred matrix, and both are dense whatever the proximity was: $1 - P$ turns every stored zero into a one, and the Gower centring leaves no zero at all. ```{r sparsify-densification} dissimilarity <- as_dissimilarity(px) c( proximity = mean(as.matrix(px) != 0), dissimilarity = mean(dissimilarity != 0), centred = mean(double_centre(dissimilarity) != 0) ) ``` So `sparsify()` is a storage format. Use it to hold a matrix between sessions or to hand it to something outside the package; `as.matrix()` spends the memory back when you want the statistics. ## Landmarks: the Nystrom approximation Compute the proximity exactly against $m \ll n$ landmark observations, then extend it to the rest by projection: $$\tilde{P} = P_{n,m}\, P_{m,m}^{-1}\, P_{m,n}.$$ Nothing of size $n \times n$ is ever formed, on the way in or on the way out. ```{r nystrom} approximation <- nystrom(forest, training, landmarks = 60, strata = training$y) approximation ``` Stratifying the landmark sample on the response is what keeps a rare class represented, and it is the argument to reach for when one class is small. ```{r nystrom-summary} summary(approximation) ``` `diagonal error` is the price. The approximation has $\tilde{P}_{ii} \ne 1$, where an exact proximity has one, and that departure is the cheapest single measure of how much the landmarks failed to span the sample. It falls as the landmarks are added: ```{r nystrom-landmarks} sapply(c(20, 60, 180), function(m) { summary(nystrom(forest, training, landmarks = m))$diagonal_error }) ``` ### What the approximation is for Unlike the sparse form, this one survives being used, because the question it answers is the geometric one. The configuration comes out of the stored factor at $O(nr^2)$ instead of $O(n^3)$, and `protest()` takes the object directly: ```{r nystrom-embedding} coordinates <- embedding(approximation, k = 2) dim(coordinates) protest(approximation, px, n_perm = 199) ``` The pairwise statistics still refuse it, for the same reason as before: there is no matrix to correlate entry by entry without building one. ```{r nystrom-refusal, error = TRUE} cka(approximation, px) ``` ## How many trees, and how much does the answer move? Neither of the above helps if the matrix is unstable, and the number of trees that settles it is a question with an answer. ```{r n-trees} required <- n_trees_required(forest, training, eps = 0.2) required attr(required, "path") ``` `autoplot()` draws the search it recorded: the criterion against the block size on log axes, with the target, the answer, and the power law fitted through the measured points. ```{r n-trees-plot, fig.width = 6, fig.height = 4} autoplot(required) ``` The answer is an integer and goes on behaving as one, so it can be handed straight back to the engine that raised the question: ```{r n-trees-integer} required + 100L ``` The criterion falls as a power of the block size, so a target set too low is not reached by any ensemble you would fit. When that happens the result is `NA` and the projection says what it would take: ```{r n-trees-unreachable} unreachable <- n_trees_required(forest, training, eps = 0.01) c(answer = unreachable, projected_trees = attr(unreachable, "projected")) ``` `stability()` asks the same question of replicates you already hold, and makes no assumption about where they came from: ```{r stability} replicates <- lapply(1:4, function(i) { as_proximity(randomForest(y ~ ., data = training, ntree = 250), newdata = training) }) agreement <- stability(replicates) agreement ``` Its plot shows every one of the $R(R-1)/2$ comparisons, with the median and the percentile interval marked. There is no band around a curve here, because there is no curve: the object holds one set of dependent agreements, and the spread of them is the whole of what it can say. ```{r stability-plot, fig.width = 6, fig.height = 3.4} autoplot(agreement) ``` ## Choosing between them | $n$ | Strategy | Cost | |---|---|---| | $\le 5{,}000$ | Dense, exact | Nothing given up | | $5{,}000$ to $50{,}000$ | Sparse for storage, dense for the statistics | Small proximities lost, memory spent again on use | | $> 50{,}000$, geometry wanted | Nystrom | Rank-$m$ approximation, no pairwise statistics | | $> 50{,}000$, a scalar wanted | Streaming | Exact, but time in place of memory | ## Streaming: never allocate it at all When the quantity of interest is a scalar, a Mantel correlation or a CKA, the matrix does not have to exist. It has to be *traversed*: every entry is needed once, and nothing needs two of them at the same time. `proximity_stream()` holds what the matrix is made of instead of the matrix, and the statistics manufacture it a block of rows at a time. ```{r stream} set.seed(2) shallow <- randomForest(y ~ ., data = training, ntree = 500, maxnodes = 8) streamed <- proximity_stream(forest, training) streamed_shallow <- proximity_stream(shallow, training) streamed ``` Read the last line before going further: at $n = 600$ with 500 trees this stream holds *more* than the matrix it stands for. The reason is in the next section, and the object reports it rather than letting you assume otherwise. The answer is the dense answer. Across the cells of `inst/simulations/streaming-cost.R` the two paths agreed to within $1.8 \times 10^{-13}$ on the statistic, on every one of the permuted statistics, and on the count of usable pairs. ```{r stream-agrees} c(streamed = cka(streamed, streamed_shallow), dense = cka(px, as_proximity(shallow, newdata = training))) ``` ### The shape this feature was promised in was wrong Earlier versions of this vignette said that `mantel_test()` would gain a `streaming` argument. It cannot have one. `mantel_test()` takes a `proximity`, and a `proximity` is a matrix that has already been built; a matrix that has already been built cannot be traversed *instead of* built. The saving is only available at the moment the object is constructed, so it belongs to the type of the input and not to a flag on the function. That is why the feature arrived as a constructor rather than as an argument. ### What it costs, which is not nothing Streaming buys memory with time, and the exchange rate is worse than it looks. The stored object is $O(nB)$ against the matrix's $O(n^2)$, but $O$ hides the constants, and below a certain size the stream is the *larger* object. Measured across fifteen cells, the leaf indicator costs 13.1 bytes per observation per tree, so the two are equal at $n \approx 1.64B$ and the stream is bigger below it: at $n = 100$ with 500 trees it holds eight times what the matrix would. It was the larger object in 8 of those 15 cells. `print()` shows both numbers, so this is checkable rather than something to reason about. Above the crossover the saving grows linearly, and at the sizes this vignette is about it is the whole game: carrying the measured constant out to $n = 100{,}000$ with 500 trees puts the indicator at some 655 MB where the matrix would need 80 GB. The time is the other half. The dense path builds the matrix once and then indexes it; the streaming path manufactures every entry each time it wants one, and a permutation test wants all of them once per permutation. On the same grid the streamed alignment took up to 6 times the dense one and the streamed Mantel test up to 12.8 times, which at $n = 800$ with 200 trees was 0.089 seconds per permutation. `n_perm` is therefore the knob that decides whether the answer arrives at all, and the block size is not: the same alignment took 2.09 seconds a row at a time and 0.069 seconds in one block, a factor of 30, for an answer that agreed to $2 \times 10^{-13}$ throughout. Blocks are for fitting in memory. The default divides a 64 MB budget by $n$. ### Two things it refuses `method = "spearman"` and the partial variant both error rather than approximate. A Spearman correlation ranks the pairs against each other, and a rank is a statement about every pair at once, so it cannot be accumulated from blocks that have been thrown away. The partial variant needs the regression fitted before the residuals can be correlated, which is two traversals and a permutation carried through both. Neither is built, and both say so.