# ℓ0-SPARSE CANONICAL CORRELATION ANALYSIS


**Ofir Lindenbaum**
Faculty of Engineering
Bar Ilan University
Ramat Gan, Israel
ofirlin@gmail.com


**Moshe Salhov[∗]** **& Amir Averbuch**
School of Computer Science
Tel Aviv University
Tel Aviv, Israel

ABSTRACT


**Yuval Kluger**
School of Medicine
Yale University
New Haven, CT, USA
yuval.kluger@yale.edu


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.

1 INTRODUCTION

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 modalities X ∈ R[D][x][×][N] and Y ∈ R[D][y][×][N], which are centered and have N samples with D[x] and
_D[y]_ features, respectively, CCA seeks canonical vectors a ∈ R[D][x] _, and b ∈_ R[D][y] _, such that, u =_
**_a[T]_** **_X, and v = b[T]_** **_Y, will maximize the sample correlations between the canonical variates, i.e._**

**_a[T]_** **_XY_** **_[T]_** **_b_**
max _ρ(a[T]_ **_X, b[T]_** **_Y ) =_** _._ (1)
**_a, b=0̸_** _∥a[T]_ **_X∥2∥b[T]_** **_Y ∥2_**

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 > D[x], D[y]. However, if we have

_∗M.S. is co affiliated with Playtika Israel_


-----

Figure 1: Illustration of ℓ0-DCCA. Data from two observed views propagate through stochastic
gates (defined in Eq. 4). The gates output is fed into two neural sub-nets that have a shared loss (see
Eq. 5). We compute this shared loss based on the neural sub-nets outputs (with dimension d = 3 in
this example). Our shared loss combines a total correlation term with a differentiable regularization
term which induces sparsity in the input variables (by sparsifying the gates).

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 < D[x], D[y] (once the number of
selected variables is smaller than N 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 sub-_
optimal 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 of X 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.

2 SPARSE CCA

The problem in Eq. 1 has a closed-form solution based on the eigendecomposition of
**_C_** _[−]x_ [1][C] _[xy][C]_ _y[−][1][C]_ _[yx]_ [and][ C] _y[−][1][C]_ _[yx][C]_ _x[−][1][C]_ _[xy][, where][ C]_ _[x][,][ C]_ _[y]_ [are within view sample covariance ma-]
trices and C _xy, C_ _yx are cross-view sample covariance matrices. However, if N is smaller than_
the number of input variables (D[x] or D[y]), 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 in X and Y . The problem can be formulated as

mina, b _−_ _ρ(a[T]_ **_X, b[T]_** **_Y ) + λ[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.

3 PROBABILISTIC REFORMULATION OF SPARSE CCA

The sparse CCA problem formulated in Eq. 2 becomes intractable for large D[x] or D[y], 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; Jana et al., 2021) were proposed for sparse supervised
learning. Here, we adopt these ideas by rewriting the canonical vectors as α = θ[x] _⊙_ **_s[x]_** and
**_β = θ[y]_** _⊙_ **_s[y], where ⊙_** denotes the Hadamard product (element wise multiplication), and θ[x] _∈_
R[D][x] _, θ[y]_ _∈_ R[D][y] . The vectors s[x] _∈{0, 1}[D][x]_, s[y] _∈{0, 1}[D][y]_ are Bernoulli random vectors with
parameters π[x] = (π1[x][, ..., π]D[x] _[x]_ [)][ and][ π][y][ = (][π]1[y][, ..., π]D[y] _[y]_ [)][. 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_

**_π[x], πmin[y], θ[x], θ[y][E]_** _−_ _ρ(α[T]_ **_X, β[T]_** **_Y ) + λ[x]∥s[x]∥0 + λ[y]∥s[y]∥0_** _,_ (3)
 

_where α = θ[x]_ _⊙_ **_s[x]_** _and β = θ[y]_ _⊙_ **_s[y]_** _and the expectation is taken with respect to the random_
_Bernoulli variables s[x]_ _and s[y]_ _(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∥s[x]∥0 = _πi[x][, and][ E][∥][s][y][∥][0][ =][ P][ π]i[y][. 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).

[P]

3.1 CONTINUOUS RELAXATION FOR SPARSE CCA

Due to the discrete nature of s[x] and s[y], 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; Lindenbaum et al., 2021), we use Gaussian-based relaxation for the Bernoulli
random variables.

then shifting it byEach relaxed Bernoulli variables z µi and truncating its values using the following hard Sigmoid functioni is defined by drawing from a centered Gaussian ϵi ∼ _N_ (0, σ[2]),

zi = max(0, min(1, µi + ϵi)). (4)

Using these relaxed Bernoulli variables, we can define the gated canonical vectors as α = θ[x] _⊙_ **z[x]**
and β = θ[y] **z[y]. We incorporate these vectors into the objective defined in Eq. 3, in which the ℓ0**
_⊙_
regularization terms can be expressed as


21 2 [erf] _−_ _√[µ]2i[x]σ_

 _[−]_ [1] 





_D[x]_

P(z[x]i
_i=1_ _[≥]_ [0) =]

X


_D[x]_

_i=1_

X


E∥z[x]∥0 =


-----

where erf() is the Gaussian error function, and is defined similarly for E∥z[y]∥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 _−_ _ρ(α[T]_ **_X, β[T]_** **_Y ) + λ[x]∥z[x]∥0 + λ[y]∥z[y]∥0_** _._

Post training, we remove the stochasicity from the gates and use all variables such that z  _[x]i_ =
max(0, min(1, µ[x]i [))][ >][ 0][ (defined similarly for][ z][y][).]

3.2 GATE INITIALIZATION

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 C _xy =_ **_[XY](N_** **_[ T]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
(C _xy)_ . 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
**z[x]** and z[y] for view X and Y, respectively, our ℓ0-Deep CCA (ℓ0-DCCA) loss is defined by

E _−_ _ρ¯(f_ **( X[ˆ]** ), g( Y[ˆ] )) + λ[x]∥z[x]∥0 + λ[y]∥z[y]∥0 _,_ (5)
 

where f **( X[ˆ]** ) = f **(z[x]** **_⊙_** **_X|θ[x]) ∈_** R[d][×][N], and g( Y[ˆ] ) = g(z[y] **_⊙_** **_Y |θ[y]) ∈_** R[d][×][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 z[x] and z[y] 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.

3.4 ALGORITHM DETAILS

Denoting the centered representations for X, Y by Ψ[x], Ψ[y] _∈_ R[d][×][N] (computed using the coupled
neural sub-nets), respectively, the empirical covariance matrix between these representations can be
expressed as **_C_** _xy =_ _N1_ 1 **[Ψ][x][(][Ψ][y][)][T][ . Using a similar notations, we express the regularized empiri-]**

cal covariance matrices of− **_X and Y as_** **_C_** _x =_ _N1_ 1 **[Ψ][x][(][Ψ][x][)][T][ +][γ][I][ and][ b]C** _y =_ _N1_ 1 **[Ψ][y][(][Ψ][y][)][T][ +][γ][I][,]**

_−_ _−_

[b]

[b]


-----

where the matrix γI (γ > 0) is added to stabilize the invertibility of **_C_** _x and_ **_C_** _y. The total correla-_

_−1/2_ _−1_ _−1/2_
tion in Eq. 5 (i.e. ¯ρ(f **( X[ˆ]** ), g( Y[ˆ] ))) can be expressed using the trace of **_C_** _y_ **_C_** _yxC_ _x_ **_C_** _xyC_ _y_ .

[b] [b]

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[b] b [b] [b] [b]
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
_i[x]_ _∈{1, ..., D[x]} and i[y]_ _∈{1, ..., D[y]} such that z[x]i[x][ >][ 0][ and z][y]i[y][ >][ 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).

4 RELATED WORK

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), an ℓ1-regularized problem can be described as

**_a, b = argmin_** _−_ Cov(a[T] **_X, b[T]_** **_Y ) + τ1∥a∥1 + τ2∥b∥1_** _,_

subject to  **_a[T]_** **_X_** 2 1, **_b[T]_** **_Y_** 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).


5 EXPERIMENTAL RESULTS

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 < D[x], D[y]. 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.

5.1 SYNTHETIC EXAMPLE

To generate samples from X, Y ∈ R[D][×][N], we follow the procedure described in (Suo et al., 2017),

considering data sampling from the following distribution **_X_** _N_ ( **0** _,_ **Σx** **Σxy** ),
**_Y_** _∼_ **0** **Σyx** **Σy**
     

where Σxy = ρ0Σx(φη[T] )Σy. We study three cases for the covariance matrices Σ = Σx = Σy.
**Model I. Identity: Σ = I** _D. Model II. Toeplitz: Σi,j = ρ[|]0[i][−][j][|], i, j = 1, ..., D._
**Model III. Sparse inverse: Σi,j =** _√Σ¯Σ¯i,ii,jΣ[¯]_ _j,j_ [, where][ ¯]Σ = Γ[−][1], Γi,j = 1i=j + 0.51i=j + 0.41i=j.

For all three cases, the vectors φ, η ∈ R[D], 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


-----

Figure 2: Left: Regularization path of ℓ0-CCA on data generated from the linear model I (described
in Section 5.1). Values on the left y-axis (green) represent the sum of active gates (by expectation).
Values on the right y-axis (blue) represent the empirical correlation between the estimated representations, i.e. ˆρ = **_φ[ˆ]T X_** _T Y ˆη, 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.

|canonical vectors a and b maximizing the co (see Proposition 1 in (Suo et al., 2017)).|orrelation objective in Eq. 1 are φ and η, respectively|Col3|Col4|
|---|---|---|---|
||Model I- Identity covariance|||
|(N, Dx, Dy)|(400,800,800)|(500,600,600)|(700,1200,1200)|
|PMA (Witten et al., 2009) IP-SCCA (Mai & Zhang, 2019) SCCA-I (Hardoon & Shawe-Taylor, 2011) SCCA-II (Gao et al., 2017) mod-SCCA (Suo et al., 2017) ℓ -CCA 0|(1.170,1.170) (1.658,1.647) (1.602,1.140) (0.060,0.066) (0.056,0.062) (0.003,0.009)|(0.850,0.850) (1.051,1.051) (1.143,0.282) (0.053,0.057) (0.05,0.056) (0.002,0.002)|(1.090,1.090) (1.544,1.542) (1.160,0.181) (0.045,0.043) (0.045,0.043) (0.001,0.002)|
||Model II- Toeplitz covariance|||
|PMA (Witten et al., 2009) IP-SCCA (Mai & Zhang, 2019) SCCA-I (Hardoon & Shawe-Taylor, 2011) SCCA-II (Gao et al., 2017) mod-SCCA (Suo et al., 2017) ℓ -CCA 0|(1.038,1.067) (NA,NA) (1.382,1.357) (0.213,0.296) (0.173,0.218) (0.101,0.079)|(1.115,0.943) (NA,NA) (1.351,1.299) (0.145,0.109) (0.136,0.098) (0.098,0.072)|(1.098,0.890) (NA,NA) (1.219,1.186) (0.110,0.088) (0.109,0.086) (0.026,0.039)|
||Model III- Sparse Inverse covariance|||
|PMA (Witten et al., 2009) IP-SCCA (Mai & Zhang, 2019) SCCA-I (Hardoon & Shawe-Taylor, 2011) SCCA-II (Gao et al., 2017) mod-SCCA (Suo et al., 2017) ℓ -CCA 0|(0.930,1.050) (0.654,0.653) (1.375,0.966) (0.129,0.190) (0.092,0.149) (0.108,0.103)|(0.670,0.450) (0.092,0.091) (1.041,0.502) (0.069,0.062) (0.068,0.059) (0.026,0.036)|(0.760,0.580) (0.282,0.285) (0.985,0.364) (0.051,0.047) (0.050,0.044) (0.009,0.005)|


Model I- Identity covariance

(N, D[x], D[y]) (400,800,800) (500,600,600) (700,1200,1200)

PMA (Witten et al., 2009) (1.170,1.170) (0.850,0.850) (1.090,1.090)

IP-SCCA (Mai & Zhang, 2019) (1.658,1.647) (1.051,1.051) (1.544,1.542)

SCCA-I (Hardoon & Shawe-Taylor, 2011) (1.602,1.140) (1.143,0.282) (1.160,0.181)

SCCA-II (Gao et al., 2017) (0.060,0.066) (0.053,0.057) (0.045,0.043)

mod-SCCA (Suo et al., 2017) (0.056,0.062) (0.05,0.056) (0.045,0.043)

_ℓ0-CCA_ (0.003,0.009) (0.002,0.002) (0.001,0.002)

Model II- Toeplitz covariance

PMA (Witten et al., 2009) (1.038,1.067) (1.115,0.943) (1.098,0.890)

IP-SCCA (Mai & Zhang, 2019) (NA,NA) (NA,NA) (NA,NA)

SCCA-I (Hardoon & Shawe-Taylor, 2011) (1.382,1.357) (1.351,1.299) (1.219,1.186)

SCCA-II (Gao et al., 2017) (0.213,0.296) (0.145,0.109) (0.110,0.088)

mod-SCCA (Suo et al., 2017) (0.173,0.218) (0.136,0.098) (0.109,0.086)

_ℓ0-CCA_ (0.101,0.079) (0.098,0.072) (0.026,0.039)

Model III- Sparse Inverse covariance

PMA (Witten et al., 2009) (0.930,1.050) (0.670,0.450) (0.760,0.580)

IP-SCCA (Mai & Zhang, 2019) (0.654,0.653) (0.092,0.091) (0.282,0.285)

SCCA-I (Hardoon & Shawe-Taylor, 2011) (1.375,0.966) (1.041,0.502) (0.985,0.364)

SCCA-II (Gao et al., 2017) (0.129,0.190) (0.069,0.062) (0.051,0.047)

mod-SCCA (Suo et al., 2017) (0.092,0.149) (0.068,0.059) (0.050,0.044)

_ℓ0-CCA_ (0.108,0.103) (0.026,0.036) (0.009,0.005)


Table 1: Evaluating the estimation quality of the canonical vectors φ and η. Each pair indicates
(eφ, eη), which are the estimation errors of φ and η respectively. We compare the proposed ℓ0CCA to other sparse CCA schemes considering three types of covariance matrices for generating X
and Y, and different dimensions (N, D[x], D[y]). The description of all three covariances appears in
Section 5.1. We highlight the smallest error obtained across all methods using boldface.

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
_T_
between the extracted representations, defined by ˆρ = ρ(φ[ˆ] **_X, ˆη[T]_** **_Y ). 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.

5.2 MULTI VIEW OF SPINNING PUPPETS

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.

Figure 3: Left: two samples from the spinning puppets videos. The videos capture the rotation of 3
desk puppets. Arrows indicate the spinning direction of each puppet. We use ℓ0-CCA to identify a
sparse subset of pixels that are correlated from both videos. Right: the values of the gates for each
video z[x] and z[y]. After training, the values of the gates are binary (i.e., {0, 1}), with 372 and 403
active gates for the left and right videos, respectively. ℓ0-CCA correctly select correlated subsets of
pixels that highlight the common puppet (Bulldog) in this example.

Figure 4: The correlated ℓ0-CCA representations (visualized using MVDM (Lindenbaum et al.,
2020)) of the Yoda+Bulldog video (left) and Bulldog+Bunny (right). We superimpose each embedding with 6 images corresponding to 6 points in the embedding spaces. The transformations
are based on the information described by pixels whose gates are active (presented in Fig. 3). The
resulting embeddings are correlated with each other, with a total correlation of ¯ρ = 1.99. The structure captured by the embeddings correctly represent the angular rotation of the Bulldog, which is
the common latent parameter in this experiment.


-----

Figure 5: Images from the coupled noisy MNIST dataset. In the bottom right of both panels, we
presents the active gates (white values within a green frame). There are 277 and 258 active gates for
view I and II, respectively.
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 ≪ _D[x], D[y]._

5.3 NOISY MNIST

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 non-_
linear 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 ℓ0DCCA.

5.4 SEISMIC EVENT CLASSIFICATION

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 ∈ R[1157] 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
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).

1https://bigsoundbank.com/search?q=car


-----

Figure 6: Top: Clean sample sonograms of an explosion based on the E and N channels (left and
right, respectively). Arrows highlight the Primary (P) and Secondary (S) waves caused by the explosion. Middle: noisy sonograms generated by adding sonograms of vehicle recordings. Bottom:
the active gates for both channels. Note that the gates are active at time-frequency bins which corre
|ond to the P and S waves|(see top left panel).|Col3|Col4|Col5|Col6|Col7|Col8|Col9|Col10|
|---|---|---|---|---|---|---|---|---|---|
|Method|Noisy MNIST (LeCun et al., 2010)|||Seismic (Lindenbaum et al., 2018)|||METABRIC (Curtis et al., 2012)|||
||MI|KM (%)|SVM (%)|MI|KM (%)|SVM (%)|MI|KM (%)|SVM (%)|
|Raw Data PCA (Pearson, 1901) CCA (Chaudhuri et al., 2009) mod-SCCA (Suo et al., 2017) SCCA-HSIC (Uurtio et al., 2018) KCCA (Bach & Jordan, 2002) grad-KCCA (Uurtio et al., 2019) multiview-ICA (Richard et al., 2020) NCCA (Michaeli et al., 2016) DCCA (Andrew et al., 2013) DCCAE (Wang et al., 2015b) ℓ0-CCA (linear) ℓ0-DCCA (non-linear)|0.130 0.130 1.290 0.342 NA 0.943 NA 1.750 1.030 1.970 1.940 1.73 2.05|16.6 16.6 66.4 23.9 NA 50.2 NA 88.0 47.5 93.2 91.8 87.1 95.4|86.6 89.3 75.8 63.1 NA 85.3 NA 90.0 77.2 93.2 94.0 88.4 95.5|0.001 0.002 0.003 0.610 0.003 0.006 0.005 0.748 0.700 0.830 0.92 0.85 0.97|35.7 38.8 38.1 71.7 38.7 38.4 40.9 90.1 86.8 94.9 97.0 93.7 98.1|41.3 41.3 40.4 86.9 49.5 92.5 41.4 94.2 91.4 94.6 97.0 94.9 97.2|0.58 0.08 0.20 0.12 NA 0.35 0.50 0.74 0.72 0.79 0.68 0.82 0.88|36.5 19.2 20.7 21.0 NA 30.8 32.6 44.7 48.7 45.2 42.9 49.5 50.3|63.8 23.6 24.1 26.0 NA 61.3 47.8 62.8 63.7 72.1 69.0 70.7 74.1|


Noisy MNIST (LeCun et al., 2010) Seismic (Lindenbaum et al., 2018) METABRIC (Curtis et al., 2012)

Method MI KM (%) SVM (%) MI KM (%) SVM (%) MI KM (%) SVM (%)

Raw Data 0.130 16.6 86.6 0.001 35.7 41.3 0.58 36.5 63.8

PCA (Pearson, 1901) 0.130 16.6 89.3 0.002 38.8 41.3 0.08 19.2 23.6

CCA (Chaudhuri et al., 2009) 1.290 66.4 75.8 0.003 38.1 40.4 0.20 20.7 24.1

mod-SCCA (Suo et al., 2017) 0.342 23.9 63.1 0.610 71.7 86.9 0.12 21.0 26.0

SCCA-HSIC (Uurtio et al., 2018) NA NA NA 0.003 38.7 49.5 NA NA NA

KCCA (Bach & Jordan, 2002) 0.943 50.2 85.3 0.006 38.4 92.5 0.35 30.8 61.3

grad-KCCA (Uurtio et al., 2019) NA NA NA 0.005 40.9 41.4 0.50 32.6 47.8

multiview-ICA (Richard et al., 2020) 1.750 88.0 90.0 0.748 90.1 94.2 0.74 44.7 62.8

NCCA (Michaeli et al., 2016) 1.030 47.5 77.2 0.700 86.8 91.4 0.72 48.7 63.7

DCCA (Andrew et al., 2013) 1.970 93.2 93.2 0.830 94.9 94.6 0.79 45.2 72.1

DCCAE (Wang et al., 2015b) 1.940 91.8 94.0 0.92 97.0 97.0 0.68 42.9 69.0

_ℓ0-CCA (linear)_ 1.73 87.1 88.4 0.85 93.7 94.9 0.82 49.5 70.7

_ℓ0-DCCA (non-linear)_ **2.05** **95.4** **95.5** **0.97** **98.1** **97.2** **0.88** **50.3** **74.1**


Table 2: Evaluation of correlated embedding extracted from the Noisy MNIST, Seismic, and
METABRIC (cancer type) datasets. The representation extracted by ℓ0-DCCA leads to higher clustering accuracy (KM), and classification accuracy (SVM) compared with several baselines. MI
indicates the Mutual Information between cluster assignments and sample labels. We use NA to
denote simulations which did not converge.
5.5 CANCER SUB-TYPE CLASSIFICATION

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).

6 CONCLUSION

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.

7 ACKNOWLEDGMENT

The work of YK was supported by the National Institutes of Health R01GM131642, UM1PA05141,
U54AG076043, P50CA121974, and U01DA053628.


-----

REFERENCES

Fares Al-Shargie, Tong Boon Tang, and Masashi Kiguchi. Assessment of mental stress effects on
prefrontal cortical activities using canonical correlation analysis: an fnirs-eeg study. Biomedical
_optics express, 8(5):2583–2598, 2017._

Galen Andrew, Raman Arora, Jeff Bilmes, and Karen Livescu. Deep canonical correlation analysis.
In International Conference on Machine Learning, pp. 1247–1255, 2013.

Francis R Bach and Michael I Jordan. Kernel independent component analysis. Journal of machine
_learning research, 3(Jul):1–48, 2002._

Yoshua Bengio, Aaron Courville, and Pascal Vincent. Representation learning: A review and new
perspectives. IEEE transactions on pattern analysis and machine intelligence, 35(8):1798–1828,
2013.

Adrian Benton, Huda Khayrallah, Biman Gujral, Dee Ann Reisinger, Sheng Zhang, and Raman
Arora. Deep generalized canonical correlation analysis. In Proceedings of the 4th Workshop on
_Representation Learning for NLP (RepL4NLP-2019), pp. 1–6, 2019._

Sarah M Bernhardt, Pallave Dasari, David Walsh, Amanda R Townsend, Timothy J Price, and
Wendy V Ingman. Hormonal modulation of breast cancer gene expression: implications for
intrinsic subtyping in premenopausal women. Frontiers in oncology, 6:241, 2016.

Jia Cai, Wei Dan, and Xiaowei Zhang. l0-based sparse canonical correlation analysis with application to cross-language document retrieval. Neurocomputing, 329:32–45, 2019.

J. D. Carroll. Generalization of canonical correlation analysis to three or more sets of variables.
_Proceedings of the 76th annual convention of the American Psychological Association, 3:227–_
228, 1968.

Xiaobin Chang, Tao Xiang, and Timothy M Hospedales. Scalable and effective deep cca via soft
decorrelation. In Proceedings of the IEEE Conference on Computer Vision and Pattern Recogni_tion, pp. 1488–1497, 2018._

Kamalika Chaudhuri, Sham M Kakade, Karen Livescu, and Karthik Sridharan. Multi-view clustering via canonical correlation analysis. In Proceedings of the 26th annual international conference
_on machine learning, pp. 129–136. ACM, 2009._

Xi Chen and Han Liu. An efficient optimization algorithm for structured sparse cca, with applications to eqtl mapping. Statistics in Biosciences, 4(1):3–26, 2012.

Zhiwen Chen, Steven X Ding, Tao Peng, Chunhua Yang, and Weihua Gui. Fault detection for nongaussian processes using generalized canonical correlation analysis and randomized algorithms.
_IEEE Transactions on Industrial Electronics, 65(2):1559–1567, 2017._

Christina Curtis, Sohrab P Shah, Suet-Feung Chin, Gulisa Turashvili, Oscar M Rueda, Mark J Dunning, Doug Speed, Andy G Lynch, Shamith Samarajiwa, Yinyin Yuan, et al. The genomic and
transcriptomic architecture of 2,000 breast tumours reveals novel subgroups. Nature, 486(7403):
346, 2012.

Sarah-Jane Dawson, Oscar M Rueda, Samuel Aparicio, and Carlos Caldas. A new genome-driven
integrated classification of breast cancer and its implications. The EMBO journal, 32(5):617–628,
2013.

Chao Gao, Zongming Ma, Harrison H Zhou, et al. Sparse cca: Adaptive estimation and computational barriers. Annals of Statistics, 45(5):2074–2101, 2017.

Edward I George and Robert E McCulloch. Variable selection via gibbs sampling. Journal of the
_American Statistical Association, 88(423):881–889, 1993._

Gregory Gundersen, Bianca Dumitrascu, Jordan T Ash, and Barbara E Engelhardt. End-to-end
training of deep probabilistic cca on paired biomedical observations. In Uncertainty in artificial
_intelligence, 2019._


-----

David R Hardoon and John Shawe-Taylor. Sparse canonical correlation analysis. Machine Learning,
83(3):331–353, 2011.

Paul Horst. Generalized canonical correlations and their applications to experimental data. Journal
_of Clinical Psychology, 17(4):331–347, 1961._

Harold Hotelling. Relations between two sets of variates. Biometrika, 28(3/4):321–377, 1936.

Soham Jana, Henry Li, Yutaro Yamada, and Ofir Lindenbaum. Support recovery with stochastic
gates: Theory and application for linear models. arXiv preprint arXiv:2110.15960, 2021.

Eric Jang, Shixiang Gu, and Ben Poole. Categorical reparameterization with gumbel-softmax. arXiv
_preprint arXiv:1611.01144, 2016._

Mahdi Karami and Dale Schuurmans. Deep probabilistic canonical correlation analysis. In Pro_ceedings of the AAAI Conference on Artificial Intelligence, volume 35, pp. 8055–8063, 2021._

J. R. Kettering. Canonical analysis of several sets of variables. Biometrika, 58(3):433–451, 12
[1971. ISSN 0006-3444. doi: 10.1093/biomet/58.3.433. URL https://doi.org/10.1093/](https://doi.org/10.1093/biomet/58.3.433)
[biomet/58.3.433.](https://doi.org/10.1093/biomet/58.3.433)

Arto Klami, Seppo Virtanen, and Samuel Kaski. Bayesian canonical correlation analysis. Journal
_of Machine Learning Research, 14(4), 2013._

Lynn Kuo and Bani Mallick. Variable selection for regression models. Sankhy¯a: The Indian Journal
_of Statistics, Series B, pp. 65–81, 1998._

Yann LeCun, Corinna Cortes, and CJ Burges. Mnist handwritten digit database. ATT Labs [Online].
_Available: http://yann.lecun.com/exdb/mnist, 2, 2010._

Roy R Lederman and Ronen Talmon. Learning the geometry of common latent variables using
alternating-diffusion. Applied and Computational Harmonic Analysis, 44(3):509–536, 2018.

Ofir Lindenbaum, Yuri Bregman, Neta Rabin, and Amir Averbuch. Multiview kernels for lowdimensional modeling of seismic events. IEEE Transactions on Geoscience and Remote Sensing,
56(6):3300–3310, 2018.

Ofir Lindenbaum, Arie Yeredor, Moshe Salhov, and Amir Averbuch. Multi-view diffusion maps.
_Information Fusion, 55:127–149, 2020._

Ofir Lindenbaum, Uri Shaham, Erez Peterfreund, Jonathan Svirsky, Nicolas Casey, and Yuval
Kluger. Differentiable unsupervised feature selection based on a gated laplacian. Advances in
_Neural Information Processing Systems, 34, 2021._

Christos Louizos, Max Welling, and Diederik P Kingma. Learning sparse neural networks through
_l 0 regularization. arXiv preprint arXiv:1712.01312, 2017._

Chris J Maddison, Andriy Mnih, and Yee Whye Teh. The concrete distribution: A continuous
relaxation of discrete random variables. arXiv preprint arXiv:1611.00712, 2016.

Qing Mai and Xin Zhang. An iterative penalized least squares approach to sparse canonical correlation analysis. Biometrics, 75(3):734–744, 2019.

Roderick P McDonald. A unified treatment of the weighting problem. Psychometrika, 33(3):351–
381, 1968.

James John McKeon. Canonical analysis: Some relations between canonical correlation, factor
_analysis, discriminant function analysis, and scaling theory. Number 13. Psychometric Society,_
1966.

Tomer Michaeli, Weiran Wang, and Karen Livescu. Nonparametric canonical correlation analysis.
In International Conference on Machine Learning, pp. 1967–1976, 2016.


-----

Elena Parkhomenko, David Tritchler, and Joseph Beyene. Sparse canonical correlation analysis
with application to genomic data integration. Statistical applications in genetics and molecular
_biology, 8(1), 2009._

Karl Pearson. Principal components analysis. The London, Edinburgh, and Dublin Philosophical
_Magazine and Journal of Science, 6(2):559, 1901._

Harold Pimentel, Zhiyue Hu, and Haiyan Huang. Biclustering by sparse canonical correlation analysis. Quantitative Biology, 6(1):56–67, 2018.

Nicholas G Polson and Lei Sun. Bayesian l 0-regularized least squares. Applied Stochastic Models
_in Business and Industry, 35(3):717–731, 2019._

Hugo Richard, Luigi Gresele, Aapo Hyv¨arinen, Bertrand Thirion, Alexandre Gramfort, and Pierre
Ablin. Modeling shared responses in neuroimaging studies through multiview ica. arXiv preprint
_arXiv:2006.06635, 2020._

Moshe Salhov, Ofir Lindenbaum, Yariv Aizenbud, Avi Silberschatz, Yoel Shkolnisky, and Amir
Averbuch. Multi-view kernel consensus for data analysis. Applied and Computational Harmonic
_Analysis, 49(1):208–228, 2020._

Shai Shalev-Shwartz, Nathan Srebro, and Tong Zhang. Trading accuracy for sparsity in optimization
problems with sparsity constraints. SIAM Journal on Optimization, 20(6):2807–2832, 2010.

Omid S Solari, James B Brown, and Peter J Bickel. Sparse canonical correlation analysis via concave
minimization. arXiv preprint arXiv:1909.07947, 2019.

Xiaotong Suo, Victor Minden, Bradley Nelson, Robert Tibshirani, and Michael Saunders. Sparse
canonical correlation analysis. arXiv preprint arXiv:1705.10865, 2017.

Arthur Tenenhaus and Michel Tenenhaus. Regularized generalized canonical correlation analysis.
_Psychometrika, 76(2):257, 2011._

Arthur Tenenhaus, Cathy Philippe, Vincent Guillemot, Kim-Anh Le Cao, Jacques Grill, and Vincent
Frouin. Variable selection for generalized canonical correlation analysis. Biostatistics, 15(3):569–
[583, 02 2014. ISSN 1465-4644. doi: 10.1093/biostatistics/kxu001. URL https://doi.org/](https://doi.org/10.1093/biostatistics/kxu001)
[10.1093/biostatistics/kxu001.](https://doi.org/10.1093/biostatistics/kxu001)

Arthur Tenenhaus, Cathy Philippe, and Vincent Frouin. Kernel generalized canonical correlation
analysis. Computational Statistics and Data Analysis, 90:114–131, 2015. ISSN 0167-9473. doi:
https://doi.org/10.1016/j.csda.2015.04.004. [URL https://www.sciencedirect.com/](https://www.sciencedirect.com/science/article/pii/S0167947315001012)
[science/article/pii/S0167947315001012.](https://www.sciencedirect.com/science/article/pii/S0167947315001012)

Bruce Thompson. Canonical correlation analysis. Encyclopedia of statistics in behavioral science,
2005.

George Tucker, Andriy Mnih, Chris J Maddison, Dieterich Lawson, and Jascha Sohl-Dickstein. Rebar: Low-variance, unbiased gradient estimates for discrete latent variable models. arXiv preprint
_arXiv:1703.07370, 2017._

Viivi Uurtio, Sahely Bhadra, and Juho Rousu. Sparse non-linear cca through hilbert-schmidt independence criterion. In 2018 IEEE International Conference on Data Mining (ICDM), pp. 1278–
1283. IEEE, 2018.

Viivi Uurtio, Sahely Bhadra, and Juho Rousu. Large-scale sparse kernel canonical correlation analysis. In International Conference on Machine Learning, pp. 6383–6391. PMLR, 2019.

John P Van de Geer. Linear relations amongk sets of variables. Psychometrika, 49(1):79–94, 1984.

Sandra Waaijenborg, Philip C Verselewel de Witt Hamer, and Aeilko H Zwinderman. Quantifying
the association between gene expressions and dna-markers by penalized canonical correlation
analysis. Statistical applications in genetics and molecular biology, 7(1), 2008.


-----

Weiran Wang, Raman Arora, Karen Livescu, and Jeff Bilmes. On deep multi-view representation
learning. In Proceedings of the 32nd International Conference on International Conference on
_Machine Learning - Volume 37, ICML’15, pp. 1083–1092. JMLR.org, 2015a._

Weiran Wang, Raman Arora, Karen Livescu, and Jeff Bilmes. On deep multi-view representation
learning. In International conference on machine learning, pp. 1083–1092. PMLR, 2015b.

Weiran Wang, Raman Arora, Karen Livescu, and Nathan Srebro. Stochastic optimization for deep
cca via nonlinear orthogonal iterations. In 2015 53rd Annual Allerton Conference on Communi_cation, Control, and Computing (Allerton), pp. 688–695. IEEE, 2015c._

Weiran Wang, Xinchen Yan, Honglak Lee, and Karen Livescu. Deep variational canonical correlation analysis. arXiv preprint arXiv:1610.03454, 2016.

Ami Wiesel, Mark Kliger, and Alfred O Hero III. A greedy approach to sparse canonical correlation
analysis. arXiv preprint arXiv:0801.2748, 2008.

Ronald J Williams. Simple statistical gradient-following algorithms for connectionist reinforcement
learning. Machine learning, 8(3-4):229–256, 1992.

Daniela M Witten, Robert Tibshirani, and Trevor Hastie. A penalized matrix decomposition, with
applications to sparse principal components and canonical correlation analysis. Biostatistics, 10
(3):515–534, 2009.

Yutaro Yamada, Ofir Lindenbaum, Sahand Negahban, and Yuval Kluger. Feature selection using
stochastic gates. In Proceedings of Machine Learning and Systems 2020, pp. 8952–8963. 2020.

Mingzhang Yin, Nhat Ho, Bowei Yan, Xiaoning Qian, and Mingyuan Zhou. Probabilistic best subset
selection via gradient-based optimization. arXiv e-prints, pp. arXiv–2006, 2020.

Kosuke Yoshida, Junichiro Yoshimoto, and Kenji Doya. Sparse kernel canonical correlation analysis
for discovery of nonlinear interactions in high-dimensional data. BMC bioinformatics, 18(1):1–
11, 2017.

Yu Zhang, Guoxu Zhou, Jing Jin, Yangsong Zhang, Xingyu Wang, and Andrzej Cichocki. Sparse
bayesian multiway canonical correlation analysis for eeg pattern recognition. Neurocomputing,
225:103–110, 2017.

Mingyuan Zhou, Haojun Chen, John W Paisley, Lu Ren, Guillermo Sapiro, and Lawrence Carin.
Non-parametric bayesian dictionary learning for sparse image representations. In NIPS, volume 9,
pp. 2295–2303. Citeseer, 2009.

Wan Zhu, Longxiang Xie, Jianye Han, and Xiangqian Guo. The application of deep learning in
cancer prognosis prediction. Cancers, 12(3):603, 2020.

A APPENDIX

B ADDITIONAL EXPERIMENTAL DETAILS

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 of X 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.

Figure 7: Regularization path of ℓ0-CCA on data generated from the linear model I (described in
Section 5.1). Values on the left y-axis (green) represent the sum of active gates (by expectation).
Values on the right y-axis (blue) represent the empirical correlation between the estimated representations, i.e. ˆρ = **_φ[ˆ]T X_** _T Y ˆη, 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],
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.

Figure 8: k-means and SVM classification accuracy (left) and mutual information score (right) vs.
the number of selected features.

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

2https://scikit-learn.org/stable/modules/generated/sklearn.cross decomposition.CCA.html
3https://gist.github.com/yuyay/16ce4914683da30f87d0
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


-----

Figure 9: Classification accuracy on the noisy seismic data. Performance is evaluated using linear
SVM in the 3 dimensional embedding. Comparing performance of ℓ[0]-DCCA for different levels of
sparsity, and using linear and nonlinear activation (tanh).

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).

C LOW DIMENSIONAL SYNTHETIC EXAMPLE

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 ∈ R[D][×][N], using the following

distribution **_X_** _N_ **0** _,_ **Σx** **Σxy** _, where Σxy = ρ0Σx(φη[T]_ )Σy. With Σ = I _D._
**_Y_** _∼_ **0** **Σyx** **Σy**
   []    

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, I _Dσ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_ .

|Col1|(D = 20, N = 4000)|Col3|Col4|
|---|---|---|---|
|σ N|0|1|2|
|CCA SCCA-II (Gao et al., 2017) ℓ -CCA 0|(0.0098,0.0099) (0.0013,0.0008) (0.0005,0.0004)|(0.034,0.036) (0.0131,0.0097) (0.0065,0.0044)|(0.256,0.258) (0.6003,0.5881) (0.1075,0.0633)|


(D = 20, N = 4000)

_σN_ 0 1 2

CCA (0.0098,0.0099) (0.034,0.036) (0.256,0.258)

SCCA-II (Gao et al., 2017) (0.0013,0.0008) (0.0131,0.0097) (0.6003,0.5881)

_ℓ0-CCA_ (0.0005,0.0004) (0.0065,0.0044) (0.1075,0.0633)


Table 3: Evaluating the estimation quality of the canonical vectors φ and η. Each pair indicates
(eφ, eη), which are the estimation errors of φ and η respectively. We compare the proposed ℓ0CCA to standard CCA and to SCCA-II (Gao et al., 2017). Different columns indicate different
additive noise levels.

D CANCER SUB-TYPE CLASSIFICATION

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


-----

**Algorithm 1 ℓ0-DCCA**

**Input: Coupled data, {X, Y }, regularization parameters λ[x], λ[y], number of epochs T** .
Initialize the gate parameters: µ[x]i [= 0][.][5][ for][ i][ = 1][, . . ., D][x][, and][ µ]i[y] [= 0][.][5][ for][ i][ = 1][, . . ., D][y][.]

**for t = 1 to T do**

Sample a stochastic gate (STG) vectors z[x], z[y] defined based on equation 4.
Apply the STG to the data **_X[ˆ]_** = z[x] _⊙_ _X and_ **_Y[ˆ] = z[y]_** _⊙_ _Y ._

Compute the loss L = E _−_ _ρ¯(f_ **( X[ˆ]** ), g( Y[ˆ] )) + λ[x]∥z[x]∥0 + λ[y]∥z[y]∥0, where the calculation of
the total correlation ¯ρ is based on Section 3.4. Where f **( X[ˆ]** ) = f **(z[x]** **_X_** **_θ[x])_** R[d][×][N], and
  **_⊙_** **_|_** _∈_
**_gUpdate( Y[ˆ] ) = θ g[x](z=[y] θ⊙[x]_** **_Y |γθ[y])x ∈L,R[d][×]θ[N][y]_** .= θ[y] _γ_ _y_ _L,_

**_θ_** **_θ_**
_−_ _∇_ _−_ _∇_
**_π[x]_** = π[x] _γ_ **_π[x]_** _L,_ **_π[y]_** = π[y] _γ_ **_π[y]_** _L._
_−_ _∇_ _−_ _∇_

**end for**
Optional: apply linear CCA on f **( X[ˆ]** ), and g( Y[ˆ] )) to remove arbitrary scalings.


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.

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.

E PSEUDOCODE

In Algorithm 1 we provide a pseudocode description of the proposed approach.

F RUN TIME ANALYSIS

Our method requires computing covariance matrices for kx and ky variables with N samples at each
iteration. This scales like O(Nkx[2] [+][ Nk]y[2][)][. 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


-----

Figure 10: Run time evaluation of the proposed approach. We measure the average runtime (over
100 runs) for different values of N and D[x], D[y].

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 (Tenenhaus 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 X _[k]_ the corresponding k-th view 1 ≤ _k ≤_ _K. The Deep_
Generalized CCA that was proposed in (Benton et al., 2019) aims to solve


**_G_** **_U_** _[T]k_ _[f]_ [(][X] _[k][)][∥][F]_ _s.t._ **_GG[T]_** = I, (6)
_∥_ _−_
_k=1_

X


min
**_G∈R[n][×][d],U k∈R[d][×][d]_**


where G is the desired shared representation, and U k is a linear transformation from the k-th deep
network output, f (X _[k]), 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 z[k] for view X _[k]. The generalized_
_ℓ0-DCCA loss is defined by_


**_G_** **_U_** _[T]k_ _[f]_ [( ˆ]X _[k])_ _F + λk_ **_z[k]_** 0 _s.t._ **_GG[T]_** = I, (7)
_∥_ _−_ _∥_ _∥_ _∥_
_k=1_

X


min
**_G∈R[n][×][d],U k∈R[d][×][d],z[k]∈R[dk]_**


where f **( X[ˆ][k]) = f** **(z[k]** **_⊙_** **_X_** **_[k]|θ[x]) ∈_** R[d][×][N] and λk is the regularization parameter that control
the sparsity of k-th selected gate vector z[k]. The optimization problem in Eq 7 can be solved using
standard optimizers such as gradient descent. The initialization of G and U k plus the analysis of
the above suggested generalized gated canonical correlation are left for further research.

H PROOF OF THEOREM 1

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] _⊙_ **_s[x]_** and β = θ[y] _⊙_ **_s[y]. 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_**
solution to Eq. 3 which leads to the same objective value.αˆ = ˆa ⊙ **_sˆ[x]_** and **_β[ˆ] = b[ˆ] ⊙_** **_sˆ[y], where ˆπi[x]_** [= 1][ if][ ˆ]ai ̸= 0 and πi[x] [= 0][ otherwise, we get a feasible]

Now, we want to show that the optimal solution to Eq. 3 is also a feasible solution to Eq. 2.
_x_ _y_
Denoting the optimal solution to Eq. 3 by ˆα and **_β[ˆ], with ˆα = θ[ˆ]_** **_sˆ[x]_** and **_β[ˆ] = θ[ˆ]_** **_sˆ[y],_**
_⊙_ _⊙_
if p(s[x]|π[x]) or p(s[y]|π[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(s[x]|π[x])_**
and p(s[y]|π[y]) are not point mass densities, we will now show by contradiction that all the points in
the support of p(s[x]|π[x]) (or p(s[y]|π[y])) would lead to the same objective values. Assume without
loss of generality that {s[x]1 _[, ....,][ s][x]L[} ∈]_ [supp][[][p][(][s][x][|]π[ˆ] _[x])], if not all values in the support lead to the_
same objective value, this means that there exist s[x]l _[,][ s]m[x]_ [such that][ L][(][s]l[x][)][ < L][(][s]m[x] [)][, where]

_L(s[x]) =_ _ρ((θ[ˆ]x_ **_sx)T X, (θˆy_** **_sˆ[y])[T]_** **_Y ) + λ[x]_** **_s[x]_** 0 + λ[y] **_sˆ[y]_** 0 _,_
_−_ **_⊙_** **_⊙_** _∥_ _∥_ _∥_ _∥_

_x_ _y_
which is the objective from Eq. 3 with optimal values of **_θ[ˆ]_** _, ˆθ_ _, ˆs[y]. Now if we set ¯π[x]_ = s[x]l [we]
would obtain Esx∼p(sx| ¯π[x])[L][(][s][x][)][ <][ E]s[x]∼p(s[x]| ˆπ[x])[L][(][s][x][)][ which contradicts the assumption that][ ˆ]π[x]

is optimal. This means that all points in supp[p(s[x]|πˆ _[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.

I LIMITATIONS AND FUTURE DIRECTIONS

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) by O(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.


-----

