Year: 2021

Abstract:
Canonical Correlation Analysis (CCA) models are powerful for studying the associations between two sets of variables. The canonically correlated representations, termed canonical variates are widely used in unsupervised learning to analyze unlabeled multi-modal registered datasets. Despite their success, CCA models may break (or overfit) if the number of variables in either of the modalities exceeds the number of samples. Moreover, often a significant fraction of the variables measures modality-specific information, and thus removing them is beneficial for identifying the canonically correlated variates. Here, we propose `0-CCA, a method for learning correlated representations based on sparse subsets of variables from two observed modalities. Sparsity is obtained by multiplying the input variables by stochastic gates, whose parameters are learned together with the CCA weights via an `0-regularized correlation loss. We further propose `0-Deep CCA for solving the problem of non-linear sparse CCA by modeling the correlated representations using deep nets. We demonstrate the efficacy of the method using several synthetic and real examples. Most notably, by gating nuisance input variables, our approach improves the extracted representations compared to other linear, non-linear and sparse CCA-based models.

heading: 1 INTRODUCTION
text: Canonical Correlation Analysis (CCA) (Hotelling, 1936; Thompson, 2005), is a classic statistical method for finding the maximally correlated linear transformations of two modalities (or views). Using modalitiesX ∈ RDx×N andY ∈ RDy×N , which are centered and haveN samples withDx and Dy features, respectively, CCA seeks canonical vectors a ∈ RDx , and b ∈ RDy , such that ,u = aTX, and v = bTY , will maximize the sample correlations between the canonical variates, i.e.
max a, b6=0
ρ(aTX, bTY ) = aTXY T b
‖aTX‖2‖bTY ‖2 . (1)
To identify non-linear relations between input variables, several extensions of CCA have been proposed. Kernel methods such as Kernel CCA (Bach & Jordan, 2002), Non-parametric CCA (Michaeli et al., 2016) or Multi-view Diffusion maps (Lindenbaum et al., 2020; Salhov et al., 2020) learn the non-linear relations in reproducing Hilbert spaces. These methods have several shortcomings: they are limited to a designed kernel, and they have poor interpolation and extrapolation capabilities. Deep CCA (Andrew et al., 2013) overcomes these limitations by learning two non-linear transformations parametrized using neural networks. The model has several extension, see for example (Wang et al., 2016; Gundersen et al., 2019; Karami & Schuurmans, 2021).
CCA models have been widely used in biology (Pimentel et al., 2018), neuroscience (Al-Shargie et al., 2017), medicine (Zhang et al., 2017), and engineering (Chen et al., 2017), for unsupervised or semi-supervised learning. By extracting meaningful dimensionality reduced representations, CCA improves downstream tasks such as clustering, classification, or manifold learning. One key limitation of CCA is that it requires more samples than features, i.e., N > Dx, Dy . However, if we have more variables than samples, the estimation based on the closed-form solution of the CCA problem (in Eq. 1) breaks (Suo et al., 2017) (because the covariance matrix is low rank). Moreover, in high dimensional data, often some of the variables do not measure the phenomenon that is common to both modalities (therefore are not correlated) and thus should be omitted from the transformations. For these reasons, there has been a growing interest in studying sparse CCA models.
Sparse CCA (SCCA) models seek linear transformations which are based on a sparse subset of variables from the input modalities X and Y . Sparsifying the feature space improves interpretability, prevents overfitting, and removes the degeneracies inherent to N < Dx, Dy (once the number of selected variables is smaller thanN the covariance matrix becomes full rank). To encourage sparsity of the canonical vectors a and b, several authors (Wiesel et al., 2008; Cai et al., 2019) propose an `0 regularized variant of Eq. 1. However, these schemes are greedy and therefore may lead to suboptimal solutions. As demonstrated by (Waaijenborg et al., 2008; Parkhomenko et al., 2009; Witten et al., 2009; Hardoon & Shawe-Taylor, 2011; Suo et al., 2017) replacing the `0 norm by `1 is differentiable and leads to a sparse solution to Eq. 1. Several authors extended these ideas by considering group sparsity (Chen & Liu, 2012; Klami et al., 2013; Solari et al., 2019). However, these methods are limited to linear transformations and may lead to shrinkage of the canonical vectors due to the `1 regularizer. There has been limited work on extending these models to sparse non-linear CCA. Specifically, there are two kernel-based extensions: two-stage kernel CCA (TSKCCA) by (Yoshida et al., 2017) and SCCA based on Hilbert-Schmidt Independence Criterion (SCCA-HSIC) by (Uurtio et al., 2018). However, these models suffer from the same limitations as KCCA and do not scale to a high dimensional regime.
We present `0-CCA, a simple yet effective method for learning correlated representation based on sparse subsets of the input variables by minimizing an `0 regularized loss. Our `0 regularization relies on a recently proposed Gaussian-based continuous relaxation of Bernoulli random variables, termed gates (Yamada et al., 2020). The gates are applied to the input features to sparsify the canonical vectors. The model and gates parameters are trained jointly via gradient descent to maximize the correlation between the representations ofX and Y while simultaneously selecting only the subsets of most correlated input features. Our contributions are four-folds: (i) By modeling the transformations using two neural networks, we provide a natural solution to the problem of sparse non-linear correlation analysis. (ii) We then propose a novel initialization scheme that helps the model identify small subsets of correlated features in the challenging regime of D > N . (iii) We apply the proposed method to synthetic data and demonstrate that our approach improves the estimation of the canonical vectors compared with existing baselines. (iv) We demonstrate the efficacy of the proposed scheme to several real work applications; specifically, we demonstrate that it leads to more reliable and interpretable representations than other linear and non-linear data fusion schemes.

heading: 2 SPARSE CCA
text: The problem in Eq. 1 has a closed-form solution based on the eigendecomposition of C−1x CxyC −1 y Cyx and C −1 y CyxC −1 x Cxy , where Cx,Cy are within view sample covariance matrices and Cxy,Cyx are cross-view sample covariance matrices. However, if N is smaller than the number of input variables (Dx or Dy), the sample covariance may be rank deficient, and the closed-form solution becomes meaningless. To overcome this limitation, we consider the problem of sparse CCA. Sparse CCA reduces the number of active variables and therefore the covariance matrix becomes full rank. Sparsity also improves model interpretability but may reduce the maximal (training) correlation value obtained in Eq. 1 (Shalev-Shwartz et al., 2010).
Sparse CCA deals with identifying maximally correlated representations based on linear combinations of sparse subsets of the input variables inX and Y . The problem can be formulated as
min a, b
− ρ(aTX, bTY ) + λx‖a‖0 + λy‖b‖0, (2)
where λx and λy are regularization parameters which control the sparsify the input variables. If ‖a‖0 and ‖b‖0 are smaller than N , we can remove the degeneracy inherent to Eq. 1, and identify meaningful correlated representations based on a sparse subset of input variables.

heading: 3 PROBABILISTIC REFORMULATION OF SPARSE CCA
text: The sparse CCA problem formulated in Eq. 2 becomes intractable for large Dx or Dy , moreover, due to the discrete nature of the `0 norm, the problem is not amenable to gradient-based optimization schemes. Fortunately, as demonstrated in sparse regression, probabilistic models such as the spikeand-slab (George & McCulloch, 1993; Kuo & Mallick, 1998; Zhou et al., 2009; Polson & Sun, 2019) provide a compelling alternative. More recently, differentiable probabilistic models such as (Louizos et al., 2017; Yamada et al., 2020) were proposed for sparse supervised learning. Here, we adopt these ideas by rewriting the canonical vectors as α = θx sx and β = θy sy , where denotes the Hadamard product (element wise multiplication), and θx ∈ RDx ,θy ∈ RDy . The vectors sx ∈ {0, 1}Dx , sy ∈ {0, 1}Dy are Bernoulli random vectors with parameters πx = (πx1 , ..., πxDx) and πy = (πy1 , ..., π y Dy ). These Bernoulli variables, act as gates and sparsify the coefficients of the canonical vectors. Now, the problem in Eq. 2 can be reformulated as an expectation minimization, which is parameterized by πx and πy . Specifically, based on the following theorem, we can reformulate the problem in Eq. 2.
Theorem 3.1 The solution of the sparse CCA problem in Eq. 2 is equivalent to the solution of the following probabilistic problem
min πx, πy, θx, θy
E [ − ρ(αTX,βTY ) + λx‖sx‖0 + λy‖sy‖0 ] , (3)
where α = θx sx and β = θy sy and the expectation is taken with respect to the random Bernoulli variables sx and sy (which are parameterized by πx and πy).
Note that the expected values of the `0 norms can be expressed using the Bernoulli parameters as E‖sx‖0 = ∑ πxi , and E‖sy‖0 = ∑ πyi . The proof relies on the fact that the optimal solution to Eq. 2 is a valid solution to Eq. 3, and vice versa. The proof is presented in the Appendix, Section H, and follows the same construction as the proof of Theorem 1 in (Yin et al., 2020).

heading: 3.1 CONTINUOUS RELAXATION FOR SPARSE CCA
text: Due to the discrete nature of sx and sy , differentiating the leading term in Eq. 3 is not straightforward. Although solutions such as REINFORCE (Williams, 1992) or the straight-through estimator (Bengio et al., 2013) enable differentiating through discrete random variables, they still suffer from high variance. Furthermore, they require many Monte Carlo samples for effective training (Tucker et al., 2017). Recently, several authors (Maddison et al., 2016; Jang et al., 2016; Louizos et al., 2017) have demonstrated that using a continuous reparametrization of discrete random variables can reduce the variance of the gradient estimates. Here, following the reparametrization presented in (Yamada et al., 2020), we use Gaussian-based relaxation for the Bernoulli random variables.
Each relaxed Bernoulli variables zi is defined by drawing from a centered Gaussian i ∼ N(0, σ2), then shifting it by µi and truncating its values using the following hard Sigmoid function
zi = max(0,min(1, µi + i)). (4)
Using these relaxed Bernoulli variables, we can define the gated canonical vectors as α = θx zx and β = θy zy . We incorporate these vectors into the objective defined in Eq. 3, in which the `0 regularization terms can be expressed as
E‖zx‖0 = Dx∑ i=1 P(zxi ≥ 0) = Dx∑ i=1 ( 1 2 − 1 2 erf ( − µ x i√ 2σ )) ,
where erf() is the Gaussian error function, and is defined similarly for E‖zy‖0. To learn to model parameters θx,θy and gate parameters µx,µy we first draw realizations for the gates, then we update the parameters by applying gradient descent to minimize
E [ − ρ(αTX,βTY ) + λx‖zx‖0 + λy‖zy‖0 ] .
Post training, we remove the stochasicity from the gates and use all variables such that zxi = max(0,min(1, µxi )) > 0 (defined similarly for z y).

heading: 3.2 GATE INITIALIZATION
text: If the gates are initialized with µi = 0.5, they will approximate “fair” Bernoulli variables. This is a reasonable choice if no prior knowledge about the solution is available; however, we can utilize the closed-form solution of the CCA problem to derive a more suitable parameter initialization for the gate. Specifically, given the empirical covariance matrix Cxy = XY T
(N−1) , we denote the thresholded covariance matrix by C̄xy , with values defined as follows
(C̄xy)ij = { (Cxy)i,j , if |(Cxy)i,j | > δ 0, otherwise,
where δ is the selected threshold value based on the desired sparsity for X and Y , specifically, if we assume that r percent of the values should be zeroed, then δ is set to be the r-th percentile of |(Cxy)|. Then, we compute the leading singular vectors u and v of C̄xy . We further threshold the absolute values of these vectors using the same percentile used for C̄xy . The initial values of the parameters of the gates are then defined by µx = ū+ 0.5, and µy = v̄+ 0.5, where ū and v̄ are the thresholded versions of the absolute value of the singular vectors. Using this procedure, we increase the initial probability for all gates in the support of ū and v̄ based on the singular vectors of C̄xy .
3.3 `0-DEEP CCA
To extend our model to non-linear function estimation, we can formulate the problem of sparse non-linear CCA by modeling the transformations using deep nets as in (Andrew et al., 2013; Wang et al., 2015b). We introduce two random Bernoulli relaxed vectors into the input layers of two neural networks trained in tandem to maximize the total correlation. Denoting the random gating vectors zx and zy for viewX and Y , respectively, our `0-Deep CCA (`0-DCCA) loss is defined by
E [ − ρ̄(f(X̂), g(Ŷ )) + λx‖zx‖0 + λy‖zy‖0 ] , (5)
where f(X̂) = f(zx X|θx) ∈ Rd×N , and g(Ŷ ) = g(zy Y |θy) ∈ Rd×N are modeled as deep networks with model parameters θ = (θx,θy), and gate parameters µ = (µx,µy). The functions f and g embed the data into a d-dimensional space (this formulation could also be used for linear f and g). The functional ρ̄(f(X̂), g(Ŷ )) measures the total correlation between the two d dimensional outputs of the deep nets, this is the sum over d correlation values computed between pairs of coordinates. Exact details on the computation of this term appear in the following subsection. The regularization parameters λx, λy control the sparsity of the input variables. The vectors zx and zy are Bernoulli relaxed vectors, with elements defined based on Eq. 4.
Figure. 1 highlights the proposed architecture. We first pass both observed modalities through the gates. Then, we feed these into view-specific neural sub-nets. Finally, we minimize the shared loss term in Eq. 5 by optimizing the parameters of the gates and the neural sub-nets.

heading: 3.4 ALGORITHM DETAILS
text: Denoting the centered representations for X,Y by Ψx,Ψy ∈ Rd×N (computed using the coupled neural sub-nets), respectively, the empirical covariance matrix between these representations can be expressed as Ĉxy = 1N−1Ψ
x(Ψy)T . Using a similar notations, we express the regularized empirical covariance matrices ofX and Y as Ĉx = 1N−1Ψ x(Ψx)T +γI and Ĉy = 1N−1Ψ y(Ψy)T +γI ,
where the matrix γI (γ > 0) is added to stabilize the invertibility of Ĉx and Ĉy . The total correlation in Eq. 5 (i.e. ρ̄(f(X̂), g(Ŷ ))) can be expressed using the trace of Ĉ −1/2 y ĈyxĈ −1 x ĈxyĈ −1/2 y .
To learn the parameters of the gates µ and of the representations θ we apply full batch gradient decent to the loss in Eq. 5. Specifically, we use Monte Carlo sampling to estimate the left part of Eq. 5. This is repeated for several steps (epochs), using one Monte Carlo sample between each
gradient step as suggested by (Louizos et al., 2017) and (Yamada et al., 2020), and worked well in our experiments. After training, we remove the stochastic part of the gates and use only the variables ix ∈ {1, ..., Dx} and iy ∈ {1, ..., Dy} such that zxix > 0 and z y iy > 0. Finally, we apply linear CCA to the low dimensional extracted embedding to remove this arbitrary scaling and use the normalized generalized eigenvectors (of the closed-form CCA solution) as our canonical vectors. Empirically, we observe that our method also works well with stochastic gradient descent, as long as the batch size is not too small. Alternatively, for small batches we can use a variant of the total correlation loss as was presented in (Wang et al., 2015c) or (Chang et al., 2018). In the Appendix section G, we present a pseudo-code of `0-DCCA, evaluate the complexity of the method, and extend it for a multi-modal setting (more than two views).

heading: 4 RELATED WORK
text: The problem of of `0 sparse CCA was studied in (Wiesel et al., 2008; Cai et al., 2019). However, both rely on a greedy heuristic which iteratively adds non zero elements to the support of a and b. Alternatively, as proposed by (Suo et al., 2017), and `1-regularized problem can be described as
a, b = argmin [ − Cov(aTX, bTY ) + τ1‖a‖1 + τ2‖b‖1 ] ,
subject to ‖aTX‖2 ≤ 1, ‖bTY ‖2 ≤ 1,
where τ1 and τ2 are regularization parameters. A number of variants have been proposed to this problem (Waaijenborg et al., 2008; Parkhomenko et al., 2009; Witten et al., 2009; Hardoon & ShaweTaylor, 2011), but they all suffer from parameter shrinkage and therefore lead to a solution which is less consistent with the correct canonical vectors (see Table 1).

heading: 5 EXPERIMENTAL RESULTS
text: We validate the effectiveness of the proposed approach on a wide range of tasks. First, using synthetic data, we demonstrate that `0-CCA correctly identifies the canonical vectors in a challenging regime of N < Dx, Dy . Next, using a coupled video dataset, we demonstrate that `0-CCA can identify the common information from high dimensional data, and embed it into correlated, lowdimensional representations. Then, we use noisy images from MNIST and multi-channel seismic data to demonstrate that `0-DCCA finds meaningful representations of the data even in a noisy regime. Finally, we use `0-DCCA to improve cancer sub-type classification using high dimensional genetic measurements. In all experiments validation sets are used for tuning the hyperparameters of all baselines by maximizing the total correlation on the validation set. We refer the reader to the Appendix for a complete description of the baselines, training procedure, and parameters choice for all methods.
sentations, i.e. ρ̂ = φ̂ T XTY η̂, where φ̂ and η̂ are the estimated canonical vectors. Dashed lines indicate the correct number of active coefficients (10) and true correlation ρ (0.9). Note that for small values of λ = λx = λy , the model selects many variables and attains a higher correlation value; in this case, `0-CCA suffers from overfitting. Right: True canonical vector φ along with the estimated vectors using `0-CCA (φ̂) and CCA (â). Due to the small sample size, CCA overfits and fails to identify the correct canonical vectors.

heading: 5.1 SYNTHETIC EXAMPLE
text: To generate samples fromX,Y ∈ RD×N , we follow the procedure described in (Suo et al., 2017), considering data sampling from the following distribution ( X Y ) ∼ N( ( 0 0 ) , ( Σx Σxy Σyx Σy ) ), where Σxy = ρ0Σx(φηT )Σy . We study three cases for the covariance matrices Σ = Σx = Σy . Model I. Identity: Σ = ID. Model II. Toeplitz: Σi,j = ρ |i−j| 0 , i, j = 1, ..., D. Model III. Sparse inverse: Σi,j = Σ̄i,j√
Σ̄i,iΣ̄j,j , where Σ̄ = Γ−1, Γi,j = 1i=j + 0.51i=j + 0.41i=j .
For all three cases, the vectors φ,η ∈ RD, are sparse with 5 nonzero elements and ρ0 = 0.9. The indices of the active elements are chosen randomly with values equal to 1/ √ 5. In this setting, the canonical vectors a and b maximizing the correlation objective in Eq. 1 are φ and η, respectively (see Proposition 1 in (Suo et al., 2017)).
Using Model I, we first generate N = 400 samples, with D = 800, and estimate the canonical vectors based on CCA and `0-CCA. In Fig. 2, we present a regularization path of the proposed scheme. Specifically, we apply `0-CCA to the data described above using various values of λ = λx = λy . We present the `0 of active gates (by expectation) along with the empirical correlation between the extracted representations, defined by ρ̂ = ρ(φ̂ T X, η̂TY ). As evident from the left panel, a wide range of λ values leads to the correct number of active coefficients (10) and the correct correlation value (ρ0 = 0.9). Next, in the right panel we present the values of φ, the `0-CCA estimate (using λ = 30) of the canonical vector φ̂, and the CCA spectral based estimate of the canonical vector â. Due to the low number of samples, the solution by CCA is wrong and not sparse, while the `0-CCA solution correctly identifies the support of φ.
Next, we evaluate the estimation error of φ using eφ = 2(1− |φ T φ̂|), and eη is defined similarly. In Table 1 we present the estimation errors of φ and ρ (averaged over 100 simulations) for Models I, II and III (identity, Toeplitz and sparse inverse covariance matrices). As baselines, we compare
the performance to 5 leading sparse CCA models. As evident from these experiments, `0-CCA significantly outperforms all baselines in its ability to learn the correct canonical vectors. Here, we restrict our evaluation to five baselines since code was not available for other sparse CCA models. We have implemented (Wiesel et al., 2008) but observed that it does not converge to the correct canonical vectors (the estimation errors eφ and eρ were 2). In the Appendix, section F we provide a runtime evaluation of the method for different values of N and D.

heading: 5.2 MULTI VIEW OF SPINNING PUPPETS
text: As an illustrative example, we use a dataset collected by (Lederman & Talmon, 2018). The authors have generated two videos capturing rotations of 3 desk puppets. One camera captures two puppets, while the other captures another two, where one puppet is shared across cameras. A snapshot from both cameras appears in Fig. 3. All puppets are placed on spinning devices that rotate the dolls at different frequencies. In both videos, there is a shared underlying parameter, namely the rotation of the common bulldog. We use 400 images from each camera, where each image has 240 × 320 = 76800 pixels (using a grayscaled version of the colored image); therefore, N D, and direct application of CCA would fail. We apply the proposed scheme using λy = λx = 50, a linear activation and embedding dimension d = 2. `0-CCA converges to embedding with a total correlation of 1.99 using 372 and 403 pixels from views X and Y . The active gates are presented in the right panels of Fig. 3. In this example, the gates correctly attenuate the modality-specific information and preserve subsets of pixels capturing common information, allowing for embedding the gated videos into a shared correlated representation.
In Fig. 4, we present the coupled two-dimensional embedding of both videos. Both embeddings are correlated with each other and reveal the correct angular orientation of the Bulldog (which is
the common latent parameter in this experiment). This indicates the ability of `0-CCA to identify correlated latent variables in the regime of N Dx, Dy .

heading: 5.3 NOISY MNIST
text: We use two noisy variants of MNIST (LeCun et al., 2010) as our coupled views. The first view is created by adding noise drawn uniformly from [0, 1] to all pixels. The second view is created by placing a random patch from a natural image in the background of the handwritten digits. Random samples from these modalities are presented in Fig. 5. Each view consists of 62, 000 samples, of which we use 40, 000 for training, 12, 000 for testing and 10, 000 are used as a validation set.
We train `0-DCCA to embed the data into a correlated 10 dimensional space while selecting subsets of input pixels. Our model selects 277, and 258 pixels from both modalities respectively (see bottom right corner of Fig. 5). Next, we evaluate the quality of learned embedding by applying k-means to the stacked embedding of both views. We run k-means (with k = 10) using 20 random initializations and record the run with the smallest sum of square distances from centroids. Given the cluster assignment, k-means clustering accuracy (KM) and mutual information (MI) are measured using the true labels. Additionally, we train a Linear-SVM (SVM) model on our train and validation datasets. SVM classification accuracy is measured on the remaining test set. The embedding provided by `0-DCCA leads to higher classification and clustering results compared with several linear and nonlinear modality fusion models appear in Table 2. In the Appendix, we provide the implementation details and present experimental results demonstrating the accuracy vs. sparsity trade-off of `0- DCCA.

heading: 5.4 SEISMIC EVENT CLASSIFICATION
text: Next, we evaluate the method using a dataset of seismic events studied by (Lindenbaum et al., 2018). Here, we focus on 537 explosions which are categorized into 3 quarries. Each event is recorded using two directional channels facing east (E) and north (N); these comprise the coupled views for the correlation analysis. Following the analysis by (Lindenbaum et al., 2018), the input features are sonogram representations of the seismic signal. Sonograms are time-frequency representations with bins equally tempered on a logarithmic scale. Each sonogram z ∈ R1157 with 89 time bins and 13 frequency bins. An example of sonograms from both channels appears in the top row of Fig. 6.
We create the noisy seismic data by adding sonograms computed based on vehicle noise from 1. Examples of noisy sonograms appear in the middle row of Fig. 6. We hold out 20% of the data as a validation set, and train `0-DCCA to embed the data in 3 dimensions. In Table 2 we present the MI, k-means and SVM accuracies computed based on `0-DCCA embedding. Furthermore, we compare the performance with several other baselines. Here, the proposed scheme improves performance in all 3 metrics while identifying a subset of 17 and 16 features from channel E and N, respectively. The active gates are presented in the bottom row of Fig. 6. Our results indicate that even in the
1https://bigsoundbank.com/search?q=car
presence of strong noise, `0-DCCA correctly activates the gates in frequency bins that coincide with the energy stamps of the primary and secondary waves (P and S in the top left of Fig. 6).

heading: 5.5 CANCER SUB-TYPE CLASSIFICATION
text: Accurate classification of cancer sub-type is vital for extending the life span of patients by personalized treatments (Zhu et al., 2020). This task is challenging since the number of measured genes is typically much larger than the number of observations. Here, we use multi-modal observations from the METABRIC data (Curtis et al., 2012) and attempt to find correlated representations to improve cancer sub-type classification. The data consists of 1, 112 breast cancer patients which are annotated by 10 subtypes based on InClust (Dawson et al., 2013). We observe two modalities, namely the RNA gene expression data, and Copy Number Alteration (CNA) data. The dimensions of these modalities are 15, 709 and 47, 127, respectively. We compute the `0-DCCA 10 dimensional embedding (and all baseline embeddings) and demonstrate using k-means and SVM that the representation identified using `0-DCCA can lead to more accurate cancer sub-type classification (see Table 2).

heading: 6 CONCLUSION
text: This paper presents a method for learning sparse non-linear transformations that maximize the canonical correlations between two modalities. Our approach is realized by gating the input layers of two neural networks, which are trained to optimize their output’s total correlations. Input variables are gated using a regularization term which encourages sparsity. We further propose a novel scheme to initialize the gates based on a thresholded cross-covariance matrix. Our method can learn informative correlated representations even when the number of variables far exceeds the number of samples. Finally, we demonstrate that the proposed scheme outperforms existing algorithms for linear and non-linear canonical correlation analysis. We describe our future directions in the Appendix section I.

heading: A APPENDIX

heading: B ADDITIONAL EXPERIMENTAL DETAILS
text: In the following sections we provide additional experimental details required for reproduction of the experiments provided in the main text. All the experiments were conducted using Intel(R) Xeon(R) CPU E5-2620 v3 @2.4Ghz x2 (12 cores total).
B.1 SELECTING THE EMBEDDING DIMENSION
Choosing the embedding dimension might require application-specific considerations. For example, for clustering purposes, if the number of clusters is known (K), we might select d = K to have a sufficient number of dimensions to separate the clusters. However, for other applications, we can use a validation set to tune the embedding dimension d; precisely, by increasing d until the total correlation in the validation set stops increasing. Using this procedure we can estimate the most significant number of correlated embedding coordinates.
B.2 SYNTHETIC EXAMPLE
For the linear model we use a learning rate of 0.005 with 10, 000 epochs. The values of λx and λy are both set to 30. These values were obtained using a cross validation procedure. The standard deviation σ of the injected noise was set to 0.5 by (Yamada et al., 2020). They have selected this value as it maximized the gradient of the regularization term at initialization. Empirically, we observe that for `0-CCA smaller values of σ translated to improved convergence. Specifically, we used σ = 0.25, which worked well in our experiments.
For each data model, we run the method 100 times with different realizations ofX and Y . In Table 1, we compare the average results of `0-CCA to PMA (Witten et al., 2009), IP-SCCA (Mai & Zhang, 2019), SCCA-I (Hardoon & Shawe-Taylor, 2011),SCCA-II (Gao et al., 2017), and mod-SCCA (Suo et al., 2017). Results of SCCA-II, PMA, and mod-SCCA were simulated by (Suo et al., 2017).
To demonstrate the tuning procedure of λx and λy , we generate data based on Model-I (described in Section 5.1) and evaluate the total correlation on training and validation samples. The results are presented in Fig. 7. As demonstrated in this plot, the correlation values on unseen samples can indicate the optimal values for the regularization parameters. Namely, the model leads to lower correlation values when selecting nuisance variables or removing some of the informative variables.
sentations, i.e. ρ̂ = φ̂ T XTY η̂, where φ̂ and η̂ are the estimated canonical vectors. Dashed lines indicate the correct number of active coefficients (10) and true correlation ρ (0.9). In this example, we present the correlation values training samples (blue circles), and on unseen samples (blue asterisks). Red vertical lines indicate the range of λ’s that maximize the correlation on the validation set. Within this range of values the model selects the “correct” set of coefficients.
B.3 NOISY MNIST
In this subsection we provide additional details regarding the noisy MNIST experiment. In Fig. 8, we present the performance as a function of the number of active gates (pixels) controlled by λx = λy . The Mutual Information (MI) score, k-means, and linear Support Vector Machine (LSVM) accuracies were computed based on `0-DCCA embedding with learning rate of 0.01. Furthermore, the number of epochs (∼ 4000) was tuned by early stopping using a validation set of size 10000. To learn 10 dimensional correlated embedding, we use the same architecture as suggested by (Wang et al., 2015a) consisting of three hidden layers with 1000 neurons each. The number of dimensions in the embedding was selected based on the number of classes in MNIST. This architecture is used for `0-DCCA, DCCA (Andrew et al., 2013) and DCCAE (Wang et al., 2015b). Note that for `0-DCCA using small values of the regularization parameters λx and λy , increases the number of selected features and degrades the performance. This is duo to the fact that as more features are selected more noise is introduced into the extracted representation (of size 10). It is interesting to note that the k-means was more robust to the introduced noise than the LSVM.
The regularization parameters λx and λy balance between the correlation loss and the amount of sparsification performed by the gates. These hyper parameters are tuned using the validation set by maximizing the total correlation value. We compare `0-DCCA to CCA (Chaudhuri et al., 2009)2, mod-SCCA (Suo et al., 2017), SCCA-HSIC (Uurtio et al., 2018), KCCA (Bach & Jordan, 2002)3,
2https://scikit-learn.org/stable/modules/generated/sklearn.cross decomposition.CCA.html 3https://gist.github.com/yuyay/16ce4914683da30f87d0
NCCA (Michaeli et al., 2016) 4, multiview-ICA (Richard et al., 2020), DCCA (Andrew et al., 2013) 5 and DCCAE (Wang et al., 2015b). For all methods we use an embedding with dimension 10, and evaluate performance with k-means using 20 random initilizations, and using LSVM by performing training on the training samples and testing on the remaining samples (split defined in the main text). The hyperparameters of all methods are tuned to maximize the total correlation on a validation set. The solution obtained by mod-SCCA is based on ∼ 320 variables. In this experiment we tried to train SCCA-HSIC (Uurtio et al., 2018) 6 for over two days, but it did not converge. Furthermore, we believe that the poor performance of the kernel methods stems from the high level of noise in the input images. We note that grad-KCCA (Uurtio et al., 2019) did not converge in this experiment, therefore we did not report its performance on MNIST.
B.4 SEISMIC EVENT CLASSIFICATION
Using the seismic data, we compare the performance of `0-DCCA with a linear and non-linear activation. In this example, we use a learning rate of 0.01 with 2000 epochs. The number of neurons for the five hidden layer are: 300, 200, 100, 50, and 40 respectively, with a tanh activation after each layer. The number of dimensions in the embedding (d = 3) was selected based on the number of classes in the data. Parameters are optimized manually to maximize the correlation on a validation set. In Fig. 9, we present SVM accuracy for different levels of sparsity. The presented number of features is the average over both modalities, and SVM performance is evaluated using 5-folds cross validation. We compare `0-DCCA to CCA (Chaudhuri et al., 2009), mod-SCCA (Suo et al., 2017), SCCA-HSIC (Uurtio et al., 2018), KCCA (Bach & Jordan, 2002), NCCA (Michaeli et al., 2016), multiview-ICA (Richard et al., 2020), DCCA (Andrew et al., 2013), and DCCAE (Wang et al., 2015b). For all methods we use an embedding with dimension d = 3, and evaluate performance with k-means using 20 random initilizations, and using linear SVM by performing a 5-folds cross validation. For the kernel methods we evaluated performance by constructing a kernel using k = 5, 10, ..., 50, nearest neighbors and selected the value which maximized performance in terms of total correlation. In this example, the optimal solution obtained mod-SCCA Suo et al. (2017) is based on ∼ 200 input variables, while the solution obtained SCCA-HSIC Uurtio et al. (2018) is based on∼ 77 variables. This means that our method uses much fewer input variables and leads to substantially higher clustering and classification accuracies (see Table 3).
4https://tomer.net.technion.ac.il/files/2017/08/NCCAcode v3.zip 5https://github.com/adrianna1211/DeepCCA tensorflow 6https://github.com/aalto-ics-kepaco/scca-hsic

heading: C LOW DIMENSIONAL SYNTHETIC EXAMPLE
text: Following an interesting suggestion raised during the review process, we provide here an evaluation of the method for the regime of D N . We use a noisy version of the synthetic data model presented in Section 5.1. Specifically, we generate samples from X,Y ∈ RD×N , using the following
distribution ( X Y ) ∼ N (( 0 0 ) , ( Σx Σxy Σyx Σy )) , where Σxy = ρ0Σx(φηT )Σy . With Σ = ID.
The we use additive noise to create the two modalities, namely X̃ = X+ X , Ỹ = Y + Y . Where X and Y are drawn from N(0, IDσN ). We use D = 20 and N = 4000 to evaluate if our method is beneficial in the regime where there are more samples than variables. In table 3 we present the average errors for different noise levels (values of σN ). As evident from this table, the proposed approach can provide a more accurate estimate of the sparse canonical vectors in the regime of D N .

heading: D CANCER SUB-TYPE CLASSIFICATION
text: In the main text, we demonstrated that the embedding extracted with `0-DCCA leads to more accurate cancer sub-type classification. Here, we provide additional details on the experiment described in Section 5.5. In this example, we use a learning rate of 0.5 with 2000 epochs. The number of neurons for the 3 hidden layers are: 500, 300, 100, with a tanh activation after each layer. The number of dimensions in the embedding (d = 10) was selected based on the number of classes in the data. Parameters are optimized manually to maximize the correlation on a validation set. We compare `0-DCCA to CCA (Chaudhuri et al., 2009), mod-SCCA (Suo et al., 2017), SCCA-HSIC (Uurtio et al., 2018), KCCA (Bach & Jordan, 2002), NCCA (Michaeli et al., 2016), multiview-ICA (Richard et al., 2020), DCCA (Andrew et al., 2013), and DCCAE (Wang et al., 2015b). We use the same SVM and k-means scheme as described in Section B.4. The solution obtained by mod-SCCA is based on ∼ 940 variables.
Algorithm 1 `0-DCCA Input: Coupled data, {X,Y }, regularization parameters λx, λy , number of epochs T . Initialize the gate parameters: µxi = 0.5 for i = 1, . . . , D x, and µyi = 0.5 for i = 1, . . . , D y .
for t = 1 to T do Sample a stochastic gate (STG) vectors zx, zy defined based on equation 4. Apply the STG to the data X̂ = zx X and Ŷ = zy Y . Compute the loss L = E [ − ρ̄(f(X̂), g(Ŷ )) +λx‖zx‖0 +λy‖zy‖0 ] , where the calculation of
the total correlation ρ̄ is based on Section 3.4. Where f(X̂) = f(zx X|θx) ∈ Rd×N , and g(Ŷ ) = g(zy Y |θy) ∈ Rd×N . Update θx = θx − γ∇θxL, θ
y = θy − γ∇θyL, πx = πx − γ∇πxL, πy = πy − γ∇πyL.
end for Optional: apply linear CCA on f(X̂), and g(Ŷ )) to remove arbitrary scalings.
To provide additional intuition on why sparsity is useful in biomedical applications, we further analyze the results presented in Section 5.5 (the METABRIC dataset). In the experiment provided in this section, our method selects 169 genes from the RNA data and extracts an embedding that leads to improved identification (clustering and classification) of cancer subtypes. To understand if the variables selected by our model are interpretable, we analyze genes selected by the model. We have observed that out of the 169 selected genes, six are part of the known risk factors for breast cancer, namely the well-known PAM50 list. The PAM50 contains 50 genes known to be critical variables for classifying breast cancer subtypes Bernhardt et al. (2016). Our model has identified six genes (in an unsupervised fashion) out of the 42 possible PAM50 genes that were initially available in the METABRIC data. To demonstrate that this is a statistically significant finding, we compute the probability of selecting 169 genes out of a total of 15019 genes while having an overlap of 6 genes with another list of size 42. This probability is 5.39E-6 (the p-value), which indicates that this overlap is unlikely to have happened by chance. To conclude, our method was able to identify known markers for breast cancer subtypes classification without using any labels. Notably, one of these genes is ‘ERBB2’, which corresponds to the HER2/neu protein. ‘ERBB2’ is known to have a large copy number in patients classified to the HER2 cancer subtype. Thus the method found a strong association between the RNA and CNA of this gene.

heading: E PSEUDOCODE
text: In Algorithm 1 we provide a pseudocode description of the proposed approach.

heading: F RUN TIME ANALYSIS
text: Our method requires computing covariance matrices for kx and ky variables with N samples at each iteration. This scales like O(Nk2x + Nk2y). To demonstrate the computational efficiency of our method, we run `0-CCA for different values of N and D and evaluate the empirical computational complexity of the method. In Fig. 10 we present the average runtime over 100 runs, the data is generates following Model I from Section 5.1.
G GENERALIZED `0-BASED SPARSE CCA
The proposed `0-CCA algorithm was designed and demonstrated on coupled datasets. However, in many applications, one may have access to multiple (> 2) datasets. There are many real-life cases where more than two modalities are available such as hyper-spectral imaging, medical imaging and more. The problem of CCA was generalized decades ago to include multiple data sets by (Horst, 1961; McKeon, 1966; McDonald, 1968) and (Carroll, 1968) to name some. These approaches were later summarized in (Kettering, 1971). Since then, many extensions have been proposed such as (Van de Geer, 1984). Some of the modern extensions of multi-view CCA comprise its regularised (Tenenhaus & Tenenhaus, 2011), kernelised (Tenenhaus et al., 2015) sparse (Tenen-
haus et al., 2014) and deep (Benton et al., 2019) variants. In this section, we extend the proposed `0-DCCA to multi-view datasets.
Let K be the number of views and Xk the corresponding k-th view 1 ≤ k ≤ K. The Deep Generalized CCA that was proposed in (Benton et al., 2019) aims to solve
min G∈Rn×d,U k∈Rd×d K∑ k=1 ‖G−UTk f(X k)‖F s.t. GGT = I, (6)
where G is the desired shared representation, and Uk is a linear transformation from the k-th deep network output, f(Xk), to the shared representation space.
While the proposed method in Section 3.3 describes a solution for coupled views, this solution can be naturally extended to include multiple views by introducing a shared common space into the optimization problem in Eq 6. Denoting the random gating vectors zk for viewXk. The generalized `0-DCCA loss is defined by
min G∈Rn×d,U k∈Rd×d,zk∈Rdk K∑ k=1 ‖G−UTk f(X̂ k)‖F + λk‖zk‖0 s.t. GGT = I, (7)
where f(X̂k) = f(zk Xk|θx) ∈ Rd×N and λk is the regularization parameter that control the sparsity of k-th selected gate vector zk. The optimization problem in Eq 7 can be solved using standard optimizers such as gradient descent. The initialization of G and Uk plus the analysis of the above suggested generalized gated canonical correlation are left for further research.

heading: H PROOF OF THEOREM 1
text: We first note that the equivalence of the solutions means that if we denote the solution to Eq. 2 by (â, b̂), then there is a corresponding set of values πx,πy,θx,θy , such that (α̂, β̂) = (â, b̂), where α = θx sx and β = θy sy . The proof follows the same construction as the proof of Theorem 1 in (Yin et al., 2020). To prove the Theorem, we will show that the optimal solution to the deterministic sparse CCA problem (in Eq. 2) is a valid solution to the probabilistic sparse CCA problem (in Eq. 3) and vice versa.
First, we denote by (â, b̂) the optimal solution of the problem defined in Eq. 2, then by setting α̂ = â ŝx and β̂ = b̂ ŝy , where π̂xi = 1 if âi 6= 0 and πxi = 0 otherwise, we get a feasible solution to Eq. 3 which leads to the same objective value.
Now, we want to show that the optimal solution to Eq. 3 is also a feasible solution to Eq. 2. Denoting the optimal solution to Eq. 3 by α̂ and β̂, with α̂ = θ̂ x ŝx and β̂ = θ̂ y ŝy ,
if p(sx|πx) or p(sy|πy) are point mass densities (i.e., πx,πy ∈ {0, 1}D), then the solution is effectively deterministic and â = α̂, b̂ = β̂ would be a valid solution to problem Eq. 2. If p(sx|πx) and p(sy|πy) are not point mass densities, we will now show by contradiction that all the points in the support of p(sx|πx) (or p(sy|πy)) would lead to the same objective values. Assume without loss of generality that {sx1 , ...., sxL} ∈ supp[p(sx|π̂
x)], if not all values in the support lead to the same objective value, this means that there exist sxl , s x m such that L(s x l ) < L(s x m), where
L(sx) = [ − ρ((θ̂ x sx)TX, (θ̂ y ŝy)TY ) + λx‖sx‖0 + λy‖ŝy‖0 ] ,
which is the objective from Eq. 3 with optimal values of θ̂ x , θ̂ y , ŝy . Now if we set π̄x = sxl we would obtain Esx∼p(sx|π̄x)L(sx) < Esx∼p(sx|π̂x)L(sx) which contradicts the assumption that π̂ x is optimal. This means that all points in supp[p(sx|π̂x)] lead to the same objective value, therefore, they are also feasible solutions to Eq. 2. Now, showing that the solution to Eq. 3 is a feasible solution to Eq. 2 and vice versa completes our proof.

heading: I LIMITATIONS AND FUTURE DIRECTIONS
text: The proposed method provides an effective solution to the problems of sparse linear and nonlinear CCA. One advantage of the suggested method compared with existing sparse CCA solutions, is that it can embed data into a d > 1 dimensional space and learn the sparsity pattern simultaneously. The proposed `0-CCA problem albeit not being convex leads to an empirically stable solution across different settings. Nonetheless, the method has some limitations, specifically, tuning λ requires a cross validation procedure, which could be costly in high dimensional regime. Another caveat of the existing approach is that it lacks guarantees when trained on small batches. In the future, we aim to extend the method to enable compatibility with small batch training. Furthermore, we are currently working on evaluating how the estimation error scales with respect to the dimension, sparsity, for different structures of the covariance matrix. Furthermore, we are working on bounding the estimation error of our method. We observe that using a Gaussian based data model (similar to the one used in section 5.1), it is possible to bound the estimation error (in probability) byO(D/N+ck/N) where D is the number of variables, N is the number of measurements, k is the true sparsity of the conical vectors, and c is an independent constant. Although this bound provides some guarantee, it is ubiquitous in the context of the regime studied in our manuscript (i.e., N D). Therefore, we are currently working on improving this bound. Empirically, we observe that the method works well in regimes that existing `1 and `0 based approaches fail.

Review summaries:
The paper proposes a sparse canonical correlation analysis method based on l_0 norm. A continuous relaxation scheme is adopted for solving sparse CCA. The proposed model is then extended to nonlinear function estimation and combined with deep neural networks. Experimental results on synthetic and real-world datasets demonstrate the effectiveness of the proposed methods.

The paper proposes a new method for l0-CCA using stochastic gating which allows for an efficient algorithm and also permits a deep-version of the l0-CCA. Results are shown on synthetic and real-world datasets that highlight the superior performance of the proposed approach.   Main Contributions:   1). The paper proposes a new approach for l0-CCA that uses a continuous relaxation of a bernoulli random variable via stochastic gating.   2). The paper extends the proposed method to l0-deep-CCA which allows non-linear interactions between the two CCA views.  3). Results are shown on synthetic and real-world datasets which show the superior performance of the proposed l0-CCA methods.