Finding Linear Structure in Large Datasets with Scalable Canonical Correlation Analysis

Introduction

Canonical Correlation Analysis (CCA), first introduced in 1936 by (Hotelling, 1936), is a foundamental statistical tool to characterize the relationship between two multidimensional variables, which finds a wide range of applications. For example, CCA naturally fits into multi-view learning tasks and tailored to generate low dimensional feature representations using abandunt and inexpensive unlabeled datasets to supplement or refine the expensive labeled data in a semi-supervised fashion. Improved generalization accuracy has been witnessed or proved in areas such as regression (Kakade & Foster, 2007), clustering (Chaudhuri et al., 2009; Blaschko & Lampert, 2008), dimension reduction (Foster et al., 2008; McWilliams et al., 2013), word embeddings (Dhillon et al., 2011, 2012), etc. Besides, CCA has also been succesfully applied to genome-wide association study (GWAS) and has been shown powerful for understanding the relationship between genetic variations and phenotypes (Witten et al., 2009; Chen et al., 2012).

There are various equivalent ways to define CCA and here we use the linear algebraic formulation of (Golub & Zha, 1995), which captures the very essense of the procedure, pursuing the directions of maximal correlations between two data matrices.

Let ${\bf S_{x}^{-\frac{1}{2}}}{\bf S_{xy}}{\bf S_{y}^{-\frac{1}{2}}}={\bf U}{\bf D}{\bf V}^{\top}$ be the singular value decomposition. Then ${\bf\Phi}={\bf S_{x}^{-\frac{1}{2}}}{\bf U}$, ${\bf\Psi}={\bf S_{y}^{-\frac{1}{2}}}{\bf V}$, and ${\bf\Lambda}={\bf D}$ where ${\bf\Phi}=(\phi_{1},\cdots,\phi_{p}),{\bf\Psi}=(\psi_{1},\cdots,\psi_{p})$ and ${\bf\Lambda}=diag(\lambda_{1},\cdots,\lambda_{p})$.

The identifiability of canonical vectors $({\bf\Phi},{\bf\Psi})$ is equivalent to the identifiability of the singular vectors $({\bf U},{\bf V})$. Lemma 1.1 implies that the leading $k$ dimensional CCA subspace can be solved by first computing the whitening matrices ${\bf S_{x}^{-\frac{1}{2}}},{\bf S_{y}^{-\frac{1}{2}}}$ and then perform a $k$-truncated SVD on the whitened covariance matrix ${\bf S_{x}^{-\frac{1}{2}}}{\bf S_{xy}}{\bf S_{y}^{-\frac{1}{2}}}$. This classical algorithm is feasible and accurate when the data matrices are small but it can be slow and numerically unstable for large scale datasets which are common in modern natural language processing (large corpora, Dhillon et al. (2011, 2012)) and multi-view learning (abandunt and inexpensive unlabeled data, Hariharan & Subramanian (2014)) applications.

Throughout the paper, we call the step of orthonormalizing the columns of ${\bf X}$ and ${\bf Y}$ whitening step. The computational complexity of the classical algorithm is dominated by the whitening step. There are two major bottlenecks,

Huge matrix multiplication ${\bf X}^{\top}{\bf X},{\bf Y}^{\top}{\bf Y}$ to obtain ${\bf S_{x}},{\bf S_{y}}$ with computational complexity $O(np_{1}^{2}+np_{2}^{2})$ for general dense ${\bf X}$ and ${\bf Y}$.

Large matrix decomposition to compute ${\bf S_{x}^{-\frac{1}{2}}}$ and ${\bf S_{y}^{-\frac{1}{2}}}$ with computational complexity $O(p_{1}^{3}+p_{2}^{3})$ (Even when ${\bf X}$ and ${\bf Y}$ are sparse, ${\bf S_{x}},{\bf S_{y}}$ are not necessarily sparse)

The whitening step dominates the $k$-truncated SVD step because the top $k$ dimensional singular vectors can be efficiently computed by randomized SVD algorithms (see Halko et al. (2011) and many others).

Another classical algorithm (built-in function in Matlab) introduced in (Björck & Golub, 1973) uses a different way of whitening. It first carrys out a QR decomposition, ${\bf X}={\bf Q}_{x}{\bf R}_{x}$ and ${\bf Y}={\bf Q}_{y}{\bf R}_{y}$ and then performs a SVD on ${\bf Q}_{x}^{\top}{\bf Q}_{y}$, which has the same computational complexity $O(np_{1}^{2}+np_{2}^{2})$ as the algorithm indicated by Lemma 1.1. However, it is difficult to exploit sparsity in QR factorization while ${\bf X}^{\top}{\bf X},{\bf Y}^{\top}{\bf Y}$ can be efficiently computed when ${\bf X}$ and ${\bf Y}$ are sparse.

Besides computational issues, extra $O(p_{1}^{2}+p_{2}^{2})$ space is necessary to store two whitening matrices ${\bf S_{x}^{-\frac{1}{2}}}$ and ${\bf S_{y}^{-\frac{1}{2}}}$ (typically dense). In high dimensional applications where the number of features is huge, this can be another bottleneck considering the capacity of RAM of personal desktops (10-20 GB). In large distributed storage systems, the extra required space might incur heavy communication cost.

Therefore, it is natural to ask: is there a scalable algorithm that avoids huge matrix decomposition and huge matrix multiplication? Is it memory efficient? Or even more ambitiously, is there an online algorithm that generates decent approximation given a fixed computational power (e.g. CPU time, FLOP)?

2 Related Work

Scalability begins to play an increasingly important role in modern machine learning applications and draws more and more attention. Recently lots of promising progress emerged in the literature concerning with randomized algorithms for large scale matrix approximations, SVD, and Principal Component Analysis (Sarlos, 2006; Liberty et al., 2007; Woolfe et al., 2008; Halko et al., 2011). Unfortunately, these techniques does not directly solve CCA due to the whitening step. Several authors have tried to devise a scalable CCA algorithm. Avron et al. (2013) proposed an efficient approach for CCA between two tall and thin matrices ($p_{1},p_{2}\ll n$) harnessing the recently developed tools, Subsampled Randomized Hadamard Transform, which only subsampled a small proportion of the $n$ data points to approximate the matrix product. However, when the size of the features, $p_{1}$ and $p_{2}$, are large, the sampling scheme does not work. Later, Lu & Foster (2014) consider sparse design matrices and formulate CCA as iterative least squares, where in each iteration a fast regression algorithm that exploits sparsity is applied.

Another related line of research considers stochastic optimization algorithms for PCA (Arora et al., 2012; Mitliagkas et al., 2013; Balsubramani et al., 2013), which date back to Oja & Karhunen (1985). Compared with batch algorithms, the stochastic versions empirically converge much faster with similar accuracy. Further, these stochastic algorithms can be applied to streaming setting where data comes sequentially (one pass or several pass) without being stored. As mentioned in (Arora et al., 2012), stochastic optimization algorithm for CCA is more challenging and remains an open problem because of the whitening step.

3 Main Contribution

The main contribution of this paper is to directly tackle CCA as a nonconvex optimization problem and propose a novel Augmented Approximate Gradient (AppGrad) scheme and its stochastic variant for finding the top $k$ dimensional canonical subspace. Its advantages over state-of-art CCA algorithms are three folds. Firstly, AppGrad scheme only involves large matrix multiplying a thin matrix of width $k$ and small matrix decomposition of dimension $k\times k$, and therefore to some extent is free from the two bottlenecks. It also benefits if ${\bf X}$ and ${\bf Y}$ are sparse while classical algorithm still needs to invert the dense matrices ${\bf X}^{\top}{\bf X}$ and ${\bf Y}^{\top}{\bf Y}$. Secondly, AppGrad achieves optimal storage complexity $O(k(p_{1}+p_{2}))$, the space necessary to store the output, compared with classical algorithms which usually require $O(p_{1}^{2}+p_{2}^{2})$ for storing the whitening matrices. Thirdly, the stochastic (online) variant of AppGrad is especially efficient for large scale datasets if moderate accuracy is desired. It is well-suited to the case when computational resources are limited or data comes as a stream. To the best of our knowledge, it is the first stochastic algorithm for CCA, which partly gives an affirmative answer to a question left open in (Arora et al., 2012).

The rest of the paper is organized as follows. We introduce AppGrad scheme and establish its convergence properties in section 2. We extend the algorithm to stochastic settings in section 3. Extensive real data experiments are presented in section 4. Concluding remarks and future work are summarized in section 5. Proof of Theorem 2.1 and Proposition 2.3 are relegated to the supplementary material.

Algorithm

For simplicity, we first focus on the leading canonical pair $(\phi_{1},\psi_{1})$ to motivate the proposed algorithms. Results for general scenario can be obtained in the same manner and will be briefly discussed in the later part of this section.

To begin with, we recast CCA as an nonconvex optimization problem (Golub & Zha, 1995).

Although (1) is a nonconvex (due to the nonconvex constraint), (Golub & Zha, 1995) showed that an alternating minimization strategy (Algorithm 1), or rather iterative least squares, actually converges to the leading canonical pair. However, each update $\phi^{t+1}={\bf S_{x}^{-1}}{\bf S_{xy}}\psi^{t}$ is computationally intensive. Essentially, the alternating least squares acts like a second order method, which is usually recognized to be inefficient for large-scale datasets, especially when current estimate is not close enough to the optimum. Therefore, it is natural to ask: is there a valid first order method that solves (1)?

Heuristics borrowed from convex optimization literature give rise to a projected gradient scheme summarized in Algorithm 2. Instead of completely solving a least squares in each iterate, a single gradient step of (1) is performed and then project back to the constrained domain, which avoids inverting a huge matrix. Unfortunately, the following proposition demonstrates that Algorithm 2 fails to converge to the leading canonical pair.

If leading canonical correlation $\lambda_{1}\neq 1$ and either $\phi_{1}$ is not an eigenvector of ${\bf S_{x}}$ or $\psi_{1}$ is not an eigenvector of ${\bf S_{y}}$, then $\forall\eta_{1},\eta_{2}>0$, the leading canonical pair $(\phi_{1},\psi_{1})$ is not a fixed point of the naive gradient scheme in Algorithm 2. Therefore, the algorithm does not converge to $(\phi_{1},\psi_{1})$.

The proof is similar to the proof of Proposition 2.2 and we leave out the details here. ∎

The failure of Algorithm 2 is due to the nonconvex nature of (1). Although every gradient step might decrease the objective function, this property no longer persists after projecting to its nonconvex domain \big{\{}(\phi,\psi)\,|\,\phi^{\top}{\bf S_{x}}\phi=1,\,\psi^{\top}{\bf S_{y}}\psi=1\big{\}} (the normalization step). On the contrary, decreases triggered by gradient descent is always maintained if projecting to a convex region.

2 AppGrad Scheme

As a remedy, we propose a novel Augmented Approximate Gradient (AppGrad) scheme summarized in Algorithm 3. It inherits the convergence guarantee of alternating least squares as well as the scalability and memory efficiency of first order methods, which only involves matrix-vector multiplication and only requires $O(p_{1}+p_{2})$ extra space.

AppGrad seems unnatural at first sight but has some nice intuitions behind as we will discuss later. The differences and similarities between these algorithms are subtle but crucial. Compared with the naive gradient descent, we introduce two auxiliary variables $(\widetilde{\phi}^{t},\widetilde{\psi}^{t})$, an unnormalized version of $(\phi^{t},\psi^{t})$. During each iterate, we keep updating $\widetilde{\phi}^{t}$ and $\widetilde{\psi}^{t}$ without scaling them to have unit norm, which in turn produces the ‘correct’ normalized counterpart, $(\phi^{t},\psi^{t})$. It turns out that $(\phi_{1},\psi_{1},\lambda_{1}\phi_{1},\lambda_{1}\psi_{1})$ is a fixed point of the dynamic system $\{(\phi^{t},\psi^{t},\widetilde{\phi}^{t},\widetilde{\psi}^{t})\}_{t=0}^{\infty}$.

$\forall\,i\leq p$, let $\widetilde{\phi}_{i}=\lambda_{i}\phi_{i},\widetilde{\psi}_{i}=\lambda_{i}\psi_{i}$, then $(\phi_{i},\psi_{i},\widetilde{\phi}_{i},\widetilde{\psi}_{i})$ are the fixed points of AppGrad scheme.

To prove the proposition, we need the following lemma that characterizes the relations among some key quantities.

${\bf S_{xy}}={\bf S_{x}}{\bf\Phi}{\bf\Lambda}{\bf\Psi}^{\top}{\bf S_{y}}$

By Lemma 1.1, $\bf{S_{x}^{-\frac{1}{2}}S_{xy}S_{y}^{-\frac{1}{2}}}={\bf U}{\bf D}{\bf V}^{\top}$, where ${\bf U}={\bf S_{x}^{\frac{1}{2}}}{\bf\Phi}$, ${\bf V}={\bf S_{y}^{\frac{1}{2}}}{\bf\Psi}$ and ${\bf D}={\bf\Lambda}$. Then we have ${\bf S_{xy}}={\bf S_{x}^{\frac{1}{2}}}{\bf U}{\bf D}{\bf V}^{\top}{\bf S_{y}^{\frac{1}{2}}}={\bf S_{x}}{\bf\Phi}{\bf\Lambda}{\bf\Psi}^{\top}{\bf S_{y}}$. ∎

Substitute $(\phi^{t},\psi^{t},\widetilde{\phi}^{t},\widetilde{\psi}^{t})=(\phi_{i},\psi_{i},\widetilde{\phi}_{i},\widetilde{\psi}_{i})$ into the iterative formula in Algorithm 3.

The second equality is direct application of Lemma 2.2. The third equality is due to the fact that ${\bf\Psi}^{\top}{\bf S_{y}}{\bf\Psi}=I_{p}$. Then,

Therefore $(\widetilde{\phi}^{t+1},\phi^{t+1})=(\widetilde{\phi}^{t},\phi^{t})=(\widetilde{\phi}_{i},\phi_{i})$. A symmetric argument will show that $(\widetilde{\psi}^{t+1},\psi^{t+1})=(\widetilde{\psi}^{t},\psi^{t})=(\widetilde{\psi}_{i},\psi_{i})$, which completes the proof. ∎

The connection between AppGrad and alternating minimization strategy is not instaneous. Intuitively, when $(\phi^{t},\psi^{t})$ is not close to $(\phi_{1},\psi_{1})$, solving the least squares completely as carried out in Algorithm 1 is a waste of computational power (informally by regarding it as a second order method, the Newton Step has fast convergence only when current estimate is close to the optimum). Instead of solving a sequence of possibly unrelevant least squares, the following lemma shows that AppGrad directly targets at the least squares that involves the leading canonical pair.

Let $(\phi_{1},\psi_{1})$ be the leading canonical pair and $(\widetilde{\phi}_{1},\widetilde{\psi}_{1})=\lambda_{1}(\phi_{1},\psi_{1})$. Then,

Let $\phi^{*}=\arg\min\limits_{\phi}\frac{1}{2n}\|{\bf X}\phi-{\bf Y}\psi_{1}\|^{2}$, by optimality condition, ${\bf S_{x}}\phi^{*}={\bf S_{xy}}\psi_{1}$. Apply Lemma 2.2,

Similar argument gives $\psi^{*}=\widetilde{\psi}_{1}$ ∎

Lemma 2 characterizes the relationship between leading canonical pair $(\phi_{1},\psi_{1})$ and its unnormalized counterpart $(\widetilde{\phi}_{1},\widetilde{\psi}_{1})$, which sheds some insight on how AppGrad works. The intuition is that $(\phi^{t},\psi^{t})$ and $(\widetilde{\phi}^{t},\widetilde{\psi}^{t})$ are current estimations of $(\phi_{1},\psi_{1})$ and $(\widetilde{\phi}_{1},\widetilde{\psi}_{1})$, and the updates of $(\widetilde{\phi}^{t+1},\widetilde{\psi}^{t+1})$ in Algorithm 3 are actually gradient steps of the least squares in (2), with the unknown truth $(\phi_{1},\psi_{1})$ approximated by $(\phi^{t},\psi^{t})$. In terms of mathematics,

The normalization step in Algorithm 3 corresponds to generating new approximations of $(\phi_{1},\psi_{1})$, namely $(\phi^{t+1},\psi^{t+1})$, using the updated $(\widetilde{\phi}^{t+1},\widetilde{\psi}^{t+1})$ through the relationship $(\phi_{1},\psi_{1})=(\widetilde{\phi}_{1}/\|\widetilde{\phi}_{1}\|_{x},\,\widetilde{\psi}_{1}/\|\widetilde{\psi}_{1}\|_{y})$. Therefore, one can interpret AppGrad as approximate gradient scheme for solving (2). When $(\widetilde{\phi}^{t},\widetilde{\psi}^{t})$ converge to $(\widetilde{\phi}_{1},\widetilde{\psi}_{1})$, its scaled version $(\phi^{t},\psi^{t})$ converge to the leading canonical pair $(\phi_{1},\psi_{1})$.

The following theorem shows that when the estimates enter a neighborhood of the true canonical pair, AppGrad is contractive. Define the error metric $e_{t}=\|\Delta\widetilde{\phi}^{t}\|^{2}+\|\Delta\widetilde{\psi}^{t}\|^{2}$ where $\Delta\widetilde{\phi}^{t}=\widetilde{\phi}^{t}-\widetilde{\phi}_{1},\Delta\widetilde{\psi}^{t}=\widetilde{\psi}^{t}-\widetilde{\psi}_{1}$.

where \delta=1-\Big{(}1-\frac{2(\lambda_{1}^{2}-\lambda_{2}^{2})-L_{1}e_{0}}{2\lambda_{1}^{2}}\Big{)}^{\frac{1}{2}}>0

The theorem reveals that the larger is the eigengap $\lambda_{1}-\lambda_{2}$, the broader is the contraction region. We didn’t try to optimize the conditions above and empirically as shown in the experiments, a randomized initialization always suffices to capture most of the correlations.

3 General Rank-k𝑘k Case

Following the spirit of rank-one case, AppGrad can be easily generalized to compute the top $k$ dimesional canonical subspace as summarized in Algorithm 4. The only difference is that the original scalar normalization is replaced by its matrix counterpart, that is to multiply the inverse of the square root matrix ${\bf\Phi}^{t+1}=\widetilde{{\bf\Phi}}^{t+1}{\bf U}_{x}{\bf D}_{x}^{-\frac{1}{2}}{\bf U}_{x}^{\top}$, ensuring that $({\bf\Phi}^{t+1})^{\top}{\bf X}^{\top}{\bf X}{\bf\Phi}^{t+1}={\bf I}_{k}$.

Notice that the gradient step only involves a large matrix multiplying a thin matrix of width $k$ and the SVD is performed on a small $k\times k$ matrix. Therefore, the computational complexity per iteration is dominated by the gradient step, of order $O(n(p_{1}+p_{2})k)$. The cost will be further reduced when the data matrices ${\bf X},{\bf Y}$ are sparse.

Compared with classical spectral agorithm which first whitens the data matrices and then performs a SVD on the whitened covariance matrix, AppGrad actually merges these two steps together. This is the key of its efficiency. In a high level, whitening the whole data matrix is not necessary and we only want to whiten the directions that contain the leading CCA subspace. However, these directions are unknown and therefore for two-step procedures, whitening the whole data matrix is unavoidable. Instead, AppGrad tries to identify (gradient step) and whiten (normalization step) these directions simultaneously. In this way, every normalization step is only performed on the potential $k$ dimensional target CCA subspace and therefore only deals with a small $k\times k$ matrix.

Parallel results of Lemma 1, Proposition 2.1, Proposition 2.2, Lemma 2 for this general scenario can be established in a similar manner. Here, to make Algorithm 4 more clear, we state the fixed point result of which the proof is similar to Proposition 2.2.

Let ${\bf\Lambda}_{k}=diag(\lambda_{1},\cdots,\lambda_{k})$ be the diagonal matrix of top $k$ canonical correlations and let ${\bf\Phi}_{k}=(\phi_{1},\cdots,\phi_{k}),{\bf\Psi}_{k}=(\phi_{1},\cdots,\phi_{k})$ be the top $k$ CCA vectors. Also denote $\widetilde{{\bf\Phi}}_{k}={\bf\Phi}_{k}{\bf\Lambda}_{k}$ and $\widetilde{{\bf\Psi}}_{k}={\bf\Psi}_{k}{\bf\Lambda}_{k}$. Then for any $k\times k$ orthogonal matrix ${\bf Q}$, $({\bf\Phi}_{k},{\bf\Psi}_{k},\widetilde{{\bf\Phi}}_{k},\widetilde{{\bf\Psi}}_{k}){\bf Q}$ is a fixed point of AppGrad scheme.

The top $k$ dimensional canonical subspace is identifiable up to a rotation matrix and Proposition 2.3 shows that every optimum is a fixed point of AppGrad scheme.

4 Kernelization

Following the same logic as Proposition 2.3, a similar fixed point result can be proved. Therefore, Algorithm 4 can be directly applied to compute ${\bf W}_{\mathcal{X}},{\bf W}_{\mathcal{Y}}$ by simply replacing ${\bf X},{\bf Y}$ with ${\bf K}_{\mathcal{X}},{\bf K}_{\mathcal{Y}}$.

Stochastic AppGrad

Recently, there is a growing interest in stochastic optimization which is shown to have better performance for large-scale learning problems (Bousquet & Bottou, 2008; Bottou, 2010). Especially in the so-called ‘data laden regime’, where data is abundant and the bottleneck is runtime, stochastic optimization dominate batch algorithms both empirically and theoretically. Given these advantages, lots of efforts have been spent on developing stochastic algorithms for principal component analysis (Oja & Karhunen, 1985; Arora et al., 2012; Mitliagkas et al., 2013; Balsubramani et al., 2013). Despite promising progress in PCA, as mentioned in (Arora et al., 2012), stochastic CCA is more challenging and remains an open problem due to the whitening step.

As a gradient scheme, AppGrad naturally generalizes to the stochastic regime and we summarize in Algorithm 5. Compared with the batch version, only a small subset of samples are used to compute the gradient, which reduces the computational cost per iteration from $O(n(p_{1}+p_{2})k)$ to $O(m(p_{1}+p_{2})k)$ ($m=|\mathcal{I}|$ is the size of the minibatch). Empirically, this makes stochastic AppGrad much faster than the batch version as we will see in the experiments. Also, for large scale applications when fully calculating the CCA subspace is prohibitive, stochastic AppGrad can generate a decent approximation given a fixed computational power, while other algorithms only give a one-shot estimate after the whole procedure is carried out completely. Moreover, when there is a generative model, as shown in (Bousquet & Bottou, 2008), due to the tradeoff between statistical and numerical accuracy, fully solving an empirical risk minimization is unnecessary since the statistical error will finally dominate. On the contrary, stochastic optimization directly tackles the problem in the population level and therefore is more statistically efficient.

It is worth mentioning that the normalization step is accomplished using a sampled Gram matrix $\frac{1}{m}{\bf X}_{\mathcal{I}}^{\top}{\bf X}_{\mathcal{I}}$ and $\frac{1}{m}{\bf Y}_{\mathcal{I}}^{\top}{\bf Y}_{\mathcal{I}}$. A key observation is that when $m\in O(k)$, $(\widetilde{{\bf\Phi}}^{t+1})^{\top}(\frac{1}{m}{\bf X}_{\mathcal{I}}^{\top}{\bf X}_{\mathcal{I}})\widetilde{{\bf\Phi}}^{t+1}\approx(\widetilde{{\bf\Phi}}^{t+1})^{\top}(\frac{1}{m}{\bf X}^{\top}{\bf X})\widetilde{{\bf\Phi}}^{t+1}$ using standard concentration inequality, because the matrix we want to approximate $(\widetilde{{\bf\Phi}}^{t+1})^{\top}(\frac{1}{m}{\bf X}^{\top}{\bf X})\widetilde{{\bf\Phi}}^{t+1}$ is a $k\times k$ matrix, while generally $O(p)$ sample is needed to have $\frac{1}{m}{\bf X}_{\mathcal{I}}^{\top}{\bf X}_{\mathcal{I}}\approx\frac{1}{n}{\bf X}^{\top}{\bf X}$. As we have argued in previous section, this bonus is a byproduct of the fact that AppGrad tries to identify and whiten the directions that contains the CCA subspace simultaneously, or else $O(p)$ samples are necessary for whitening the whole data matrices.

Experiments

In this section, we present experiments on four real datasets to evaluate the effectiveness of the proposed algorithms for computing the top 20 ($k$=20) dimensional canonical subspace. A short summary of the datasets is in Table 1.

Mediamill is an annotated video dataset from the Mediamill Challenge (Snoek et al., 2006). Each image is a representative keyframe of a video shot annotated with 101 labels and consists of 120 features. CCA is performed to explore the correlation structure between the images and its labels.

MNIST is a database of handwritten digits. CCA is used to learn correlated representations between the left and right halves of the images.

Penn Tree Bank dataset is extracted from Wall Street Journal, which consists of $1.17$ million tokens and a vocabulary size of $43,000$ (Lamar et al., 2010). CCA has been successfully used on this dataset to build low dimensional word embeddings (Dhillon et al., 2011, 2012). The task here is a CCA between words and their context. We only consider the 10, 000 most frequent words to avoid sample sparsity.

URL Reputation dataset (Ma et al., 2009) is extracted from UCI machine learning repository. The dataset contains 2.4 million URLs each represented by 3.2 million features. For simplicity we only use the first 2 million samples. $38\%$ of the features are host based features like WHOIS info, IP prefix and $62\%$ are lexical based features like Hostname and Primary domain. We run a CCA between a subset of host based features and a subset of lexical based features.

Evaluation Criterion: The evaluation criterion we use for the first three datasets (Mediamill, MNIST, Penn Tree Bank) is Proportions of Correlations Captured (PCC). To introduce this term, we first define Total Correlations Captured (TCC) between two matrices to be the sum of their canonical correlations as defined in Lemma 1.1. Then, for estimated top $k$ dimensional canonical subspace $\widehat{{\bf\Phi}}_{k},\widehat{{\bf\Psi}}_{k}$ and true leading $k$ dimensional CCA subspace ${\bf\Phi}_{k},{\bf\Psi}_{k}$, PCC is defined as

Intuitively PCC characterizes the proportion of correlations captured by certain algorithm compared with the true CCA subspace. Therefore, the higher is PCC the better is the estimated CCA subspace. This criterion has two major advantages over subspace distance $\|P_{\widehat{{\bf\Phi}}_{k}}-P_{{\bf\Phi}_{k}}\|$ ($P_{\Omega}$ is projection matrix of the column space of $\Omega$). First, it is more natural and relevant considering that the goal of CCA is to capture most correlations between two data matrices. Second, when the eigengap $\Delta\lambda=\lambda_{k}-\lambda_{k+1}$ is not big enough, the top $k$ dimensional CCA subspace is ill posed while the correlations captured is well defined.

We use the output of the standard spectral algorithms as the truth $({\bf\Phi}_{k},{\bf\Psi}_{k})$ to calculate the denominator of PCC. However, for URL Reputation dataset, the number of samples and features are too large for the algorithm to compute the true CCA subspace in a reasonable amount of time and instead we only compare the numerator $\mbox{TCC}({\bf X}\widehat{{\bf\Phi}}_{k},{\bf Y}\widehat{{\bf\Psi}}_{k})$ (monotone w.r.t. PCC) for different algorithms.

Initialization We initialize $({\bf\Phi}^{0},{\bf\Psi}^{0})$ by first drawing $i.i.d.$ samples from standard Gaussian distribution and then normalize such that $({\bf\Phi}^{0})^{\top}{\bf S_{x}}{\bf\Phi}^{0}=I_{k}$ and $({\bf\Psi}^{0})^{\top}{\bf S_{y}}{\bf\Psi}^{0}=I_{k}$

Stepsize For both AppGrad and stochastic AppGrad, a small part of the training set is held out and cross-validation is used to choose the step size adaptively.

Regularization For all the algorithms, a little regularization is added for numerical stability which means we replace Gram matrix ${\bf X}^{\top}{\bf X}$ with ${\bf X}^{\top}{\bf X}+\lambda{\bf I}$ for some small positive $\lambda$.

Oversampling Oversampling means when aiming for top $k$ dimensional subspace, people usually computes top $k+l$ dimesional subspace from which a best $k$ diemsional subspace is extracted. In practice, $l=5\sim 10$ suffices to improve the performance. We only do a oversampling of 5 in the URL dataset.

2 Summary of Results

For the first three datasets (Mediamill, MNIST, Penn Tree Bank), both in-sample and out-of-sample PCC are computed for AppGrad and Stochastic AppGrad as summarized in Figure1. As you can see, both algorithms nearly capture most of the correlations compared with the true CCA subspace and stochastic AppGrad consistently achieves same PCC with much less computational cost than its batch version. Moreover, the larger is the size of the data, the bigger advantage will stochastic AppGrad obtain. One thing to notice is that, as revealed in Mediamill dataset, out-of-sample PCC is not necessarily less than in-sample PCC because both denominator and numerator will change on the hold out set.

For URL Reputation dataset, as we mentioned earlier, classical algorithms fails on a typical desktop. The reason is that these algorithms only produce a one-shot estimate after the whole procedure is completed, which is usually prohibitive for huge datasets. In this scenario, the advantage of online algorithms like stochastic AppGrad becomes crucial. Further, the stochastic nature makes the algorithm cost-effective and generate decent approximations given fixed computational resources (e.g. FLOP). As revealed by Figure 2, as the number of iterations increases, stochastic AppGrad captures more and more correlations.

Since the true CCA subspaces for URL dataset is too slow to compute, we compare our algorithm with some naive heuristics which can be carried out efficiently in large scale and catches a reasonable amount of correlation. Below is a brief description of them.

Non-Whitening (NW-CCA): directly perform SVD on the unwhitened covariance matrix ${\bf X}^{T}{\bf Y}$. This strategy is also used in (Witten et al., 2009)

Diagnoally Whitening (DW-CCA) (Lu & Foster, 2014): avoid inverting matrices by approximating ${\bf S_{x}^{-\frac{1}{2}}},{\bf S_{y}^{-\frac{1}{2}}}$ with $(\mbox{diag}({\bf S_{x}}))^{-\frac{1}{2}}$ and $(\mbox{diag}({\bf S_{y}}))^{-\frac{1}{2}}$.

Whitening the leading $m$ Principal Component Directions (PCA-CCA): First compute the leading $m$ dimensional principal component subspace and project the data matrices ${\bf X}$ and ${\bf Y}$ to the subspace, denote them ${\bf U}_{x}$ and ${\bf U}_{y}$. Then compute the top $k$ dimensional CCA subspace of the pair $({\bf U}_{x},{\bf U}_{y})$. At last, transform the CCA subspace of $({\bf U}_{x},{\bf U}_{y})$ back to the CCA subspace of orginal matrix pair $({\bf X},{\bf Y})$. Specifically for this example, we choose $m=1200$ (log(FLOP)=35, dominating the computational cost of Stochastic AppGrad) .

For all the heuristics mentioned above, SVD and PCA steps are carried out using the randomized algorithms in (Halko et al., 2011). For PCA-CCA, as the number of Principal Components ($m$) increases, more correlation will be captured but the computational cost will also increase. When $m=p$, PCA-CCA is reduced to the orginal CCA.

Essentially, all the heuristics are incorrect algorithms and try to approximately whiten the data matrices. As suggested by Figure 2, stochastic AppGrad significantly captures much more correlations.

Conclusions and Future Work

In this paper, we present a novel first order method, AppGrad, to tackle large scale CCA as a nonconvex optimization problem. This bottleneck-free algorithm is both memory efficient and computationally scalable. More importantly, its online variant is well-suited to practical high dimensional applications where batch algorithm is prohibitive and data laden regime where data is abundant and runtime is main concern.

Further, AppGrad is flexible and structure information can be easily incorporated into the algorithm. For example, if the canonical vectors are assumed to be sparse (Witten et al., 2009; Gao et al., 2014), a thresholding step can be added between the gradient step and normalization step to obtain sparse solutions while it is hard to add sparse constraint to the classical CCA formulation which is a generalized eigenvalue problem. Heuristics in (Witten et al., 2009) avoid this by simply skipping the whitening procedure (NW-CCA). (Gao et al., 2014) resorts to semidefinite programming and therefore is very slow. AppGrad with thresholding works well in simulations and we leave its theoretical properties for future research.

References

Proofs

A brief review of the notations in the main paper:

Further, we define $cos_{x}(u,v)=\frac{u^{\top}{\bf S_{x}}v}{{\!\|}u{\|_{x}\!}\,{\!\|}v{\|_{x}\!}}$, the cosine of the angle between two vectors induced by the inner product $\left\langle u,v\right\rangle=u^{\top}{\bf S_{x}}v$. Similarly, we define $cos_{y}(u,v)=\frac{u^{\top}{\bf S_{y}}v}{{\!\|}u{\|_{y}\!}\,{\!\|}v{\|_{y}\!}}$. To prove the theorem, we will repeatedly use the following lemma.

${\!\|}\Delta\phi^{t}{\|_{x}\!}\leq\frac{1}{\lambda_{1}}\sqrt{\frac{2}{1+cos_{x}(\phi^{t},\phi_{1})}}{\!\|}\Delta\widetilde{\phi}^{t}{\|_{x}\!}$ and ${\!\|}\Delta\psi^{t}{\|_{y}\!}\leq\frac{1}{\lambda_{1}}\sqrt{\frac{2}{1+cos_{y}(\psi^{t},\psi_{1})}}{\!\|}\Delta\widetilde{\psi}^{t}{\|_{y}\!}$

Proof of Lemma 6.1 Notice that $cos_{x}(\widetilde{\phi}^{t},\widetilde{\phi}_{1})=cos_{x}(\phi^{t},\phi_{1})$, then

Also notice that $\|\phi^{t}\|_{x}=\|\phi_{1}\|_{x}=1$, which implies $cos_{x}(\phi^{t},\phi_{1})=1-\|\phi^{t}-\phi_{1}\|_{x}^{2}/2=1-\|\Delta\phi^{t}\|_{x}^{2}/2$. Further

2 Proof of Theorem 2.1

Without loss of generality, we can always assume $cos_{x}(\widetilde{\phi}^{t},\widetilde{\phi}_{1}),cos_{y}(\widetilde{\psi}^{t},\widetilde{\psi}_{1})\geq 0$ because the canonical vectors are only identifiable up to a flip in sign and we can always choose $\widetilde{\phi}_{1},\widetilde{\psi}_{1}$ such that the cosines are nonnegative. Apply simple algebra to the gradient step $\widetilde{\phi}^{t+1}=\widetilde{\phi}^{t}-\eta({\bf S_{x}}\widetilde{\phi}^{t}-{\bf S_{xy}}\psi^{t})$, we have

By Lemma 2.2, $\eta({\bf S_{x}}\widetilde{\phi}_{1}-{\bf S_{xy}}\psi_{1})=\eta({\bf S_{x}}\widetilde{\phi}_{1}-\lambda_{1}{\bf S_{x}}\phi_{1})=0$, which implies

The last inequality uses the assumption that $\lambda_{max}({\bf S_{x}}),\lambda_{max}({\bf S_{y}})\leq L_{1}$. By Lemma6.1, $\|\Delta\psi^{t}\|_{y}\leq\frac{\sqrt{2}}{\lambda_{1}}\|\Delta\widetilde{\psi}^{t}\|_{y}$. Hence, $\|{\bf S_{xy}}\Delta\psi^{t}\|\leq\sqrt{2L_{1}}\|\Delta\widetilde{\psi}^{t}\|_{y}$. Also notice that $\|{\bf S_{x}}\Delta\widetilde{\phi}^{t}\|\leq\|{\bf S_{x}^{\frac{1}{2}}}\|\|{\bf S_{x}^{\frac{1}{2}}}\Delta\widetilde{\phi}^{t}\|\leq L_{1}^{\frac{1}{2}}\|\Delta\widetilde{\phi}^{t}\|_{x}$, then

Now, we are going to bound $(\Delta\widetilde{\phi}^{t})^{T}{\bf S_{xy}}\Delta\psi^{t}$. Because ${\bf S_{y}^{\frac{1}{2}}}{\bf\Psi}$ is an orthonormal matrix (orthogonal if $p=p_{1}$) and ${\bf S_{y}^{\frac{1}{2}}}\psi_{t}$ is a unit vector, there exisit coefficients $\alpha_{1},\cdots,\alpha_{p},\alpha_{\perp}$ and unit vector $\psi_{\perp}\in ColSpan({\bf S_{y}^{\frac{1}{2}}}{\bf\Psi})^{\perp}$ such that ${\bf S_{y}^{\frac{1}{2}}}\psi_{t}=\sum_{i=1}^{p}\alpha_{i}{\bf S_{y}^{\frac{1}{2}}}\psi_{i}+\alpha_{\perp}{\bf S_{y}^{\frac{1}{2}}}\psi_{\perp},\sum_{i=1}^{p}\alpha_{i}^{2}+\alpha_{\perp}^{2}=1$. Therefore,

By definition, $1-\alpha_{1}=1-cos_{y}(\psi^{t},\psi_{1})=\frac{\|\Delta\psi^{t}\|^{2}_{y}}{2}$. Further by Lemma 6.1,

Notice that $\sqrt{a}+\sqrt{b}\leq\sqrt{2(a+b)}$, we have

By definition, \delta=1-\frac{1}{\lambda_{1}}\Big{(}\frac{L_{1}}{2}\|\Delta\widetilde{\psi}^{0}\|^{2}+\frac{L_{1}}{2}\|\Delta\widetilde{\phi}^{0}\|^{2}+\lambda_{2}^{2}\Big{)}^{\frac{1}{2}} and $\eta=\frac{\delta}{6L_{1}}$. Substitute in (6) with $t=0$,

3 Proof of Proposition 2.3

Substitute $({\bf\Phi}^{t},{\bf\Psi}^{t},\widetilde{{\bf\Phi}}^{t},\widetilde{{\bf\Psi}}^{t})=({\bf\Phi}_{k},{\bf\Psi}_{k},{\bf\Phi}_{k}{\bf\Lambda}_{k},{\bf\Psi}_{k}{\bf\Lambda}_{k}){\bf Q}$ into the iterative formula in Algorithm 4.

The second equality is direct application of Lemma 1. The third equality is due to the fact that ${\bf\Psi}^{\top}{\bf S_{y}}{\bf\Psi}=I_{p}$. Then,

Therefore $({\bf\Phi}^{t+1},\widetilde{{\bf\Phi}}^{t+1})=({\bf\Phi}^{t},\widetilde{{\bf\Phi}}^{t})=({\bf\Phi}_{k},{\bf\Phi}_{k}{\bf\Lambda}_{k}){\bf Q}$. A symmetric argument will show that $({\bf\Psi}^{t+1},\widetilde{{\bf\Psi}}^{t+1})=({\bf\Psi}^{t},\widetilde{{\bf\Psi}}^{t})=({\bf\Psi}_{k},{\bf\Psi}_{k}{\bf\Lambda}_{k}){\bf Q}$, which completes the proof.