{'fraction_non_alphanumeric': 0.03161303949869746, 'fraction_numerical': 0.009364218827008379, 'mean_word_length': 4.634272114240381, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Combining Structured and Unstructured Randomness in Large Scale PCA

Principal Component Analysis (PCA) is a ubiquitous tool with many applications in machine learning including feature construction, subspace embedding, and outlier detection. In this paper, we present an algorithm for computing the top principal components of a dataset with a large number of rows (examples) and columns (features). Our algorithm leverages both structured and unstructured random projections to retain good accuracy while being computationally efficient. We demonstrate the technique on the winning submission the KDD 2010 Cup.


Introduction
PCA [8], and the related Singular Value Decomposition (SVD), are versatile tools in machine learning because they can compress a high dimensional dataset to a small number of dimensions. This compression can uncover hidden patterns in the data [11], reduce noise, or facilitate the application of algorithms that need not scale well to high dimensional data. In a nutshell, if X ∈ R n×p is the data matrix having n examples and p features, whose SVD is X = UΣV ⊤ , then UΣ are the principal components and V is called the loadings. Furthermore, computing the principal components of a new example x can be done via V ⊤ x. Classic algorithms for PCA cannot scale to large datasets which has led to the recent interest in randomized SVD algorithms [4,6,9]. When the data is nominally p-dimensional but distributed close to a k-dimensional subspace, these algorithms just need two streaming passes over the data plus an additional O(pk 2 ) processing to compute the top k singular values and vectors with very good accuracy. In addition they require O(pk) memory as working storage, and in a distributed context require O(pk) communication between worker nodes.
As the volumes of datasets increase, the limits of the randomized algorithms are being stretched. Datasets with millions of examples and millions of dimensions are common nowadays and by industrial standards they are considered small. As an example, when p = 10 8 and k = 300, the memory requirements of randomized algorithms are about 223 GB. However, existing randomized algorithms reveal the actual top loadings and principal components, whereas in machine learning we typically care only about the mapping R p → R k that projects a high dimensional vector x into the top principal component space. In this paper, we propose an algorithm to efficiently and accurately approximate this mapping with reduced storage requirements. We do this by first applying a structured random projection to the data vectors to a lower dimension d ≪ p and then relying on the existing randomized algorithms for truncated SVD. With high probability, our resulting projection is onto a subspace very close to the top principal component subspace. Our implementation interleaves the application of structured random projection and the steps of randomized SVD which leads to minimal overhead. In this manner, both working storage and communication requirements are reduced from O(pk) to O(dk).
Herein we focus on hashing as our way of computing structured random projections. Hashing preserves sparsity, which fits well with the sparse datasets we use in the experiments as well as other modern datasets arising from text and social graphs. Other structured random projections such as subsampled fast transforms (e.g. Hadamard, Hartley) can be used, and should be used for dense data.
A natural question is whether both the structured randomness we employ and the unstructured randomness 1 already existent in randomized SVD algorithms are necessary. Couldn't we just project upfront down to O(k) dimensions? However, projecting directly to O(k) dimensions loses a lot of accuracy. In this sense our algorithm interpolates between two extremes: a fast but crude upfront (structured) random projection followed by PCA and a slow but accurate randomized SVD algorithm that requires orthogonalization of a p × k matrix. In practice, randomized SVD is constrained by the space to store, and time to orthogonalize, a p×k matrix. We identify a parameter d, the size of upfront random projection, that should match the available hardware, while causing low distortion. For a commodity machine d = 10 5 to d = 10 6 is typical.

Relation to Prior Work
Randomized algorithms for numerical linear algebra have recently gained much attention in theoretical computer science [5]. In this work we are particularly interested in fast algorithms for truncated SVD such as those developed in Halko et. al. [4]. Using randomization to speed up SVD goes back at least a decade [7], when Papadimitriou et. al. made similar arguments to ours regarding the impact of random projections without orthogonalization on SVD. More recent advances have led to very sharp bounds [4], and practically useful algorithms [9,6,10,3] at least for datasets with either not too many rows or not too many columns. Unfortunately in machine learning we are often faced with datasets where both the number of examples and the number of features exceed the limits of the randomized SVD algorithms in today's hardware. The use of structured randomness allows us to reduce the number of features to a size that the randomized SVD algorithms can handle, without distorting the final embedding too much. The single hash utilized here is a computationally convenient technique. There exist more complicated hashing-based dimensionality reduction techniques with superior inner product preservation guarantees [1] for which analogous arguments hold.

The Algorithm
Let X ∈ R n×p be the data matrix. We assume that the features have zero empirical mean. 2 The principal components can be computed via the SVD X = UΣV ⊤ where U ∈ R n×n and V ∈ R p×p are orthogonal matrices and Σ ∈ R n×p is a (rectangular) diagonal matrix with Σ ii arranged in nonascending order. Truncating by retaining the top k singular values and corresponding vectors yields the best (in Frobenius norm [2]) rank-k approximation of X: where Σ † indicates the Moore-Penrose pseudo-inverse of Σ, and Theorem 1 will show that this mapping can be approximated by our algorithm.
Obtaining V k and Σ k could be done with any SVD algorithm, however, once n or p is large, only the randomized SVD algorithms are practical. The randomized algorithms work in two phases. In the first phase they probe the range of the input matrix with a random matrix Ω. They potentially perform multiple passes over the data, though here we will only assume one pass. Next, they orthogonalize the image of Ω under the input matrix and project onto that basis in the second pass. Even though these algorithms have been previously adapted for PCA [10,3], they assume that the orthogonalization step can be done efficiently. This is only true if either n or p is not too large, but not both.
We can easily eliminate the dependence on n by looking at the empirical covariance matrix 1 n X ⊤ X, whose top eigenvectors are V k . We can then apply a randomized SVD algorithm on this matrix. A two-pass randomized algorithm with orthogonalization of columns has space complexity O(pk) and time complexity O(pk 2 ); furthermore computing the image of Ω under the empirical covariance is Ω d×k ← Random Gaussian matrix (Ω ij ∼ N (0, 1)) 3: Y d×k ← (XH) ⊤ (XH) Ω / n ⊲ First data-parallel pass 4: Z d×k ← (XH) ⊤ (XH) Q / n ⊲ Second data-parallel pass 6: data parallel as can be seen by 1 This procedure produces V k and Σ k which can be used in a subsequent pass over the data to produce the (whitened) principal components. In practice, if the orthonormal basis for the column (feature) space fits into main memory, the algorithm is very fast, and so is suitable for large data sets where the number of features is modest, up to circa p = 10 6 on current commodity hardware.
For datasets with hundreds of millions of features, such as the adjacency matrix of an online social network, the space complexity associated with orthogonalizing the approximate basis is impractical. We therefore propose using structured randomness, without explicit materialization of the projection matrix and without orthogonalization, to reduce the number of features into the zone where randomized SVD algorithms are viable. Although not materialized, we represent the structured randomness as a matrix H ∈ R p×d . By interleaving the structured randomness with the randomized algorithm, we arrive at Algorithm 1. This mainly follows the algorithm from [4] except at step 6 where we use that the covariance is symmetric to write it as in line 3. We also compute the spectral decomposition of a k × k matrix instead of the SVD of a d × k matrix as it is sufficient to extract the loadings and the singular values.
For sparse data, e.g. text or social graphs, we find a hash based structured randomness [13] to be computationally convenient and empirically effective. Conceptually, this scheme multiplies the data by a hashing matrix H ∈ R p×d which is determined by two hash functions h : {1, . . . , p} → {1, . . . , d} and ξ : {1, . . . , p} → {±1}, with H ij | ξ,h = ξ(i)1 h(i)=j . For analytical purposes, the hash functions h and ξ are considered drawn uniformly at random from a universal family. For dense data (not presented here), subsampled Hadamard or Hartley transforms can be used.

Analysis
Our proposed algorithm is a composition of dimensionality reduction via structured randomness and established randomized SVD techniques. Because the error properties of randomized SVD techniques are well understood, we will focus on the impact of dimensionality reduction on an exact truncated SVD. Let XH =ŨΣṼ ⊤ be the SVD of the projected version of the data matrix, and let U 1 ,Σ 1 , andṼ 1 be the truncated SVD components as follows where k < min{d, n}. Although our algorithm manipulates a transformed version of the empirical covariance matrix H ⊤ X ⊤ XH, it is more convenient to consider the transformed version of the Gram matrix XHH ⊤ X ⊤ . We do this because the Gram matrix has the same nonzero eigenvalues as the covariance matrix and the corresponding eigenvectors are the whitened principal components. In particular, any mapping H which approximately preserves inner products will lead to a Gram matrix that is close to the original one in Frobenius norm. In particular, hash based structured randomness approximately preserves inner products [13].
Since the perturbation of the Gram matrix has small Frobenius norm, we can apply a classic theorem by Wedin [12], which states that the left singular subspaces associated with the original and perturbed matrices will be close (in the sense of a small canonical angle) if two things are true: first, if there is a gap between the k th largest singular value of the perturbed matrix and the (k + 1) th largest singular value of the original matrix; and second, if the difference between the Gram matrices has small norm. This reasoning leads to the following theorem.
Theorem 1 (Subspace Approximation). Let X be a data matrix with n rows, and letX = XH be the hashed data matrix. Let XX ⊤ andXX ⊤ have spectral decompositions UΣ 2 U ⊤ andŨΣ 2Ũ⊤ respectively, conformally partitioned as in (1). Let Φ be the matrix of canonical angles between the column spaces R(U 1 ) of U 1 and R(Ũ 1 ) ofŨ 1 . Suppose that there are numbers α, γ > 0 such that miñ where σ(Σ) denotes the set of diagonal values. Let r(X) denote the set of rows of X. Define If d ≥ 144 log(n/δ)/ǫ 2 and η ≤ ǫ/(18 2 log(n/δ) log(d/δ)), then with probability at least 6δ with respect to the uniform distribution over 1 H ⊤ X ⊤ will be in a space close to the whitened principal components obtained from the exact projection U 1 An important condition is a gap between the k th original singular value and the (k + 1) th perturbed singular value. As can be seen in the full proof, the Frobenius norm of the difference between the Gram matrices need only be bounded with high probability once. Then Wedin's deterministic error bound can be applied for all k. Therefore, if there is a large spectral gap at any k and we truncate the PCA at r ≥ k, then the subspaces induced by the top k whitened components will be close.

Experiment
The winning submission to the 2010 KDD Cup was a linear model developed using extensive feature engineering [14]. We focus on the winning design matrix because it is publicly available, large enough to prohibit exact decomposition, but small enough to admit standard two-pass randomized decomposition. Specifically we used the kdda dataset 3 which consists of 8.4 million training examples, 20.2 million features, and 314 million non-zero entries.
We used REDSVD [6] in sparse PCA mode as a baseline, which implements a two-pass randomized decomposition with a materialized Gaussian matrix. Computing the top 40 principal components consumes 17.5 gigabytes of RAM and takes 1591 seconds to compute the loadings and factor scores. On the same machine, algorithm 1 with d = 10 6 uses 780 megabytes of RAM and takes 1030 seconds to compute the loadings and factor scores. The decrease in memory usage is expected (due to p/d ≈ 20), while the decrease in computation time is mostly attributable to the orthogonalization step being faster due to working with smaller vectors.
To motivate the use of PCA in this setting, we augmented the design matrix by interacting the raw feature values with the first few principal components as computed by HPCA, but otherwise did not alter the baseline training procedure of the winning submission. With only 5 components the improvement in performance exceeds the difference between the winning and 2 nd place score.

Conclusion
In this paper we analyze theoretically and empirically the composition of structured and unstructured randomness with established randomized SVD techniques for the purpose of computing PCA components. Theoretically we leverage inner-product preservation guarantee to show the resulting PCA components are close to those computed by fully unstructured randomness. Empirically, the resulting algorithm is so highly scalable that experiments with some of the largest publicly available matrices were easily implemented on a commodity laptop.