# A Jacobi-Davidson Method for Large Scale Canonical Correlation Analysis  
[PERSON]  
1 College of Computer and Information Science, Fujian Agriculture and Forestry University,  
Fuzhou 350002, China; [EMAIL_ADDRESS]  
\[\max\ \operatorname{trace}(X^{\mathbb{T}}CY),\quad\text{subject to}\quad X^{ \mathbb{T}}(A+\kappa_{a}I_{m})X=I_{k}\text{ and }Y^{\mathbb{T}}(B+\kappa_{b}I_{n})Y=I_{k},\]  
where \(\kappa_{a}\) and \(\kappa_{b}\) are called regularization parameters and they usually are chosen to maximize the cross-validation score [11]. In other words, \(A\) and \(B\) are replaced by \(A+\kappa_{a}I_{m}\) and \(B+\kappa_{b}I_{n}\) to keep the invertible of \(A\) and \(B\), respectively. Hence, in this paper, by default we assume \(A\) and \(B\) are both positive definite and \(m\geq n\) unless explicitly stated otherwise. As shown in [4], the optimization problem (3) can be equivalently transformed into solving the following generalized eigenvalue problem  
\[Kz:= \begin{bmatrix}0&C\\ C^{\mathbb{T}}&0\end{bmatrix}\begin{bmatrix}x\\ y\end{bmatrix} =\lambda \begin{bmatrix}A&0\\ 0&B\end{bmatrix}\begin{bmatrix}x\\ y\end{bmatrix} =\lambda Mz, \tag{4}\]  
where the positive definiteness of the matrices \(A\) and \(B\) implies \(M\) being positive definite. The generalized eigenvalue problem (4) is referred as the Canonical Correlation Generalized Eigenvalue Problem (CCGEP) in this paper. Define \(J:= \begin{bmatrix}R_{\lambda}^{-1}&0\\ 0&R_{\text{s}}^{-1}\end{bmatrix} \) where  
\[A=R_{\lambda}^{\mathbb{T}}R_{\lambda}\quad\text{and}\quad B=R_{\text{s}}^{ \mathbb{T}}R_{\text{s}} \tag{5}\]  
are their Cholesky decomposition. It is easy to verify that  
\[J^{\mathbb{T}} \begin{bmatrix}0&C\\ C^{\mathbb{T}}&0\end{bmatrix} JJ^{-1} \begin{bmatrix}x\\ y\end{bmatrix} =\lambda J^{\mathbb{T}} \begin{bmatrix}A&0\\ 0&B\end{bmatrix} JJ^{-1} \begin{bmatrix}x\\ y\end{bmatrix} \]  
gives rise to  
\[ \begin{bmatrix}0&\widetilde{C}\\ \widetilde{C}^{\mathbb{T}}&0\end{bmatrix}\begin{bmatrix}p\\ q\end{bmatrix} =\lambda \begin{bmatrix}p\\ q\end{bmatrix} , \tag{6}\]  
where  
\[\widetilde{C}=R_{\lambda}^{-\mathbb{T}}CR_{\text{s}}^{-1},\quad p=R_{\lambda} x\quad\text{and}\quad q=R_{\text{s}}y, \tag{7}\]  
and it implies  
\[\widetilde{C}q=\lambda p\quad\text{and}\quad\widetilde{C}^{\mathbb{T}}p= \lambda q.\]  
It means that the eigenpairs of (4) can be obtained by computing the singular values and the associated left and right singular vectors of \(\widetilde{C}\). This method works well when the sample size \(d\) and feature dimension \(m\) and \(n\) are of moderate size but it will be very slow and numerically unstable for large-scale datasets which are ubiquitous in the age of \"Big Data\" [12]. For large-scale datasets, the equivalence between (4) and (6) makes it possible for us to simply adapt the subspace type algorithms for calculating the partial singular values decomposition, such as Lanczos type algorithms [13; 14] and Davidson type algorithms [15; 16], and then translate them for CCGEP (4). However, in practice, the decompositions of the sample covariance matrices \(A\) and \(B\) are usually unavailable in large scale matrix cases. The reason is that the decomposition is too expensive to compute explicitly for large scale problems, and may destroy the sparsity and some structural information. Furthermore, sometime sample covariance matrices \(A\) and \(B\) should never be explicitly formed, such as in online learning systems. Meanwhile, in [17], it is suggested to solve CCGEP (4) by considering the large scale symmetric positive definite pencil \(\{K,M\}\). Some subspace type numerical algorithms also have been generalized to computing partial eigenpairs of \(\{K,M\}\), see [18; 19]. However, these generic algorithms do not make use of the special structure in (4), and they usually are less efficient than custom-made algorithms. Therefore, existing algorithms either can not avoid the covariance matrices decomposition, or do not consider the structure of CCGEP.