Randomized Nonlinear Component Analysis

Introduction

Principal Component Analysis (Pearson, 1901) and Canonical Correlation Analysis (Hotelling, 1936) are two of the most popular multivariate analysis methods. They have played a crucial role in a vast array of applications since their conception a century ago.

Principal Component Analysis (PCA) rotates a collection of correlated variables into their uncorrelated principal components (also known as factors or latent variables). Principal components owe their name to the following property: the first principal component captures the maximum amount of variance in the data; successive components account for the maximum amount of remaining variance in dimensions orthogonal to the preceding ones. PCA is commonly used for dimensionality reduction, assuming that core properties of a high-dimensional sample are largely captured by a small number of principal components.

Canonical Correlation Analysis (CCA) computes linear transformations of a pair of random variables such that their projections are maximally correlated. Analogous to principal components, the projections of the pair of random variables are mutually orthogonal and ordered by their amount of explained cross-correlation. CCA is widely used to learn from multiple modalities of data (Kakade & Foster, 2007), an ability particularly useful when some of the modalities are only available at training time but keeping information about them at testing time is beneficial (Chaudhuri et al., 2009; Vapnik & Vashist, 2009).

The applications of PCA and CCA are ubiquitous. Some examples are feature extraction, time-series prediction, finance, medicine, meteorology, chemometrics, biology, neurology, natural language processing, speech recognition, computer vision or multimodal signal processing (Jolliffe, 2002).

Despite their success, an impediment of PCA and CCA for modern data analysis is that both reveal only linear relationships between the variables under study. To overcome this limitation, several nonlinear extensions have been proposed for both PCA and CCA. For PCA, these include Kernel Principal Component Analysis or KPCA (Schölkopf et al., 1999) and Autoencoder Neural Networks (Baldi & Hornik, 1989; Hinton & Salakhutdinov, 2006). For CCA, common extensions are Kernel Canonical Correlation Analysis or KCCA (Lai & Fyfe, 2000; Bach & Jordan, 2002) and Deep Canonical Correlation Analysis (Andrew et al., 2013). However, these solutions tend to have rather high computational complexity (often cubic in the sample size), are difficult to parallelize or are not accompanied by theoretical guarantees.

In a separate strand of recent research, randomized strategies have been introduced to construct features that can help reveal nonlinear patterns in data when used in conjunction with linear algorithms (Rahimi & Recht, 2008; Le et al., 2013). For basic tasks such as regression or classification, using nonlinear random features incurs little or no loss in performance compared with exact kernel methods, while achieving drastic savings in computational complexity (from cubic to linear in the sample size). Furthermore, random features are amenable to simple implementation and clean theoretical analysis.

The main contribution of this paper is to lay the foundations for nonlinear, randomized variants of PCA and CCA. Therefore, we dedicate attention to studying the spectral properties of low-rank kernel matrices constructed as sums of random feature dot-products. Our analysis relies on the recently developed matrix Bernstein inequality (Mackey et al., 2014). With little additional effort, our analysis extends to other popular multivariate analysis tools such as linear discriminant analysis, spectral clustering and the randomized dependence coefficient.

We demonstrate the effectiveness of the proposed randomized methods by experimenting with several real-world data and comparing against the state-of-the-art Deep Canonical Correlation Analysis (Andrew et al., 2013). As a novel application of the presented methods, we derive an algorithm to learn using privileged information (Vapnik & Vashist, 2009) and a scalable strategy to train nonlinear autoencoder neural networks. Additional numerical simulations are provided to validate the tightness of the concentration bounds derived in our theoretical analysis. Lastly, the presented methods are very simple to implement; we provide R source code at:

There has been a recent stream of research in kernel approximations via randomized feature maps since the seminal work of Rahimi & Recht (2008). For instance, their extensions to dot-product kernels (Kar & Karnick, 2012) and polynomial kernels (Hamid et al., 2014); the development of advanced sampling techniques using Quasi-Monte-Carlo methods (Yang et al., 2014) or their accelerated computation via fast Walsh-Hadamard transforms (Le et al., 2013).

The use of randomized techniques for kernelized component analysis methods dates back to (Achlioptas et al., 2002), where three kernel sub-sampling strategies were suggested to speed up KPCA. On the other hand, (Avron et al., 2013) made use of randomized Walsh-Hadamard transforms to adapt linear CCA to large-scale datasets. The use of non-linear random features is more scarce and has only appeared twice in previous literature. First, Lopez-Paz et al. (2013) defined the dependence statistic RDC as the largest canonical correlation between two sets of copula random projections. Second, McWilliams et al. (2013) used the Nyström method to define a randomized feature map and performed CCA to achieve state-of-the-art semi-supervised learning.

Random Nonlinear Features

We start our presentation by recalling a few key aspects of nonlinear random features.

Consider the class $\mathcal{F}_{p}$ of functions whose weights decay faster than some sampling distribution $p$; formally:

Kernel machines, Gaussian processes, AdaBoost, and neural networks are models that fit within this function class.

for a suitable loss function $c(\hat{y},y)$ that penalizes departure of $f_{m}(\bm{x})$ from the true label $y$; for us, the least-squares loss $c(\hat{y},y)=(\hat{y}-y)^{2}$ will be most convenient.

Using the (precomputed) nonlinear random features $\bm{z}(\bm{X})$ ultimately transforms the nonconvex optimization of (3) into a least-squares problem of the form

This form remarkably simplifies computation (in practice, we solve a regularized version of it), while incurring only a bounded error. Theorem 1 formalizes this claim.

Solving (5) takes $O(ndm+m^{2}n)$ operations, while testing $t$ points on the fitted model takes $O(tdm)$ operations. Recent techniques that use subsampled Hadamard randomized transforms (Le et al., 2013) allow faster computation of the random features, yielding $O(n\log(d)m+m^{2}n)$ operations to solve (5) and $O(t\log(d)m)$ to test $t$ new points.

It is of special interest that randomized algorithms are in many cases more robust than their deterministic analogues (Mahoney, 2011) because of the implicit regularization induced by randomness.

where $p(\bm{w})$ is set to be the inverse Fourier transform of $k$ and $b_{i}\sim\mathcal{U}(0,2\pi)$ (Rahimi & Recht, 2008)—e.g., the Gaussian kernel $k(\bm{x},\bm{y})=\exp(-s\|\bm{x}-\bm{y}\|_{2}^{2})$ can be approximated using $\bm{w}_{i}\sim\mathcal{N}(\bm{0},2s\bm{I})$.

The focus of this paper is on building scalable kernel component analysis methods which not only exploit these approximations but are also accompanied by theoretical guarantees.

Importantly, our analysis extends straight-forwardly to features constructed using the Nyström method (Williams & Seeger, 2001) when its basis are bounded and sampled at random. Recent theoretical and empirical evidence suggest the superiority of the Nyström method when compared to the aforementioned Fourier features (Yang et al., 2012). However, we did not experience large differences (Section 5), while random Fourier features are faster to compute and do not need to be stored at test time (Le et al., 2013).

Principal Component Analysis (PCA)

For a centered data matrix (zero mean columns) $\bm{X}$, PCA requires computing the (full) singular value decomposition

where $\bm{\Sigma}$ is a diagonal matrix containing the singular values of $\bm{X}$ in decreasing order. The principal components are computed via the linear transformation $\bm{X}\bm{F}$.

Nonlinear variants of PCA are also known; notably,

Kernel PCA (Schölkopf et al., 1999) uses the kernel trick to embed data into a high-dimension Reproducing Kernel Hilbert Space, where regular PCA is performed. Computation of the principal components reduces to an eigenvalue problem, which takes $O(n^{3})$ operations.

Autoencoders (Hinton & Salakhutdinov, 2006) are artificial neural networks configured to learn their own input. They are trained to learn compressed representations of data. The transformation computed by a linear autoencoder with a bottleneck of size $r<d$ is the projection into the subspace spanned by the first $r$ principal components of the training data (Baldi & Hornik, 1989).

We propose RPCA, a nonlinear randomized variant of PCA. We may view RPCA as a low-rank approximation of KPCA when the latter is equipped with a shift-invariant kernel. RPCA may be thus understood as linear PCA on a randomized nonlinear mapping of the data. Schematically,

The computational complexity is $O(n^{3})$ for KPCA, $O(d^{2}n)$ for PCA and $O(m^{2}n)$ for RPCA. PCA and RPCA are both linear in the sample size $n$. When using nonlinear features as in (2), PCA loadings are no longer linear transformations but approximations of nonlinear functions belonging to the function class $\mathcal{F}_{p}$ described in Section 2.

To analyze this convergence we appeal to the recently proposed Matrix Bernstein Inequality. In the theorem below and henceforth $\|\bm{X}\|$ denotes the operator norm.

The convergence rate of RPCA to its exact kernel counterpart KPCA is expressed by the following theorem, which actually invokes the Hermitian matrix version of Theorem 3, and hence depends on $d$ instead of $2d$, and uses matrix squares when defining the variance $\sigma^{2}$.

We follow a derivation similar to Tropp (2012).

Next, taking all summands $\bm{E}_{i}$ together we obtain

Observe that random features and kernel evaluations are upper-bounded by $1$; thus both $B$ and $\|\bm{K}\|$ are upper-bounded by $n$, yielding the bound (7). ∎

To obtain a characterization in relative-error, we can divide both sides of (7) by $\|\bm{K}\|$. This results in a bound that depends on $n$ logarithmically (since $\|\bm{K}\|=O(n)$). Moreover, additional information may be extracted from the tail-probability version of Theorem . Please refer to Section 5.1 for additional discussion on this aspect.

Before closing this section, we mention a byproduct of our above analysis.

Extension to Spectral Clustering. Spectral clustering (Luxburg, 2007) uses the spectrum of $\bm{K}$ to perform dimensionality reduction before applying $k$-means. Therefore, the analysis of RPCA may be easily extended to obtain a randomized and nonlinear variant of spectral clustering.

Canonical Correlation Analysis (CCA)

where $\bm{C}_{XY}$ is the covariance $\text{cov}(\bm{X},\bm{Y})$, while the diagonal terms $\gamma\bm{I}$ act as regularization.

In another words, CCA processes two different views of the same data (i.e., speech audio signals and paired speaker video frames) and returns their maximally correlated linear transformations. This is particularly useful when the two views are available at training time, but only one of them is available at test time (Kakade & Foster, 2007; Chaudhuri et al., 2009; Vapnik & Vashist, 2009).

Several nonlinear extensions of CCA have been proposed:

Kernel Canonical Correlation Analysis or KCCA (Lai & Fyfe, 2000; Bach & Jordan, 2002) uses the kernel trick to derive a nonparametric, nonlinear regularized CCA algorithm. Its exact computation takes time $O(n^{3})$.

Deep Canonical Correlation Analysis or DCCA (Andrew et al., 2013) feeds the pair of input variables through a deep neural network. Transformation weights are learnt using gradient descent to maximize the correlation of the output mappings.

The computational complexity is $O(n^{3})$ for KCCA, $O((p^{2}+q^{2})n)$ for CCA and $O((m_{x}^{2}+m_{y}^{2})n)$ for RCCA. CCA and RCCA are both linear in the sample size $n$.

As with PCA, we are interested in characterizing the convergence rate of RCCA to its exact kernel counterpart KCCA as $m_{x}$ and $m_{y}$ grow. The solution of KCCA is the eigensystem of the matrix $\bm{R}^{-1}\bm{L}$, where,

and $(\gamma_{x},\gamma_{y})$ are positive regularizers mandatory to avoid spurious $\pm 1$ correlations (Bach & Jordan, 2002). Theorem 4 characterizes the convergence rate of RCCA to KCCA. Let $\hat{\bm{R}}^{-1}$ and $\hat{\bm{L}}$ be the approximations to (8) and (9) obtained by using $m$ random features; that is

As the matrices are block-diagonal, we have

We analyze the first term of the maximum; the latter can be analyzed analogously. Let $\hat{\bm{A}}:=(\hat{\bm{K}}_{x}+\gamma_{x}\bm{I})^{-1}$ and $\bm{A}:=(\bm{K}_{x}+\gamma_{x}\bm{I})^{-1}$. Define the individual error terms

Recall that the $m_{x}+m_{y}$ random features are sampled i.i.d. and that the data matrices $\bm{X}$, $\bm{Y}$ are constant. Therefore, the random matrices $\bm{\hat{}}{\bm{K}}_{x}^{(1)},\ldots,\bm{\hat{}}{\bm{K}}_{x}^{(m_{x})},\bm{\hat{}}{\bm{K}}_{y}^{(1)},\ldots,\bm{\hat{}}{\bm{K}}_{y}^{(m_{y})}$ are i.i.d. random variables. Hence, their expectations factorize:

and the norm of this deviation is bounded as

The inequality follows by applying Hölder twice after using the triangle inequality. We now turn to the issue of computing the variance, which is defined as

Consider first second argument of the maximum above, for which we expand an individual term in the summand:

An invocation of Jensen on the definition of $\sigma^{2}$ along with the two bounds above yields the worst-case estimate

We may now appeal to the matrix Bernstein inequality (Theorem 3) to obtain the bound

Before concluding this section, we briefly comment on two easy extensions of our above result.

Extension to RDC.

The Randomized Dependence Coefficient or RDC (Lopez-Paz et al., 2013) is defined as the largest canonical correlation of RCCA when performed on the copula transformation of the data matrices of $\bm{X}$ and $\bm{Y}$. Our analysis applies to the further understanding of RDC.

Experiments

We investigate the performance of RCCA in multiple experiments with real-world data against state-of-the-art algorithms. Section 5.3 provides a novel algorithm based on RCCA to perform learning using privileged information (Vapnik & Vashist, 2009). Section 5.4 introduces the use of RPCA as a tool to train autoencoders in a scalable manner.

We set our random (Fourier) features to approximate the Gaussian kernel, as described in the second paragraph of Section 2.1. We also compare to the Nyström method, set to construct an $m-$dimensional feature space formed by the evaluations of the Gaussian kernel on $m$ random points from the training set (Yang et al., 2012). Gaussian kernel widths $\{s_{x},s_{y}\}$ are set using the median heuristic.

Figure 1 depicts the value of the norms from equations (7, 12) as the parameters $\{n,m,\gamma\}$ vary, when averaged over a total of $1000$ random samples $\{\bm{X}_{i},\bm{Y}_{i}\}_{i=1}^{1000}$. The simulations agree with the presented theoretical analysis: the sample size $n$ and regularization parameter $\gamma$ exhibit a linear effect, while increasing the number of random features $m$ induces an $O(m^{-1/2})$ reduction in error (the closest function $O(m^{-1/2})$ is overlaid in red for comparison).

2 Canonical Correlation Analysis

We compare three variants of CCA on the task of learning correlated features from two modalities of the same data: linear CCA, state-of-the-art Deep CCA (Andrew et al., 2013) and the proposed (Fourier and Nyström based) RCCA. We were unable to run exact KCCA on the proposed datasets due to its cubic complexity; other low-rank approximations such as the one of Arora & Livescu (2012) were shown inferior to DCCA, and hence omitted in our analysis.

We replicate the two experiments presented in Andrew et al. (2013). The task is to measure performance as the accumulated correlation between the canonical variables associated with the largest training canonical correlations on some unseen test data. The participating datasets are MNIST and XRMB, which are introduced in the following.

Learn correlated representations between the left and right halves of the MNIST images (LeCun & Cortes, 1998). Each image has a width and height of 28 pixels; therefore, each of the two views of CCA consists on 392 features. 54000 random samples are used for training, 10000 for testing and 6000 to cross-validate the parameters of (D)CCA.

X-Ray Microbeam Speech Data.

Learn correlated representations of simultaneous acoustic and articulatory speech measurements (Westbury, 1994). The articulatory measurements describe the position of the speaker’s lips, tongue and jaws for seven consecutive frames, yielding a 112-dimensional vector at each point in time; the acoustic measurements are the MFCCs for the same frames, producing a 273-dimensional vector for each point in time. 30000 random samples are used for training, 10000 for testing and 10000 to cross-validate the parameters of (D)CCA.

Summary of Results.

Table 1 shows the sum of the largest canonical correlations (corr.) obtained by each CCA variant and their running times (minutes, single 1.8GHz core) on the MNIST and XRMB test sets. Given enough random projections ($m=m_{x}=m_{y}$), RCCA is able to explain the most amount of test correlation while running drastically faster than DCCARunning times for DCCA correspond to a single cross-validation iteration of its ten hyper-parameters. DCCA has 2 layers for MNIST and 8 layers for XRMB.. Moreover, when using random features (i) the number of weights required to be stored at test time for RCCA is up to two orders of magnitude lower than for DCCA and (ii) the use of Fastfood multiplications (Le et al., 2013) allows much faster model evaluation.

Parameter Selection.

No parameters were tuned for RCCA: the kernel widths were heuristically set and CCA regularization is implicitly provided by the use of randomness (thus set to $10^{-8}$). The number of random features $m$ can be set to the maximum value that fits within the available (training or test time) computational budget. On the contrary, previous state-of-the-art DCCA has ten parameters (two autoencoder parameters for pretraining, number of hidden layers, number of hidden units and CCA regularizers for each view), which were cross-validated using the grids described in Andrew et al. (2013). Cross-validating RCCA parameters did not significantly improve performance.

If desired, further speed improvements for RCCA could be achieved by distributing the computation of covariance matrices over several CPUs or GPUs, and by making use of truncated SVD routines (Baglama & Reichel, 2006).

3 Learning Using Privileged Information

In Vapnik’s Learning Using Privileged Information (LUPI) paradigm (Vapnik & Vashist, 2009) the learner has access to a set of privileged features or information $\bm{X}_{\star}$, exclusive of training time. These features are understood as helpful high-level “teacher explanations” about each of the training samples. The challenge is to build algorithms able to extract information from this privileged features at training time in order to build a better classifier at test time. We propose to use RCCA to construct a highly correlated subspace between the regular features $\bm{X}$ and the privileged features $\bm{X}_{\star}$, accessible at test time through a nonlinear transformation of $\bm{X}$.

We perform 14 random training/test partitions of $1000$ samples each. Each partition groups a random subset of $10$ animals as class “” and a second random subset of $10$ animals as class “$1$”. Hence, each experiment is a different, challenging binary classification problem. Figure 2 shows the test classification accuracy of a linear SVM when using as features the images’ SURF descriptors or the RCCA “semi-privileged” features. As a side note, directly using the high-level attributes yields $100\%$ accuracy. The cost parameter of the linear SVM is cross-validated on the grid $[10^{-4},\ldots,10^{4}]$. We observe an average improvement of $14\%$ in classification when using the RCCA basis instead of the image features alone. Results are statistically significant respect to a paired Wilcoxon test on a $95\%$ confidence interval. The SVM+ algorithm (Vapnik & Vashist, 2009) did not improve on the regular SVM using SURF descriptors.

4 Randomized Autoencoders

Figure 3 shows the reconstruction of unseen MNIST and CIFAR-10 images after being compressed with RPCA. The number random projections was set to $m=2000$. The number of latent dimensions was set to $d=20$ for MNIST, and $d=40$ (first row) or $d=100$ (second row) for CIFAR-10. Training took under 200 seconds for each full dataset.

We thank the anonymous reviewers for their numerous comments, and the fruitful discussions had with Yarin Gal, Mark van der Wilk and Maxim Rabinovich. Lopez-Paz is supported by Obra Social “la Caixa”.

References