Nonparametric Canonical Correlation Analysis

Introduction

A common task in data analysis is to reveal the common variability in multiple views of the same phenomenon, while suppressing view-specific noise factors. Canonical correlation analysis (CCA) (Hotelling, 1936) is a classical statistical technique that targets this goal. In CCA, linear projections of two random vectors are sought, such that the resulting low-dimensional vectors are maximally correlated. This tool has found widespread use in various fields, including recent application to natural language processing (Dhillon et al., 2011), speech recognition (Arora and Livescu, 2013), genomics (Witten and Tibshirani, 2009), and cross-modal retrieval (Gong et al., 2014).

One of the shortcomings of CCA is its restriction to linear mappings, since many real-world multi-view datasets exhibit highly nonlinear relationships. To overcome this limitation, several extensions of CCA have been proposed for finding maximally correlated nonlinear projections. In kernel CCA (KCCA) (Akaho, 2001; Melzer et al., 2001; Bach and Jordan, 2002; Hardoon et al., 2004), these nonlinear mappings are chosen from two reproducing kernel Hilbert spaces (RKHS). In deep CCA (DCCA) (Andrew et al., 2013), the projections are obtained from two deep neural networks that are trained to output maximally correlated vectors. Nonparametric CCA-type methods, which are not limited to specific function classes, include the alternating conditional expectations (ACE) algorithm and its extensions (Breiman and Friedman, 1985; Balakrishnan et al., 2012; Makur et al., 2015). Nonlinear CCA methods are advantageous over linear CCA in a range of applications (Hardoon et al., 2004; Melzer et al., 2001; Wang et al., 2015b). However, existing nonlinear CCA approaches are very computationally demanding, and are often impractical to apply on large data.

Interestingly, the problem of finding the most correlated nonlinear projections of two random variables has been studied by Lancaster (1958) and Hannan (1961), long before the derivation of KCCA, DCCA and ACE. They characterized the optimal projections in the population setting, without restricting the solution to an RKHS or to have any particular parametric form. However, these theoretical results have not inspired practical algorithms.

In this paper, we revisit Lancaster’s theory, and use it to devise a practical algorithm for nonparametric CCA (NCCA). Specifically, we show that the solution to the nonlinear CCA problem can be expressed in terms of the singular value decomposition (SVD) of a certain operator, which is defined via the population density. Therefore, to obtain a practical method, we estimate the density from training data and use the estimate in the solution. The resulting algorithm reduces to solving an eigenvalue system with a particular kernel that depends on the joint distribution between the views. While superficially similar to other eigenvalue methods, it is fundamentally different from them and in particular has crucial advantages over KCCA. For example, unlike KCCA, NCCA does not involve computing the inverse of any matrix, making it computationally feasible on large data where KCCA (even using approximation techniques) is impractical. We elucidate this and other contrasts in Sec. 3 below. We show that NCCA achieves state-of-the art performance, while being much more computationally efficient than KCCA and DCCA.

In certain situations, nonlinearity is needed for one view but not for the other. In such cases, it may be advantageous to constrain the projection of the second view to be linear. An additional contribution of this paper is the derivation of a closed-form solution to this partially linear CCA (PLCCA) problem in the population setting. We show that PLCCA has essentially the same form as linear CCA, but with the optimal linear predictor term in CCA replaced by an optimal nonlinear predictor in PLCCA. Thus, moving from the population setting to sample data entails simply using nonlinear regression to estimate this predictor. The resulting algorithm is efficient and, as we demonstrate on realistic data, sometimes matches DCCA and significantly outperforms CCA and KCCA.

Background

To facilitate the analogy with partially linear CCA (Sec. 3.2), we note that the CCA solution can also be expressed in terms of the optimal predictor (in the mean squared error sense) of $X$ from $Y$, given by $\hat{X}=\boldsymbol{\Sigma}_{xy}\boldsymbol{\Sigma}_{yy}^{-1}Y$, and its covariance $\boldsymbol{\Sigma}_{\hat{x}\hat{x}}=\boldsymbol{\Sigma}_{xy}\boldsymbol{\Sigma}_{yy}^{-1}\boldsymbol{\Sigma}_{yx}$. Specifically, $\mathbf{U}$ corresponds to the eigenvectors of $\mathbf{K}=\mathbf{T}\mathbf{T}^{\top}=\boldsymbol{\Sigma}_{xx}^{-1/2}\boldsymbol{\Sigma}_{\hat{x}\hat{x}}\boldsymbol{\Sigma}_{xx}^{-1/2}$, and, by algebraic manipulation, the optimal projections can be written as

where $\mathbf{D}$ is a diagonal matrix with the top $L$ eigenvalues of $\mathbf{K}$ on its diagonal.

In the more recently proposed DCCA approach (Andrew et al., 2013), $\boldsymbol{\mathcal{A}}$ and $\boldsymbol{\mathcal{B}}$ are the families of functions that can be implemented using two deep neural networks of predefined architecture. As a parametric method, DCCA scales better than approximate KCCA for large datasets (Wang et al., 2015b).

Lancaster (1958) studied a variant of problem (3), where $\boldsymbol{\mathcal{A}}$ and $\boldsymbol{\mathcal{B}}$ are the families of all measurable functions. This setting may seem too unrestrictive. However, it turns out that in the population setting, the optimal projections are well-defined even without imposing smoothness in any way. Lancaster characterized the optimal (possibly nonlinear) mappings $f_{i}$ and $g_{i}$ for one-dimensional $X$ and $Y$ ($D_{x}=D_{y}=1$). In particular, he showed that if $X,Y$ are jointly Gaussian, then the optimal projections are Hermite polynomials. Eagleson (1964) extended this analysis to the Gamma, Poisson, binomial, negative binomial, and hypergeometric distributions. Hannan (1961) gave Lancaster’s characterization a functional analysis interpretation, which confirmed its validity also for multi-dimensional views.

Lancaster’s population solution has never been used for devising a practical CCA algorithm that works with sample data. Here, we revisit Lancaster’s result, extend it to a semi-parametric setting, and devise practical algorithms that work with sample data. Clearly, in the finite-sample setting, it is necessary to impose smoothness. Our approach to imposing smoothness is different from KCCA, which formulates the problem as one of finding the optimal smooth solution (in an RKHS) and then approximates it from samples. Here, we first derive the optimal solution among all (not necessarily smooth) measurable functions, and then approximate it by using smoothed versions of the true densities, which we estimate from data. As we show, the resulting algorithm has significant advantages over KCCA.

Nonparametric and partially linear CCA

We treat the following two variants of the nonlinear CCA problem (3): (i) Nonparametric CCA in which both $\boldsymbol{\mathcal{A}}$ and $\boldsymbol{\mathcal{B}}$ are the sets of all (nonparametric) measurable functions; (ii) Partially linear CCA (PLCCA), in which $\boldsymbol{\mathcal{A}}$ is the set of all linear functions $\mathbf{f}(\mathbf{x})=\mathbf{W}^{T}\mathbf{x}$, and $\boldsymbol{\mathcal{B}}$ is the set of all (nonparametric) measurable functions $\mathbf{g}(\mathbf{y})$. We start by deriving closed-form solutions in the population setting, and then plug in an empirical estimate of $p(\mathbf{x},\mathbf{y})$.

whereFormally, $s(\mathbf{x},\mathbf{y})$ is the Radon-Nikodym derivative of the joint probability measure w.r.t. the product of marginal measures, assuming the former is absolutely continuous w.r.t. the latter.

where $\delta_{ij}$ is Kronecker’s delta function.

When $\mathcal{S}$ is a compact operator, the solution to problem (6) can be expressed in terms of its SVD (see e.g., (Bolla, 2013, Proposition A.2.8)). Specifically, in this case $\mathcal{S}$ possesses a discrete set of singular values $\sigma_{1}\geq\sigma_{2}\geq\dots$ and corresponding left and right singular functions $\psi_{i}\in\mathcal{H}_{x},\phi_{i}\in\mathcal{H}_{y}$, and the maximal value of the objective in (6) is precisely $\sigma_{1}+\ldots+\sigma_{L}$ and is attained with

This shows that NCCA resembles other spectral dimensionality reduction algorithms, in that the projections are the eigenfunctions of some kernel. However, in NCCA, the kernel is not specified by the user. From (8), we see that $k(\mathbf{x},\mathbf{x}^{\prime})$ corresponds to the inner product between $s(\mathbf{x},\cdot)$ and $s(\mathbf{x}^{\prime},\cdot)$ (equivalently $p(\mathbf{y}|\mathbf{x})/p(\mathbf{y})$ and $p(\mathbf{y}|\mathbf{x}^{\prime})/p(\mathbf{y})$). Therefore, as visualized in Fig. 1, in NCCA $\mathbf{x}$ is considered similar to $\mathbf{x}^{\prime}$ if the conditional distribution of $Y$ given $X=\mathbf{x}$ is similar to that of $Y$ given $X=\mathbf{x}^{\prime}$.

A sufficient condition for $\mathcal{S}$ to be compact is that it be a Hilbert-Schmidt operator, i.e., that

2 Partially linear CCA (PLCCA)

The above derivation of NCCA can be easily adapted to cases in which $\boldsymbol{\mathcal{A}}$ and $\boldsymbol{\mathcal{B}}$ are different families of functions. As an example, we next derive PLCCA, in which $\boldsymbol{\mathcal{A}}$ is the set of all linear functions of $X$ while $\boldsymbol{\mathcal{B}}$ is still the set of all (nonparametric) measurable functions of $Y$.

Having determined the optimal $\mathbf{f}(\mathbf{x})=\mathbf{W}^{\top}\mathbf{x}$, we can compute the optimal $\mathbf{g}(\mathbf{y})$ using the following lemmaA simpler version of this lemma, in which $\mathbf{f}(\mathbf{x})=\mathbf{y}$ and $\mathbf{g}$ is linear, appeared in Eldar and Oppenheim (2003). The proof of Lemma 3.1 is provided in the Supplemntary Material and follows closely that of (Eldar and Oppenheim, 2003, Theorem 1)..

Substituting $\mathbf{f}(\mathbf{x})=\mathbf{W}^{\top}\mathbf{x}=\mathbf{U}^{\top}\boldsymbol{\Sigma}_{xx}^{-1/2}\mathbf{x}$ into (12), we obtain that the partially linear CCA projections are

where $\mathbf{D}$ is the diagonal $L\times L$ matrix that has the top $L$ eigenvalues of $\mathbf{K}_{\text{PL}}$ on its diagonal.

Comparing (13) with (2), we see that PLCCA has the exact same form as CCA. The only difference is that here $\hat{X}$ is the optimal nonlinear predictor of $X$ from $Y$ (a nonlinear function of $Y$), whereas in CCA, $\hat{X}$ corresponded to the best linear predictor of $X$ from $Y$ (a linear function of $Y$).

3 Practical implementations

The NCCA and PLCCA solutions require knowing the joint probability density $p(\mathbf{x},\mathbf{y})$ of the views. Given a set of training data $\{(\mathbf{x}_{i},\mathbf{y}_{i})\}_{i=1}^{N}$ drawn independently from $p(\mathbf{x},\mathbf{y})$, we can estimate $p(\mathbf{x},\mathbf{y})$ and plug it into our formulas. There are many ways of estimating this density. We next present the algorithms resulting from using one particular choice, namely the kernel density estimates (KDEs)

where $w(t)\propto e^{-t/2}$ is the Gaussian kernel, and $\sigma_{x}$ and $\sigma_{y}$ are the kernel widths of the two views.

We note that, theoretically, KDEs suffer from the curse of dimensionality, and use of other density estimation methods is certainly possible. However, we make two important observations. First, real-world data sets often have low-dimensional manifold structure, and the KDE accuracy is affected only by the intrinsic dimensionality. As shown in (Ozakin and Gray, 2009), if the data lies on an $r$-dimensional manifold, then the KDE converges to the true density at a rate ofThis requires normalizing the KDE differently, but the scaling cancels out in $s(\mathbf{x},\mathbf{y})=p(\mathbf{x},\mathbf{y})/p(\mathbf{x})p(\mathbf{y})$. $\mathcal{O}(n^{-\frac{4}{r+4}})$. Indeed, KDEs have been shown to work well in practice in relatively high dimensions (Georgescu et al., 2003), as is also confirmed in our experiments. Second, the NCCA algorithm resulting from working with KDEs involves the same Gaussian affinity matrices used in (Gaussian kernel) KCCA. Thus, intuitively, the amount of smoothness required for obtaining accurate results in high dimensions is similar for NCCA and KCCA. Nevertheless, NCCA has a clear advantage over KCCA in terms of both performance and computation.

The NCCA implementation, with the specific choice of Gaussian KDEs, is given in Algorithm 1. If the input dimensionality is too high, we first perform PCA on the inputs for more robust density estimates. To make our algorithm computationally efficient, we truncate the Gaussian affinities $\mathbf{W}^{x}_{ij}$ to zero if $\mathbf{x}_{i}$ is not within the $k$-nearest neighbors of $\mathbf{x}_{j}$ (similarly for view 2). This leads to a sparse matrix $\SS$, whose SVD can be computed efficiently.

To obtain out-of-sample mapping for a new view 1 test sample $\mathbf{x}$, we use the Nyström method (Williams and Seeger, 2001), which avoids recomputing SVD. Specifically, recall that the view 1 projections are the eigenfunctions of the positive definite kernel $k(\mathbf{x},\mathbf{x}^{\prime})$ of (8). Computing this kernel function between $\mathbf{x}$ and the training samples leads to (notice the corresponding view 2 input of $\mathbf{x}$ is not needed)

Thus, applying the Nyström method, the projections of $\mathbf{x}$ can be approximated as

for $i=1,\dots,L+1$, where $\sigma_{i}$ is the $i$th singular value of $\SS$. The second equality follows from substituting (16) and using the fact that $\mathbf{f}_{i}$ and $\mathbf{g}_{i}$ are singular vectors of $\SS$. Note again that since the affinity matrices are sparse, the mappings are computed via fast sparse matrix multiplication.

Notice that NCCA is not equivalent to KCCA with any kernel. KCCA requires two kernels, each of which only sees one view; the NCCA kernel (8) depends on both views through their joint distribution. In terms of practical implementation, our KDE-based NCCA solves a different eigenproblem and does not involve any full matrix inverses. Indeed, both methods compute the SVD of the matrix $\mathbf{Q}_{x}^{-1}\mathbf{W}^{x}\mathbf{W}^{y}\mathbf{Q}_{y}^{-1}$. However, in NCCA, $\mathbf{Q}_{x},\mathbf{Q}_{y}$ are diagonal matrices containing the sums of rows/columns of $\mathbf{W}^{x}/\mathbf{W}^{y}$, whereas in KCCA, $\mathbf{Q}_{x}=\mathbf{W}^{x}+r_{x}\mathbf{I}$, $\mathbf{Q}_{y}=\mathbf{W}^{y}+r_{y}\mathbf{I}$, for some positive regularization parameters $r_{x},r_{y}$. Moreover, in NCCA this factorization gives the projections, whereas in KCCA it gives the coefficients in the RKHS.

An additional key distinction is that NCCA does not require regularization in order to be well defined. In contrast, KCCA must use regularization, as otherwise the matrix it factorizes collapses to the identity matrix, and the resulting projections are meaningless. This is due to the fact that KCCA attempts to estimate covariances in the infinite-dimensional feature space, whereas NCCA is based on estimating probability densities in the primal space.

The resulting computational differences are striking. The number of training samples $N$ is often such that the $N\times N$ matrices in either NCCA or KCCA cannot even be stored in memory. However, these matrices are sparse, with only $kN$ entries if we retain $k$ neighbors. Therefore, in NCCA the storage problem is alleviated and matrix multiplication and eigendecomposition are $O(kN^{2})$ operations instead of $O(N^{3})$. In KCCA, one cannot take advantage of truncated kernel affinities, because of the need to compute the inverses of kernel matrices, which are in general not sparse, so direct computation is often infeasible in terms of both memory and time. Low-rank KCCA approximations (as used in our experiments below) with rank $M$ have a time complexity $O(M^{3}+M^{2}N)$, which is still challenging with typical ranks in the thousands or tens of thousands.

Related work

Several recent multi-view learning algorithms use products or sums of single-view affinity matrices, diffusion matrices, or Markov transition matrices. The combined kernels constructed in these methods resemble our matrix $\SS=\mathbf{W}^{x}\mathbf{W}^{y}$. Such an approach has been used, for example, for multi-view spectral clustering (de Sa, 2005; Zhou and Burges, 2007; Kumar et al., 2011), metric fusion (Wang et al., 2012), common manifold learning (Lederman and Talmon, 2014), and multi-view nonlinear system identification (Boots and Gordon, 2012). Note, however, that in NCCA the matrix $\SS$ corresponds to the product $\mathbf{W}^{x}\mathbf{W}^{y}$ only when using a separable Gaussian kernel for estimating the joint density $p(\mathbf{x},\mathbf{y})$. If a non-separable density estimate is used, then the matrix $\SS$ no longer resembles the previously proposed multi-view kernels. Furthermore, although algorithmically similar, NCCA arises from a completely different motivation: It maximizes the correlation between the views, whereas these other methods do not.

Experiments

In the following experiments, we compare PLCCA/NCCA with linear CCA, two kernel CCA approximations using random Fourier features (FKCCA, Lopez-Paz et al. (2014)) and Nyström approximation (NKCCA, Williams and Seeger (2001)) as described in Wang et al. (2015b), and deep CCA (DCCA, Andrew et al. (2013)).

We begin with the 2D synthetic dataset ($1000$ training samples) in Fig. 2(a,b), where samples of the two input manifolds are colored according to their common degree of freedom. Clearly, a linear mapping in view 1 cannot unfold the manifold to align the two views, and linear CCA indeed fails (results not shown). We extract a one-dimensional projection for each view using different nonlinear CCAs, and plot the projection $g(\mathbf{y})$ vs. $f(\mathbf{x})$ of test data (a different set of $1000$ random samples from the same distribution) in Fig. 2(c-f). Since the second view is essentially a linear manifold (plus noise), for NKCCA we use a linear kernel in view 2 and a Gaussian kernel in view 1, and for DCCA we use a linear network for view 2 and two hidden layers of $512$ ReLU units for view 1. Overall, NCCA achieves better alignment of the views while compressing the noise (variations not described by the common degree of freedom). While DCCA also succeeds in unfolding the view 1 manifold, it fails to compress the noise.

The University of Wisconsin X-Ray Micro-Beam (XRMB) corpus (Westbury, 1994) consists of simultaneously recorded speech and articulatory measurements. Following Andrew et al. (2013) and Lopez-Paz et al. (2014), the acoustic view inputs are 39D Mel-frequencey cepstral coefficients and the articulatory view inputs are horizontal/vertical displacement of 8 pellets attached to different parts of the vocal tract, each then concatenated over a 7-frame context window, for speaker ’JW11’. As in (Lopez-Paz et al., 2014), we randomly shuffle the frames and generate splits of $30K$/$10K$/$11K$ frames for training/tuning/testing, and we refer to the result as the ’JW11-s’ setup (random splits better satisfy the i.i.d. assumption of train/tune/test data than splits by utterances as in (Andrew et al., 2013)). We extract $112D$ projections with each algorithm and measure the total correlation between the two views of the test set, after an additional $112D$ linear CCA. As in prior work, for both FKCCA and NKCCA we use rank-$6000$ approximations for the kernel matrices; for DCCA we use two ReLU (Nair and Hinton, 2010) hidden layers of width $1800$/$1200$ for view 1/2 respectively and run stochastic optimization with minibatch size $750$ as in (Wang et al., 2015a) for $100$ epochs. Kernel widths for FKCCA/NKCCA, learning rate and momentum for DCCA, kernel widths and neighborhood sizes for NCCA/PLCCA are selected by grid search based on total tuning set correlation. Sensitivity to their values is mild over a large range; e.g., setting the kernel widths to 30-60% of the sample $L_{2}$ norm gives similarly good results. For NCCA/PLCCA, input dimensionalities are first reduced by PCA to 20% of the original ones (except that PLCCA does not apply PCA for view 2 in order to extract a $112D$ projection). The total correlation achieved by each algorithm is given in Table 1. We also report the running time (in seconds) of the algorithms (measured with a single thread on a workstation with a 3.2GHz CPU and 56G main memory), each using its optimal hyperparameters, and including the time for exact $15$-nearest neighbor search for NCCA/PLCCA. Overall, NCCA achieves the best canonical correlation while being much faster than the other nonlinear methods.

We now demonstrate the algorithms on a noisy MNIST dataset, generated identically to that of Wang et al. (2015b) but with a larger training set. View 1 inputs are randomly rotated images ($28\times 28$, gray scale) from the original MNIST dataset (LeCun et al., 1998), and the corresponding view 2 inputs are randomly chosen images with the same identity plus additive uniform pixel noise. We generate $450K$/$10K$/$10K$ pairs of images for training/tuning/testing (Wang et al. (2015b) uses a $50K$-pair training set). This dataset satisfies the multi-view assumption that given the label, the views are uncorrelated, so that the most correlated subspaces should retain class information and exclude the noise. Following Wang et al. (2015b), we extract a low-dimensional projection of the view 1 images with each algorithm, run spectral clustering to partition the splits into $10$ classes (with clustering parameters tuned as in (Wang et al., 2015b)), and compare the clustering with ground-truth labels and report the clustering accuracy. We also train a one-vs.-one linear SVM (Chang and Lin, 2011) on the projections with highest cluster accuracy for each algorithm (we reveal labels of 10% of the training set for fast SVM training) and report the classification error rates. The tuning procedure is as for XRMB except that we now select the projection dimensionality from $\{10,20,30\}$. For NCCA/PLCCA we first reduce dimensionality to $100$ by PCA for density estimation and exact nearest neighbor search, and use a randomized algorithm (Halko et al., 2011) to compute the SVD of the $450K\times 450K$ matrix $\SS$; for RKCCA/NKCCA we use an approximation rank of $5000$; for DCCA we use $3$ ReLU hidden layers of $1500$ units in each view and train with stochastic optimization of minibatch size $4500$. Clustering and classification results on the original 784D view 1 inputs are recorded as the baseline. Table 2 shows the clustering accuracy and classification error rates on the test set, as well as training run times, and Figure 3 shows t-SNE embeddings (van der Maaten and Hinton, 2008) of several algorithms with their optimal hyper-parameters. NCCA and DCCA achieve near perfect class separation.

Several points are worth noting regarding the experiments. First, the computation for NCCA and PLCCA is dominated by the exact kNN search; approximate search (Arya et al., 1998; Andoni and Indyk, 2006) should make NCCA/PLCCA much more efficient. Second, we have not explored the space of choices for density estimates; alternative choices, such as adaptive KDE (Terrell and Scott, 1992), could also further improve performance. Our current choice of KDE would seem to require large training sets for high-dimensional problems. Indeed, with less training data we do observe a drop in performance, but NCCA still outperforms KCCA; for example, using a 50K subset of the MNIST training set—an order of magnitude less data—the classification error rates when using FKCCA/NKCCA/DCCA/NCCA are 5.9/5.2/2.9/4.7%.

Conclusion

We have presented closed-form solutions to the nonparametric CCA (NCCA) and partially linear CCA (PLCCA) problems. As opposed to kernel CCA, which restricts the nonparametric projections to lie in a predefined RKHS, we have addressed the unconstrained setting. We have shown that the optimal nonparametric projections can be obtained from the SVD of a kernel defined via the pointwise mutual information between the views. This leads to a simple algorithm that outperforms KCCA and matches deep CCA on multiple datasets, while being more computationally efficient than either for moderate-sized data sets. Future work includes leveraging approximate nearest neighbor search and alternative density estimates.

Appendix A Proof of Lemma 3.1

Our objective is the sum of correlations in all $L$ dimensions. Let us consider the correlation in the $j$th dimension. From the Cauchy-Schwartz inequality, we have

Thanks to Nathan Srebro, Ryota Tomioka, and Yochai Blau for fruitful discussions. This research was supported by NSF grant IIS-1321015. The opinions expressed in this work are those of the authors and do not necessarily reflect the views of the funding agency.

References