Stochastic Canonical Correlation Analysis: A Riemannian Approach
Authors: Zihang Meng
Rudrasis Chakraborty
Vikas Singh
Venue:  (2021)
Keywords: CCA
streaming
differential geometry
DeepCCA
fairness

Abstract:  We present an efficient stochastic algorithm (RSG+) for canonical correlation analysis (CCA) derived via a differential geometric perspective of the underlying optimization task. We show that exploiting the Riemannian structure of the problem reveals natural strategies for modified forms of manifold stochastic gradient descent schemes that have been variously used in the literature for numerical optimization on manifolds. Our developments complement existing methods for this problem which either require $O(d^3)$ time complexity per iteration with $O(\frac{1}{\sqrt{t}})$ convergence rate (where $d$ is the dimensionality) or only extract the top $1$ component with $O(\frac{1}{t})$ convergence rate. In contrast, our algorithm achieves $O(d^2k)$ runtime complexity per iteration for extracting top $k$ canonical components with $O(\frac{1}{t})$ convergence rate. We present our theoretical analysis as well as experiments describing the empirical behavior of our algorithm, including a potential application of this idea for training fair models where the label of protected attribute is missing or otherwise unavailable.

Under review as a conference paper at ICLR 2021
STOCHASTIC CANONICAL CORRELATION ANALYSIS: A RIEMANNIAN APPROACH
Anonymous authors Paper under double-blind review
ABSTRACT
We present an efficient stochastic algorithm (RSG+) for canonical correlation analysis (CCA) derived via differential geometric perspective of the underlying optimization task. We show that exploiting the Riemannian structure of the problem reveals natural strategies for modified forms of manifold stochastic gradient descent schemes that have been variously used in the literature for numerical optimization manifolds. Our developments complement existing methods for this problem which either require (ad)) time complexity per iteration with O convergence
rate (where d is the dimensionality) or only extract the top 1 component with ( convergence rate. In contrast, our algorithm achieves (dP))) runtime complexity per iteration for extracting top & canonical components with O(2) convergence rate. We present our theoretical analysis as well as experiments describing the empirical behavior of our algorithm, including potential application of this idea for training fair models where the label of protected attribute is missing or otherwise unavailable.
INTRODUCTION
Canonical correlation analysis (CCA) is a popular method for evaluating correlations between two sets of variables. is commonly used in unsupervised multi-view learning, where the multiple views of the data may correspond to image, text, audio and so on Rupnik & Shawe-TTaylor (2010); Chaudhuri et (2009) : Luo al. (015)). Classical CCA formulations have also been extended to leverage advances in representation learning, for example, (ndrew et al. (2013) showed how the CCA can be interfaced with deep neural networks enabling modern use cases. Many results over the last few years have used CCA or its variants for problems including measuring representational similarity in deep neural networks Morcos et al. (2018), speech recognition Couture et al. (2019), etc.
The goal in CCA is to find linear combinations within two random variables X and Y which have maximum correlation with each other. Formally, the CCA problem is defined in the following way. Given pair of random variables, d, variatt random variable X and dy-variate random variable Y, with unknown joint probability distribution, find the projection matrices U E Rd.x and V E Rd xk with k < min (dy dy). such that the correlation is maximized:
maximize trace Ex" [TTY] V) s.t. E X U = Ix, E Y
(1)
Here, X,Y Y are samples of X and Y respectively. The objective function in ① is the expected cross- correlation in the projected space and the constraints specify that different canonical components should be decorrelated. Let Cx Ex [X and Cy Ey Ex TY] be the covariance matrices, and CXY = ( YY] denote cross-covariance. Let us define the whitened covariance T == Ox 1/1122xr Cxr 1/2 and t (and v& contains the top-k left (and right) singular vectors of T. It is 1/2 known Golub & Zha 1992) that the optimum of 1 is achieved at U -x2/22 V*
In practice, we may be given two views of N samples X E R Nxd, and Y E RR A natural approach to solving CCA is based on on the following sequence of steps. We first compute the empirical covariance and cross-covariance matrices, namely, Čx = 1/xX"X, Ōr = 1/NYTy and Čxx I/xXTY. We then calculate the empirical whitened cross-covariance matrix T, finally. compute U*, V* by applying k-truncated SVD to T.

---

Under review as a conference paper at ICLR 2021
Runtime and memory considerations. The above procedure is simple but is only feasible when the data matrices are small. But in most modern applications, not only are the datasets large but also the dimension d (let d max dz, d of each sample can be quite high, especially if representations are being learned using deep neural network models. As a result, the computational footprint of the foregoing algorithm can be quite high. This has motivated the study of stochastic optimization routines for solving CCA. Observe that in contrast to the typical settings where stochastic optimization schemes are most effective, the CCA objective does not decompose over samples in the dataset. Many efficient strategies have been proposed in the literature: for example, Ge et al. (2016); WWan et al. (2016) present Empirical Risk Minimization (ERM) models which optimize the empirical objective. More recently. Gao et al. (2019); [Bhatia et al. (2018) ; Arora et al. (2017) describe proposals that optimize the population objective. To summarize the approaches succinctly, if we are satisfied with identifying the top 1 component of CCA, effective schemes are available by utilizing either extensions of the Oja's rule (1982) to generalized eigenvalue problem Bhatia et al. (2018) or the alternating SVRG algorithm Gao et al. (2019). Otherwise, stochastic approach must make of explicit whitening operation which involves a cost of de for each iteration AAroa et al. (0177).
Observation. Most approaches either directly optimize [ instead reparametrized or regularized form Ge et al. (2016) Zhl & Li (2016); AAooa et al. (2017). Often, the search space for U and V corresponds the entire R (ignoring the constraints for the moment). But if the formulation could be cast in a form which involved approximately writing U and V as product of several matrices with nicer properties, we may obtain specialized routines which are tailored to exploit those properties. Such a reformulation is not difficult to derive- - where the matrices used to express U and V can be identified as objects that live in well studied geometric spaces. Then, utilizing the geometry of the space and borrowing relevant tools from differential geometry leads to an efficient approximate algorithm for top-k CCA which optimizes the population objective streaming fashion.
Contributions. (a) First, we re- pe-aameeterize the top-k CCA problem as an optimization problem on specific matrix manifolds, and show that it is equivalent to the original formulation in equation [ (b) Informed by the geometry of the manifold, we derive stochastic gradient descent algorithms for solving the ---arammtterzed problem with (d² cost per iteration and provide convergence rate guarantees. (c) This analysis provides direct mechanism to obtain an upper bound on the number of iterations needed to guarantee an error w.r.t. the population objective for the CCA problem. (d) The algorithm works in a streaming manner it easily scales large datasets and we do not need to assume access to the full dataset at the outset. (e) We present empirical evidence for both the standard CCA model and the DeepCCA setting Andrew et al. (2013), describing advantages and limitations.
2 STOCHASTIC CCA: REFORMULATION, ALGORITHM AND ANALYSIS
The formulation of Stochastic CCA and the subsequent optimization scheme will seek to utilize the geometry of the feasible set for computational gains. Specifically, we will use the following manifolds (please see Absil et al. (2007) for more details):
(a) Stiefel: St(p, n). The manifold consists of n x p, with < n, column orthonormal matrices, i.c., St(p, n) = {x R=*P| = Ip)]-
(b) Grassmanian: Gr(p, r2). The manifold consists of p-dimensional subspaces in R", with p < n.
(c) Rotations: SO(n) m). the manifold/group consists of n X n special orthogonal matrices, i.c.,
SO(n) = X E R x R=x"IxTX = xxT = In, det( X)= = 11-
We summarize certain geometric properties/operations for these manifolds in the Appendix but have been leveraged in recent works for other problems also Li et al. (2020): Rezende et al. (2020).
Let us recall the objective function for CCA as given in (1. We denote X E Nx as the matrix consisting of the samples (x:) drawn from a zero mean random variable X ~ X' and Y E R. Xdy denotes the matrix consisting of samples {y:] drawn from a zero mean random variable Y y. For notational and formulation simplicity, we assume that d, = d, = d in the remainder of the paper although the results hold for general d' and dy.
Let Cx. Cy be the covariance matrices of X, Y. Also, let Cxy be the cross-correlation matrix between X and Y. Then, we can write the CCA objective as
max U,V
F F= trace VX)
subject to UTCUU= It VTC,V= =
(2)
2

---

Under review as a conference paper at ICLR 2021
) is the matrix consisting of {u,} (Vv))) where (fu,) (vy}) are the canonical directions. The constraints in equation 2 are called whitening constraints.
Here, E R (V E R Rd
Let us define matrices Ü, V E Rdxk which lie on the Stiefel manifold, St(k, d). Also, let Su. S, E R denote upper triangular matrices and Qu, Qu E SO(k). We can rewrite the above equation and the constraint as follows.
A Reformulation for CCA
max S.V.S..S.RR.. U UuUQuS.: V=VQ.
F trace (UTCxyV)
(3a)
subject to
C²UXU - CyV =
(3b)
Ú. E d); Qui Qu SO(k); S,,, S, is upper triangular
Here, we will maximize 3a with respect to U, V, Su, Su, Qui and Qu satisfying equation 3b
Main adjustment from to 3): In (2) while U and V should decorrelate Cx and Cy respectively, the optimization/search unrestricted and treats them as arbitrary matrices. In contrast, equation additionally decomposes U and V as U = UQuS and V = VQ..,, with the components as structured matrices. Hence, the optimization is regularized.
The above adjustment raises two questions: (i) does there exist a non-empty feasible solution set for (3) (ii) if solution to (3) can be found (which we will describe later), how "good" is this solution for the CCA objective problem, i.e., for ?
Existence of feasible solution: We need to evaluate if the constraints in 3b can be satis-
fied at all. Observe that by using Ü to be the top-k principal directions of X. S,, be the 1/ -top-k eigen values of C'x and Qu to be any orthogonal matrix, we can easily satisfy the "whiten- ing constraint" and hence UQ, is feasible solution of U in 3 and similarly for V. From this starting point, we can optimize the objective while ensuring that we maintain feasibility.
Is the solution for equation 3 good approximation for equation 2:: We can show that under
some assumptions, the estimator for canonical correlation, i.c., solution of equation B is consistent, i... solves equation [ We will state this formally shortly.
Before characterizing the propertics of a solution for equation we first provide some additional intuition behind equation and describe how it helps computationally.
Intuition behind the decomposition U = UQ: A key observation is the following. Recall that by performing principal component analysis (PCA), the resultant projection matrix will exactly satisfy the decorrelation condition needed for the "whitening constraint" in equation 2 (projection matrix consists of the eigen-vectors of x). A natural question to ask is: Can we utilize streaming PCA algorithm to help us obtain an efficient streaming CCA algorithm? Let us assume that our estimate for canonical correlation directions, i.c., solutions of equation 3 lies in the principal subspace calculated above. If so, we can use the decomposition U = Ü A., (analogously for V), where contains the principal directions, i.e., St(k, d) and A, is full rank k x matrix containing the coefficients of the span. But maintaining the full rank constraint during optimization is hard, so we further decompose A, into A., Q, with Qu E SO(k); S, is upper triangular. Additionally, we ensure the diagonal of S,, to be non-zero to maintain full-rank of S., During optimization, we can maintain the non-zero diagonal entries by optimizing the log of the diagonal entries instead.
Why equation 3 helpss First, we note that CCA seeks to maximize the total correlation under the constraint that different components are decorrelated. The difficult part in the optimization is to ensure decorrelation, which leads to higher complexity in existing streaming CCA algorithms. On the contrary, equation R we separate equation Cataaiooonnto into finding the PCs and finding the linear coefficients for the span of principal directions. Then, by utilizing efficient streaming PCA algorithm, lower complexity can be achieved. We will defer describing the specific details of the optimization itself until the next sub-section. First, will show formally why substituting equation with equation 3a equation i3 sensible under some assumptions.
3

---

Under review as a conference paper at ICLR 2021
2.1 How TO USE THE REFORMULATION IN EQUATION 3P
We first start by stating some mild assumptions needed for the analysis.
Assumptions: (a) The random variables X ~ (N Mxp) and Y ~ (0 2nn) with Dx 3 cla and Σy cla for some 0. (b) The samples X and Y drawn from X' and y respectively have mean. (c) For given &o I d, Σn and En have non-zero top-k eigen values.
A high-level solution to optimize F in equation 3 summarized earlier.
Recall the following scheme which we briefly
(a) Initialize Ü, V. E St(k: d) as the top-k eigen vectors of Cx (//N)*XX and Cy
(1/))TTY respectively; Initialize Qu and Qu to be random SO(k) matrices; (h) Set S., and S, be diagonal matrices with the diagonal entries to be the square root of the inverse of the top-k eigen values (to satisfy upper-triangula property):
Observe that with this initialization, the constraints in equation 3b are satisfied. With a feasible solution for U and V in hand, we may optimize equation 3a while satisfying equation 3b The specific details of how this is done is not critical this point long as we assume that a suitable numerical optimization scheme exists and can be implemented.
With the component matrices, we can construct the solution U UQ.S and V VQ..sp
Why the solution makes sense? We now show how the presented solution, assuming access to an effective numerical procedure, approximates the CCA problem presented in equation 2 We formally state the result in the following theorem with sketch of proof (appendix includes the full proof) by first stating the following proposition and definition.
Definition 1. A random variable X is called sub-Gaut Gaussian if the norm given by |X|| inf (d > E (exp (mace(X² " x1/a2) & 2} is finite. Let U E BRdx then XU sub-Gaussian Ver
shynin 2017
Proposition (Reiß et al. (2020). Let X be a random variable which follows a sub-Gaussian distribution. Let X be the approximation of X E N (samples frawn.from X) with the top-k principal vectors. Let Cx be the covariance of X. Also, assume that x is the 8th eigen value of Cx exforr = 1,... d
and X; 2 deti for all i. Then, the PCA reconstruction error, denoted by & = X 811 (in the Frobenius norm sense) can be upper bounded as follows
1/12/2
& S min 22 All,
where A = Cx - Üx-
d AA+1
The aforementioned proposition suggests that the error between the data matrix X and the recon- structed data matrix X using the top-k principal vectors is bounded.
CCA objective denoted by F and F. The following theorem states that we can bound the error, = || F|l (proof is in the Appendix). The proof
Recall from and 3 that the value
includes hoper-bounding E by the reconstruction error of the data projected on the principal directions using Prop. [
Theorem 1. Using the hypothesis and assumptions above, the approximation error E = = Fl1 is bounded and goes to zero while the whitening constraints in equation 3D are satisfied.
Now, the only unresolved issue an optimization scheme for equation 3a that keeps the constraints in equation satisfied by leveraging the geometry of the feasible set.
2.2 How TO NUMERICALLY OPTIMIZE
(A
SATISPYING CONSTRAINTS (3)
Overview. We now describe how to maximize the formulation in equation 3a equation 3b with
respect to Ü. V. Qu. Qr. S2 and S.. We will first compute top-k principal vectors to get U and V. Then, we will use a gradient update rule to solve for Q.. Qu, S, and S, to improve the objective. Since all these matrices are "struccturdd", care must be taken to ensure that the matrices remain
4

---

Under review as a conference paper at ICLR 2021
on their respective manifolds which is where the geometry of the manifolds will offer desirable properties. We re-purpose a Riemannian stochastic gradient descent (RSGD) to do this, so call our algorithm RSG+. Of course, more sophisticated Riemannian optimization techniques can be substituted in. For instance, different Riemannian optimization methods are available in Absil et al. (2007) and optimization schemes for many manifolds are offered in PyManOpt Boumal et al. 2014))
The algorithm block is presented in Algorithm I (a direct implementable block for the algo- rithm including the expression for gradients is presented in the Appendix A.3).. Let Fpri = trace ()U) + trace (V CyV) be the contribution from the principal directions which we used to ensure the "whitening constraint". Let Fean trace (UT CxyV) be the contribution from the canonical correlation directions. The algorithm consists of four main blocks denoted by different colors, namely (a) the Red block deals with gradient calculation of the objective function where we calculate the top-k principal vectors (denoted by Fpo)) with respect to Ü. V: (b) the Green block describes calculation of the gradient corresponding to the canonical directions (denoted by Femn) with respect to U, V. Su. S.. Qu and Qu: (c) the Gray block combines the gradient computation from both Fpri and F Faan with respect to unknowns U, V, S.s Su, Qu and Qui and finally (d) the Blue block performs batch update of the canonical directions Feam using Riemannian gradient updates.
Gradient calculations. The gradient update for Ü, V is divided into two parts (a) The (Red block) gradient updates the "principal" directions (denoted by Vū:Fmi and V..Fme. which is specifically designed to satisfy the whitening constraint. Since this requires updating the principal subspaces, so, the gradient descent needs proceed the manifold of subspaces, i.e., on the Grassmannian.
(b) The (green block) gradient from the objective function in equation 3 is denoted by Vi and V Fcan. In order to ensure that the Riemannian gradient update for Ú and V stays on the manifold St(k,, d). we need to make sure that the gradients, i.e., Vū Fean and V Fean lies in the tangent space of St(k, d). In order to do that, we need to first calculate the Euclidean gradient and then project on to the tangent space of St(k,d) d).
The gradient updates for Qu, Qu, Sw, S.p are given in the green block, denoted by VQ. Fem. VQ. Fean. Vs, and V Fea Note that unlike the previous step, this gradient only has components from canonical correlation computation. As before, this step requires first computing the Euclidean gradient and then projecting on to the tangent space of the underlying Riemannian manifolds involved, i.e., SO(k) k) and the space of upper triangular matrices.
Finally, we get the gradient to update the canonical directions by combining the gradients which is shown in gray block. With these gradients we can perform a batch update as shown in the blue block. Using convergence results presented next in Propositions 2 this scheme can be shown (under some assumptions) to approximately optimize the CCA objective in equation 2
We can now move to the convergence properties of the algorithm. We present two results stating the asymptotic proof of convergence for top-k principal vectors and canonical directions in the algorithm.
Proposition 2 (Chakraborty et al. (2020)). (Asymptotically) If the samples, X, are drawn from Gaussian distribution, then the gradient update rule presented in Step A returns an orthonormal basis - the top-k principal vectors of the covariance matrix Cx.
Proposition 3. (Bonnabel (2013) Consider a connected Riemannian manifold M with injectivity radius bounded from below by I > O. Assume that the sequence of step sizes (7) satisfy the condition (a) Σ 2 < do (b) Σ oc. Suppose {Ai} lie in a compact set K C -M. We also suppose that
3D > O such that, gA,
F. A
F
< D. Then VA, F + 0 and + cc.
Notice that in our problem, the manifold M can be Gr(p, n). St(p, n) or SO(p). Hence all the assumptions in Proposition (3 are satisfied if we guarantee the step sizes satisfy the aforementioned condition. One example of the step sizes that satisfies the property is 7
III
2.3 CONVERGENCE RATE AND COMPLEXITY OF THE RSG+ ALGORITHM
In this section, we describe the convergence rate and complexity of the algorithm proposed in Algorithm Observe that the key component of Algorithm is Riemannian gradient update. Let
5

---

Under review as a conference paper at ICLR 2021
A, be the generic entity needed to be updated in the algorithm using the Riemannian gradient update As+1 E. (~h V. where "/r is the step size at time step t. Also assume {At} C .M for a Riemannian manifold M. The following proposition states that under certain assumptions, the Riemannian gradient update has convergence rate of O (i).
Proposition 4. (Nemirovski et al. (2009); Bécigneul & Ganea (2018)) Let {Ac} lie inside a geodesic ball of radius less than the minimum of the injectivity radius and the strong convexity radius of M. Assume M to be geodesically complete Riemannian manifold with sectional curvature lower bounded by K < O. Moreover, assume that the step size {re diverges and the squared step size converges. Then, the Riemannian gradient descent update given by A ExpA. (~.
with bounded V A. F, i.e., | Fl1 S C' co some C 2 0, converges in the rate of O (t).
All Riemannian manifolds we used, i.c. Gr(k; d). St(kk d) and SO(k) are geodesically complete, and these manifolds have non-negative sectional curvatures, i.e. lower bounded by = 0. Now, as long as the Riemannian updates lie inside the geodesic ball of radius less than the minimum of injectivity and convexity radius, the convergence rate for RGD applies in setting.
Running time. To evaluate time complexity, we must look at the main compute-hese hhaae steps needed.
The basic modules are Exp and Exp maps for St(k, d). Gr(k, d) and SO(k) manifolds (see Table in the appendix). Observe that the complexity of these modules is influenced by the complexity of svd needed for the Exp map for the St and Gr manifolds. Our algorithm involves structured matrices of size X k and k k, so any matrix operation should not exceed cost of O (maax(d k. k²))), since in general > k. Specifically, the most expensive calculation is SVD of matrices of size d x k, which is (dd)) ()), see Golub Reinsch 1971). All other calculations are dominated by this term.
EXPERIMENTS
We first evaluate RSG+ for extracting top-k canonical components on three benchmark datasets and show that it performs favorably compared with Arora et al. (017)). Then, show that RSG+ also fits into feature learning in DeepCCA Andrew et al. (2013) and can scale to large feature dimensions where the non-stochastic method fails to. Finally we show that RSG+ can be used improve fairness of deep neural networks without needing labels of protected attributes during training.
Algorithm 1: Riemannian SGD based algorithm (RSG+) to compute canonical directions
Input: X E R. NXdx Y FERRR Rxd, &
Output: U V RR V R R
1 Initialize U, V.. Q.., S.; Partition X, into batches size B. Batch jth denoted by X, and Yj
2 for (1...., do
Gradient for top-k principal vectors: calculating V Vp.
1. Partition X, (Y₂) into L (L = _ blocks size d, x (dy x k);
2. Let the 1th block be denoted by Z/ (Z,'):
3. Orthogonalize each block and let the orthogonalized block be denoted by 27 (2
4. Let the subspace spanned by each 27 (and 27) 27 Gr(k, d.) (and 2/ Gr(k, d,));
5. Update Wy Fgn and V Fpi based on 2 and 2"" respectively.
Gradient from equation 3 calculating V Fen., Vp Vq,, Fen, Vq. Vs, Fen. Vs, F Calculation of Riemannian gradients of Ü, Qu. Qu S. and S. from equation i.e. objective from CCA.
F where, A c (Ū, V. Q.. S., The final gradients of Ü and V is a combination of gradient from objective for principal vectors and CCA; On the other hand, the gradients for Qu, Q. S., S. is from only CCA objective;
Gradient to update canonical directions: calculating V F
Batch update canonical directions
A Exp,
(AF where. is a generic entity: A E (Ū, V. Qu. Qu, Su,
end
s U UQ.S. and == VQQ...::

---

Under review as a conference paper at ICLR 2021
3.1 CCA ON FIXED DATASETS
Datasets and baseline. We conduct experiments three benchmark datasets (MNIST LeCun et al. (2010), Mediamill Snock et al. (006) and CIFAR-10 CIFAR-10(rrihhvskkk (2009) to evaluate the performance of RSG+ extract top-k canonical components. To best knowledge, Arora et al. (2017) is the only previous work which stochastically optimizes the population objective in a streaming fashion and can extract top-k components, so we compare our RSG+ with the matrix stochastic gradient
(MSG) method proposed in |rooa et al. (2017) (There are two methods proposed in Arora et al.
(2017) and we choose MSG because it performs better in the experiments of Arora al. 2017)) The details about the three datasets and how we process them are as follows:
WNISTLeCuus et al. (2010): MNIST contains grey-scale images of size 28 x 28. We use its full training set containing 60K images. Every image is split into left/right half, which are used as the two views. Mediamill Snoek et al. (2006): Mediamill contains around 25.8K paired features of videos and corresponding commentary of dimension 120, 101 respectively. CIFAR-10 Krizhevsky (2009): CIFAR-10 contains 60K 32 x 32 color images. Like MNIST, we split the images into left/right half and use them two views.
Evaluation metric. We choose to use Proportion of Correlations Captured (PCC) which is widely used Ma et al. (2015); Ge et al. (2166), partly due to its efficiency, especially for relatively large datasets. Let U E Rd, V E Rdy denote the estimated subspaces returned by RSG+, and U E Rd,xx x V E x denote the true canonical subspaces (all for top-k). The PCC is defined as PCC TCC(XÚ,YV) TCC(XU-,YYS).
where TCC is the sum of canonical correlations between two matrices.
Performance. See A.44 for the implementation deails. The performance in terms of PCC as function of # of seen samples (coming in streaming way) shown in Fig. and the runtime is reported in 10 times faster. One case where our RSG+ underperform Arora et al. (2017) is when the top-k eigenvalues dominated by the top-/ eigenvalues with < (Fig. lb on Mediamill dataset, the top-4 eigenvalues of the covariance matrix in view are: 8.61, 2.99, 1.15, 0.37. The first eigenvalue is dominantly large compared with the rest and our RSG+ performs better for A = and worse than Arora et al. (2017) for k 2, We also plot the runtime of RSG+ under different data dimension (set dg dy
A.5 Our RSG+ captures more correlation than MSG Arora et al. (2017) while being
and number of total samples sampled from joint gaussian distribution A.5
We implemented the method from Yger et al. (2012) and conduct experiments on the three datasets above. The results are shown in Table We tune the step size between [0.0001, 0.1] and / = 0.99 as used in their paper. On MNIST and MEDIAMILL, the method performs comparably with ours except & = 4 case on MNIST where it does not converge well. Since this algorithms also has d3 complexity, the runtime is 100x more than ours on MNIST and 2 more on Mediamill. On CIFAR 10, we fail to find a suitable step size for convergence.
3.2 CCA FOR DEEP FEATURE LEARNING
Background and motivation. A deep neural network (DNN) extension of CCA was proposed by Andrew et al. (2013) and has become popular in the multi view representation learning tasks. The idea is to learn deep neural network as the mapping from original data space to latent space where the canonical correlations are maximized. We refer the reader to Andrew et al. (2013) for details of the task. Since deep neural networks usually trained using SGD on mini-batches this requires
MNIST
MEDIAMILL
CIFAR
3.7
07.
06.
2.3
0.5
0.4
8.8
ba
0:3
-o-3 bad
0.2
-- had
0.1
0153 kut
43 and
355
of aarrgaee ssen
# carpplas seen
Missspllle seen
(a) on MNIST
(b) on Mediamill
(c) on CIFAR
Figure 1: Performance three datasets in terms of PCC function of # of seen samples.
7

---

Under review as a conference paper at ICLR 2021
getting estimate of CCA objective at every iteration in streaming fashion, thus our RSG+ can be natural choice here. We conduct experiments on a noisy version of MNIST dataset to evaluate RSG+.
Dataset. We follow Wang et al. (2013a) to construct a noisy ver- sion of MNIST: View is ran- domly sampled image which is first rescaled to [0, 1] and then rotated by a random angle from
Table 1: Results of Yger et al.
2012) (on CIFAR-10, our imple-
mentation of Y ger et al. (2012)
faces convergence issues).
MNIST
Mediamill
Performance
1:
k = 2
k 4
k
=
k 4
PCC
0.93
0.81
0.53
0.55
0.61
0.51
Time (s)
575.88
536.46
540.91
41.89
28.66
28.76
[-2 -1. View 2 is randomly sampled from the same class as view Then we add independent
uniform noise from [0, 1] to each pixel. Finally the image is truncated into [0, 1] to form the view 2.
Implementation details. We use a simple 2-layer MLP with ReLU nonlinearity, where the hidden dimension in the middle is 512 and the output feature dimension is d E {100, 500, 1000). After the network is trained on the CCA objective, we use a linear Support Vector Machine (SVM) to measure classification accuracy on output latent features. Andrew et al. (2013) uses the closed form CCA objective on the current batch directly, which costs Ocd)) memory and time for every iteration.
Table 2: Results of feature learning on MNIST.
N/A means fails to yield a result on our hardware.
(ccuracy(%5) d = 100
d = 500
= 1000
DeepCCA
80.57
N/A
N/A 86.39
Ours
79.79
84.09
Performance.
Table 2 shows that we get similar performance when d = 100 and can scale to large
latent dimensions d = 1000 while the batch method Andrew et al. (2013) encounters numerical difficulty on our GPU resources and Pytorch Paszke al. (2019) platform in performing eigen- decomposition of d x d matrix when d = 500, and becomes difficult if d is larger than 1000.
3.3 CCA FOR FAIRNESS
Table 3: Fairness results on CelebA. We applied CCA on three different layers in Resnet- respec- tively. See A.8 for positions of conv 0, 1,
Background and motivation. Fairness is be- coming an important issue to consider in the design of learning algorithms. A common strat- egy to make an algorithm fair is remove the influence of one/more protected attributes when training the models, see Lokhande al.] (2020). Most methods assume that the labels of pro- tected attributes are known during training but this may not always be possible. CCA enables considering slightly different setting, where we may not have par-sample protected attributes which may be sensitive or hard to obtain for third- sard-pparties Price & Cohen 2019). On the other hand, we assume that model trained to predict the protected attribute labels has been trained and is provided. For example, if the protected attribute is gender, we only assume that well trained classifier which is trained to predict gender from the samples is available rather than sample-wise gender values themselves. We next demonstrate that fairness of the model, using standard measures, can be improved via constraints on correlation values from CCA.
Accuracy(%%
DEO(%)
DDP(%)
Unconstrained Ours-convi Ours-cony1
76.3
22.3
4.8
76.5
17.4
1.4
77.7
15.3
3.2
Ours-conv2
75.9
22.0
2.8
Dataset. CelebA Wang et al. 20155) consists of 200K celebrity face images from the internet. There are up 40 labels, each of which is binary-valued. Here, we follow Lokhande et al. (2020) focus on the attactiveness attribute (which we want to train classifier to predict) and the gender is treated as "protected" since it may lead an unfair classifier according to Lokhande et al. (2020).
Method. Our strategy is inspired by Morcos et al. (2018) which showed that canonical correlations can reveal the similarity in neural networks: when two networks (same architecture) are trained using different labels/schemes for example, canonical correlations can indicate how similar their features are. Our observation is the following- Consider a classifier that is trained on gender (the protected attribute). and another classifier that is trained on attractiveness, if the features extracted by the latter model share a high similarity with the one trained to predict gender, then it is more likely that the latter model is influenced by features in the image pertinent to gender, which will lead to an unfairly biased trained model. We show that by imposing a loss on the canonical correlation between the network being trained (but we lack per-sample protected attribute information) and well trained classifier pre-trained on the protected attributes, we can get a more fair model. This may enable training fairer models in settings which would otherwise be difficult.
8

---

Under review as a conference paper at ICLR 2021
Implementation details. To simulate the case where only have a pretrained network protected attributes, we train a Resnet- He et al. (2016) on gender attribute, and when we train the classifier to predict attractiveness, we add loss using the canonical correlations between these two networks entropy TCcA where the first term is the standard cross entropy
on intermediate layers: Taoual = entoppy term and the second term is the canonical correlation. See A. for more details of training/evaluation.
Results. We choose two commonly used error metrics for fairness: difference in Equality of Opportunity Hardt et al. (2016) (DEO), and difference in Demographic Parity Yao & Huang (2017) (DDP), See appendix A.6 for more detailed explaination of the two metrics. We conduct experiments by applying the canonical correlation loss on three different layers in Resnet-18. In Table we can see that applying canonical correlation loss generally improves the DEO and DDP metrics (lower is better) over the standard model (trained using cross entropy loss only). Specifically, applying the loss on early layers like convi and convl gets better performance than applying at a relatively late layer like conv2. Another promising aspect of our approach is that is can easily handle the case where the protected attribute is continuous variable (as long as a well trained regression network on the protected attribute is given) while other methods like Lokhande et al. (2020); Zhang et al. (2018) need to first discretize the variable and then enforce constraints which can be much more involved.
4 RELATED WORK
Stochastic CCA: There has been much interest in designing scalable and provable algorithms for CCA: Ma et al. (2015) proposed the first stochastic algorithm for CCA, while only local convergence
proven for non-stochastic version. Wang et al. (2016) designed algorithm which uses alternating SVRG combined with shift-and-inee and-inverrt crndconiii,,,,ng, with global convergence. These stochastic methods, together with Ge et al. (016) Allen-Zhu Li (016), which reduce CCA problem to generalized eigenvalue problem and solve it by performing efficient power method, all belongs to the methods that try to solve empirical CCA problem, it can be seen as ERM approxiamtion of the priginal population objective, which requires solving numerical optimization of the empirical CCA objective on a fixed data set. These methods usually assume the access to the full dataset in the beginning, which is not very suitable for many practical applications where data tend to come in streaming way. Recently, there are increasingly interest in considering population CCA problem Arora et al. (2177) Gao et al. 2019). The main difficulty in population setting is we have limited knowledge about the objective unless we know the distribution of X and Y. Arora et al. (2017) handles this problem by deriving an estimation of gradient of population objecitve whose error can be properly bounded so that applying proximal gradient to a convex relexed objective will provably converge. Gao al. (019) provides tightened analysis of the time complexity of the algorithm in Wang et al. (2016), and provides sample complexity under certain distribution. The problem we are trying to solve in this work is the same as that in Arooa al. (2017); Gao et al. (2019): to optimize the population objective of CCA in streaming fashion.
Riemannian Optimization: Riemannian optimization is a generalization of standard Euclidean optimization methods to smooth manifolds, which takes the following form: Given f .M- - R. solve min (z),, where M is a Riemannian manifold. One advantage is that it provides a nice way to express many constrained optimization problems as unconstrained problems. Applications include matrix and tensor factorization [ssteva et al. (201T). Tan et al. (2014), PCA Edelman et al. (1998). CCA Yger al.) (2012), and so on. Y Ygrr& et al. (2012) rewrites CCA formulation as Riemannian optimization on Stiefel manifold. In our work, we further explore the ability of Riemannian optimization framework, decomposing the linear space spanned by canonical vectors into products of several matrices which lie in several different Riemannian manifolds.
5 CONCLUSIONS
In this work, we presented a stochastic approach (RSG+) for the CCA model based on the observation that the solution of CCA can be decomposed into product of matrices which lie on certain structured spaces. This affords specialized numerical schemes and makes the optimization more efficient. The optimization is based on Riemannian stochastic gradient descent and we provide a proof for its O convergence rate with kd time complexity per iteration. In experimental evaluations, we find that our RSG+ behaves favorably relative to the baseline stochastic CCA method in capturing the correlation in the datasets. We also show the use of RSG+ in the DeepCCA setting showing feasibility when scaling to large dimensions well as in an interesting use case in training fair models.

---

Under review as a conference paper at ICLR 2021
REFERENCES
P-A Absil, Robert Mahony, and Rodolphe Sepulchre. Riemannian geometry of grassmann manifolds with a view on algorithmic computation. Acta Applicandae Mathematica, 80(2):19--220, 2004.
Pierre-Antoine Absil, Robert E. Mahony, and Rodolphe Sepulchre. Optimization algorithms on matrix manifolds. 2007.
Zeyuan Allen Zhu and Yuanzhi Li. Doubly accelerated methods for faster cca and generalized cigendecomposition. In JCML, 2016.
Galen Andrew, Raman Arora, Jeff Bilmes, and Karen Livescu. Deep canonical correlation analysis. In Interaational conference on machine learning, pp. 1247-1255, 2013.
Raman Arora, Teodor Vanislavov Marinov, Poorya Mianjy, and Nati Srebro. Stochastic approximation for canonical correlation analysis. In Advances Neural Information Processing Systems, pp- 1775-47844 2017.
Gary Bécigneul and Octavian- -ugen Ganea. Riemannian adaptive optimization methods. arXiv preprint ArXiv: 1810.00760, 2018.
Kush Bhatia, Aldo Pacchiano, Nicolas Flammarion, Peter L Bartlett, and Michael I Jordan, Gen-oja: Simple & efficient algorithm for streaming generalized eigenvector computation. Advances in Neural Information Processing Systems, pp. 7016-7025, 2018.
Silvere Bonnabel. Stochastic gradient descent riemannian manifolds. IEEE Transactions Automatic Control, 58(9):2117-2229, 2013.
William M Boothby. An introduction differentiable manifolds and Riemannian geometry. Academic press, 1986.
Nicolas Boumal, Bamdev Mishra, P-A Absil, and Rodolphe Sepulchre. Manopt, matlab toolbox for optimization manifolds. The Journal of Machine Learning Research, 15(1):1455-145) 2014.
Rudrasis
Chakraborty, Liu Yang, Soren Hauberg, and Baba Vemuri. Intrinsic grassmann averages for online linear, robust and nonlinear subspace learning. IEEE Transactions Pullern Analysis and Machine Intelligence, 2020.
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 learming, pp. 129-1 36. 2009.
Heather D. Couture, Roland Kwitt, J. S. Marron, Melissa A. Troester, Charles M. Perou, and Marc Niethammer. Deep multi-view view learning via task-optimal cca. ArXiv, abs//997.007799 2019.
Michele Donini, Luca Oneto, Shai Ben-David, John S Shawe-Taylor, and Massimiliano Pontil. Empirical risk minimization under fairness constraints. In Advances in Neural Information Processing Systems, pp. 2791-2801 2018.
Alan Edelman, Tomás A. Arias, and Steven Thomas Smith. The geometry of algorithms with orthogonality constraints. SIAM J. Matrix Anal. Appl., 20:303-353, 1998.
Chao Gao, Dan Garber, Nathan Srebro, Jialei Wang, and Weiran Wang. Stochastic canonical correlation analysis. Journal of Machine Learning Research, 20(167):1 2019.
Rong Ge, Chi Jin, Pranecth Netrapalli, Aaron Sidford, et al. Efficient algorithms for large-scale gen- eralized eigenvector computation and canonical correlation analysis. In International Conference on Machine Learning, pp. 2751-2750, 2016.
Gene H Golub and Christian Reinsch. Singular value decomposition and least squares solutions. In Linear Algebra, pp. 134-151. Springer, 1971.
Gene H. Golub and Hongyuan Zha. The canonical correlations of matrix pairs and their numerical computation. 1992.
10

---

Under review as a conference paper at ICLR 2021
Moritz Hardt, Eric Price, and Nati Srebro. Equality of opportunity in supervised learning. In Advances in neural information processing systems, pp. 3315-3323, 2016.
Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun. Deep residual learning for image recognition. In Proceedings of the IEEE conference on computer vision and pattern recognition, pp. 770-778, 2016.
Sigurdur Helgason. Differential geometry and symmetric spaces, volume 341. American Mathemati- cal Soc., 2001.
Mariya Ishteva, Pierre-Antoine Absil, Sabine Van Huffel, and Lieven De Lathauwer. Best low multilinear rank approximation of higher-order tensors, based on the riemannian trust-region scheme. SIAM J. Matrix Anal. Appl., 32:115-135, 2011.
Tetsuya Kaneko, Simone Fiori, and Toshihisa Tanaka. Empirical arithmetic averaging over the compact stiefel manifold. IEEE Transactions on Signal Processing, 61(4):883-894 2012.
Alex Krizhevsky. Learning multiple layers of features from tiny images. Technical report, 2009,
Yann LeCun, Corinna Cortes, and CJ Burges. Mnist handwritten digit database. ATT Labs [Online]. 2, 2010.
Available: http://ww..ccn..o//eexi/iiiee
Jun Li, Li Fuxin, and Sinisa Todorovic. Efficient riemannian optimization on the stiefel manifold via the cayley transform. arXiv preprint arXiv:2002.01113, 2020.
Vishnu Suresh Lokhande, Aditya Kumar Akash, Sathya N Ravi, and Vikas Singh. Fairalm: Augmented lagrangian method for training fair models with little regret. arXiv preprint arXiv:2004.01355, 2020.
Yong Luo, Dacheng Tao, Kotagiri Ramamohanarao, Chao Xu, and Yonggang Wen. Tensor canonical correlation analysis for multi-view dimension reduction. IEEE transactions on Knowledge and Data Engineering, 27(11):31 (111)33111-31244 2015.
Zhuang Ma, Yichao Lu, and Dean Foster. Finding linear structure in large datasets with scalable canonical correlation analysis. In ICML, 2015.
Ari Morcos, Maithra Raghu, and Samy Bengio. Insights on representational similarity in neural networks with canonical correlation. In Advances in Neural Information Processing Systems, pp. 5727-5736, 2018
Arkadi Nemirovski, Anatoli B. Juditsky, Guanghui Lan, and Alexander Shapiro. Robust stochastic approximation approach to stochastic programming. SIAM J. Optimization, 19:1574-1609, 2009.
Erkki Oja. Simplified neuron model as a principal component analyzer. Journal of Mathematical Biology, 15:267-273, 1982.
Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, et al. Pytorch: An imperative style, pigh-performance deep learning library. In Advances Neural Information Processing Systems, pp. 0024-8035, 2019.
W Nicholson Price and Glenn Cohen. Privacy the age of medical big data. Nature medicine, 25 (1):37-43, 2019.
Markus ReiB. Martin Wahl, al. Nonasymptotic upper bounds for the reconstruction error pea. Annals of Statistics, 1098)
1098-09981111 1123,, 2020.
Danilo Jimenez Rezende, George Papamakarios, Sébastien Racanière, Michael S Albergo, Gurtej Kanwar, Phiala E Shanaban, and Kyle Cranmer. Normalizing flows on tori and spheres. arXiv preprint urXiv:2002.04428 2020.
Jan Rupnik and John Shawe Taylor. Multi-view canonical correlation analysis. In Conference on Data Mining and Data Warehouses (SiKDD 20/0), pp. 1-4, 2010.
11

---

Under review as a conference paper at ICLR 2021
Cees GM Snoek, Marcel Worring, Jan C Van Gemert, Jan-Mark Geusebroek, and Amoll WM Smeulders. The challenge problem for automated detection of 101 semantic concepts in multimedia. In Proceedings of the 14th ACM international conference on Multimedia, pp. 421-4300 2006.
Raghav Subbarao and Peter Meer. Nonlinear mean shift over riemannian manifolds. International journal of computer vision, 84(1):1, 2009.
Mingkui Tan, Ivor Wai-Hung Tsang, Li Wang. Bart Vandereycken, and Sinno Jialin Pan. Riemannian pursuit for big matrix recovery. In /CML, 2014.
Roman Vershynin. Four lectures on probabilistic methods for data science. The Mathematics of Data, IAS/Park City Mathematics Series, pp.
Weiran Wang, Raman Arora, Karen Livescu, and Jeff Bilmes. On deep multi-viww representation learning. In International Conference on Machine Learning, 1083-1092, 2015a.
Weiran Wang. Raman Arora, Karen Livescu, and Nathan Srebro. Stochastic optimization for deep cca nonlinear orthogonal iterations. In 2015 53rd Annual Allerton Conference on Communicution, Control, and Computing (Allerton), pp. 688-695. IEEE, 2015b.
Weiran Wang. Jialei Wang, Dan Garber, and Nati Srebro. Efficient globally convergent stochastic optimization for canonical correlation analysis. Advances in Neural Information Processing Systems, pp. 766-774, 2016.
Sirui Yao and Bert Huang. Beyond parity: Fairness objectives for collaborative filtering. In Advances in Neural Information Processing Systems, pp. 2921-2930 2017.
Florian Yger, Maxime Berar, Gilles Gasso, and Alain Rakotomamonjy. Adaptive canonical correlation analysis based on matrix manifolds. In ICML, 2012.
Brian Hu Zhang, Blake Lemoine, and Margaret Mitchell. Mitigating unwanted biases with adversarial learning. In Proceedings of the 2018 AAAI/ACM Conference AI, Ethics, and Society, pp. 335- 340, 2018
12

---

Under review as a conference paper at ICLR 2021
A APPENDIX
A.1 A BRIEF REVIEW OF RELEVANT DIFFERENTIAL GEOMETRY CONCEPTS
To make the paper self-contained, we briefly review certain differential geometry concepts. We only include condensed description needed for our algorithm and analysis and refer the interested reader to |Bothby (196)) for a comprehensive and rigorous treatment of the topic.
Riemannian Manifold: A Riemannian manifold, M. (of dimension m)) is defined as a (smooth) topo- logical space which is locally diffeomorphic the Euclidean space R Additionally, M is equipped with a Riemannian metric g which can be defined as 9x M x M - R, where TxM M is the tangent space at X of M, see Fig. 2
= Exp,
Y Exp,it!)
M
Figure 2: Schematic description of an exemplar manifold (M) and the visual illustration of Exp and Exp map.
If X E EM, the Riemannian Exponential map at x, denoted by Expx M + .M is defined
as 711) where
[0 + M. We can find by solving the following differential equation:
(0) = x, (vto E 1))
d7 dt
= U. In general Expx is not invertible but the inverse Expx
t=to u c M - Tx is defined only uf (X),, where is called the injectivity radius Boothby (186) of M. This concept will be useful to define the mechanics of gradient descent on the manifold.
In reformulation, we will shortly make decomposing U and V into a product of several matrices. (a) St(p, n): the manifold consists of n x column orthonormal matrices (b) Gr(p, n): the manifold consists of p-dimensional subspaces in R (c) (nn.). ). the manifold/group consists of n x n special orthogonal matrices, i.e., space of orthogonal matrices with determinant 1.
of the following manifolds, specifically, when
Differential Geometry of SO(n): SO(n) is a compact Riemannian manifold, hence by the Hopf- Rinow thcorem, it is also geodesically complete manifold Helgason 2001). Its geometry is well
understood and we recall a few relevant concepts here and refer the reader to Helgason (000T) for details. (0(n) has a Lie group structure and the corresponding Lie algebra, so(n), defined as, so(n) = {W E R I = W). In other words, so(n) (the set of Left invariant vector fields with associated Lie bracket) is the set of n x n anti-symmetric matrices. The Lie bracket, [.]. operator on so sn) is defined as the commutator, i.e., for U, V E s0(1). [ UV - VU. Now, we
can define a Riemannian metric on SO(n) as follows: (U, = trace (UTV) where, U,V E
Tx (S)()))), X E (O(n). Note that, it can be shown that this is a bi-invariant Riemannian metric. Under this bi-invariant metric, now we define the Riemannian exponential and inverse exponential map as follows. Let, X. Y E SO(n), U E Tx (SO(n)). Then, Expx (Y)) = X log(X and Expx (U) = X exp( UU where, exp, log are the matrix exponential and logarithm respectively.
Differential Geometry of the Stiefel manifold: The set of all full column rank (n p) dimensional real matrices form a Stiefel manifold, St(p, n), where n 2 p. A compact Stiefel manifold is the of all column orthonormal real matrices. When p < n. St(p,n) n) can be identified with SO(n)/SO(5 p).Note that, when we consider the quotient space, SO(n) (S)(n)S((nnn)) p), we assume that SO(n p) FF(SO(n (p) is a subgroup of SO(n), where, F SO(n p) - SO(n) defined 4
by X
is an isomorphism from SO(n p) to F(SO(n p)).
X
Differential Geometry of the Grassmannian Gr(p, n): The Grassmann manifold (or the Grass- mannian) is defined as the set of all p-dimensional linear subspaces in R and denoted by Gr(p, n). where p E Z+, E Z+ 2 p. Grassmannian is symmetric space and can be identified with the quotient space (0(n)/5 /S)//(( (O(p) x O(n p)). where (O(p) x O(n p)) is the set of all n x n matrices whose top left p X and bottom right n p
p submatrices are orthogonal and all other entries are 0, and overall the determinant is 1. A point X E Gr(p. n) can be specified by a basis, X. We say that X Col(X) if X is a basis of X, where Col(.) is the column span operator. It is easy to see that the general linear group GL(p) acts isometrically, freely and properly on St(p, n). Moreover, Gr(p, n.) can be identified with the quotient space St(p, (p) Hence, the projection
13

---

Under review as a conference paper at ICLR 2021
map II St(p, n) + Gr(p, n) is a Riemannian submersion, where [I( _ Col(X). Moreover, the triplet (St(p, n), II, Gr(p, n)) is a fiber bundle.
At every point X E St(p, n.), we can define the vertical space, Vx StSp, n) to be Ker(Il,x) ) Further, given gSt we define the horizontal space, Hxx to be the -ottotogoa complement of Vx. Now, from the theory of principal bundles, for every vector field U on Gr(p, r)), define the horizontal lift of U to be the unique vector field U on St(p, n) for which Ux Hxx and TxxXx = Unxx), for all X E St(p, n). As, IT is a Riemannian submersion, the isomorphism 1Lx|Hx Hx Gr(p, n) is an isometry from (Hx,9x) to (Tn(x))r(p)) n),, xn(())))))))))))))))))) So, (i(x) is defined as:
si (U Vu(x)) - (Ux, Vx) =
(4)
where, U.V E Tnxx n) and xx, Unxx) Vx Vi(xy, Ux Hx and Vx E Hxx
We covered the exponential map and the Riemannian metric above, and their explicit formulation for manifolds listed above is provided for easy reference Table 4
9x (U. V)
Exp x (U)
Expx (Y) -)- X0--XX(-)-- X)
St(p. x) Kaneko et al. 2012
trace (U
UV
ÜSVT svd(X + U))
Gr m) Absil et al.
2004
trace
II. (U)"(U)I,"")))
DVT,
(( -1 -8,
ÜSVT svd X U))
X=
n Subbarao Meer
2009
Trace (XUXXIV
Xexpm (TU)
Xlogm (X*))
Table Explicit forms for some operations we need. IT() returns X's column space; IT, is IT's diffcrential.
A.2 PROOF OF THEOREM 1
We first restate the assumptions from section 3.1:
Assumptions: (a) The random variables X ~ (O x.,) and Y ~ NN(,, E.y) with E, I cla and M, cla for some > 0. (b) The samples X and Y drawn from X' and y respectively have zero mean. (c) For given & < d, Mg and M, have non-zero top-k eigen values.
Let F be the trace value solution for Eq. (2), and F be the trace value solution for Eqs. equation 3a equation 3b) we next restate Theorem and give its proof:
Theorem. Under the assumptions and notations above, the approximation error E = ||
Fl is
bounded and goes to zero while the whitening constraints in equation 3D 3 are satisfied.
Proof. Let Qu, Sa, Q., S, be the solutions for Eqs. equation 3a and equation 3b Ü. V be matrices consisting of top-k eigen vectors of ('/N)XTX X and (1/N)Y respectively, V. be solutions for equation 2 Let X, XUQ,S, and Y, VVQ.S,.. Also let X,, = XU and Y, YV. Observe that mean of X,, Y,, Xx and Y, are zero. Moreover the sample covariance of X,, and Y, are given by UT CxU and V CyC respectively. Thus by the constraint in equation 2 x = I and Y. Y. Tr. Let these covariance matrices be denoted by C(X..) and C(V.) respectively. Analogously the sample covariance of X, and Y, are given by ST Q UT respectively. Let these covariance matrices be denoted by C(X,) and C(Y,) respectively.
UTTTTTCUUQ CxUQuS, and ST Cy
Using Def. know Xu, X,, Xx, and Y. follow uub-Gaussian distributions.
Let F trace (UTCxrV) which can be rewritten as trace (x Moreover, let F
trace ST) ST
which similarly can be rewritten by F trace [TTY.)
Consider the approximation error between the objective functions as E F Fl. We can rewrite E = trace (X
trace 1 y J. Due to von Neumann's 's trace inequality and Cauchy-Schwarl
14

---

Under review as a conference paper at ICLR 2021
inequality, have
E = trace
[2 - -x.1, X²Yo)
trace
(Xu (X/- (xu)xu)
Xu) /
(P.-K.))
(using Von Neumann's 's trace inequality)
Σ 0 Xx)o, -V,
1 Xx -xa)
||l|| (V, - -Y.) ||
(using Cauchy-Sewarte 's inequality)
(A.1)
where (A) denote the i singular value of matrix A and Now, using Proposition I we get
|| denotes the Frobenius norm.
2
2
113
Xa - Xx)) 11e 5 min
[222/ Ag 12,
X
A
1++1
2 , 112
(Y, (Y.-)x.) 111 min
|22/| 1ylllyylla
(A.2)
2
x
where, A = C(X,) - c(X,), - C(Y)). C(X.) and C(Y,) respectively. Now, assume that C(X,) = It and C(Y,) = Ir as X,, and Y, are
Here X and X are the eigen values of
solutions of Eq. equation 2 Furthermore assume X K+1 2 A and X A 2 A for some O. Then, we can rewrite Eq. equation A.I as
2 - C(xx)l A
2 C( min V2 C((Y))l..
E < min [V2k c(x.)ll,,
A
0. Observe that, the limiting conditions for C(X,) and C(Y)) can be satisfied by the "whitening" constraint. In other words, as C(X,) = 1g and
And C(X.) + Is or C(Y,) In. E +
C(Y.) = I,, C(X.) and C(Y,) converge to C(X,) and C(Y,), the approximation error goes to zero.
A.3 RSG+ ALGORITHM
Here we show our algorithm with more details about the gradients in every step in Alg
A.4 IMPLEMENTATION DETAILS OF CCA ON FIXED DATASET
Implementation details. On all three benchmark datasets, we only passed the data once for both RSG+ and MSG Arora et al. (2017) and MSG results. We conducted experiments on different dimensions of target space: A = 1, 2, 4. The choice of & is motivated by the fact that the spectrum of the datasets decays quickly. Since our RSG+ processes data in small blocks, we let data come in mini-batches (mini-batch size was set to 100).
the code from Arora al. (2017) produce
A.5 RUNTIME OF RSG+ AND BASELINE METHODS
The runtime comparison of RSG+ and MSG reported Table 5 Our algorithm is 5-10 times faster.
We also plot the runtime of our algorithm under different data dimension (set d, d, d) and number of total samples sampled from joint gaussian distribution in Fig. 3
A.6 ERROR METRICS FOR FAIRNESS
Equality of Opportunity (EO) Hardt et al. (2016): A classifier h is said to satisfy EO the prediction is independent of the protected attribute (in our experiment is binary variable where
15

---

Under review as a conference paper at ICLR 2021
Algorithm 2: Riemannian SGD based algorithm (RSG+) to compute canonical directions
Input: E R Nxd, E R Nxdy, &
Output: U Rd xk V E R
Initialize U, V, Qu, Su:
Partition data X, Y into batches size B. Let batch be denoted by X, and Y
for E {1, _ }
do
Gradient for top-k principal vectors: calculating Vi.Ppi Vp
1. Partition X, (Y3) into (I. = =L #. I blocks of size d. x (dy x k);
2. Let the 1th block be denoted by Zi (2?"):
orthogonalized block denoted by 2 (2,"): 4. the subspace spanned by each 2i (and 2 ) 2, Gr(k; da) (and 2 E Gr(k; dg)):
3. Orthogonalize each block and
ViFm
Σ Exp (i))
ViFri Σ Exp (zi)
Gradient from equation
calculating Vij Viram. Vq, Fem, Vq. Fems Vss Feam Vs, Fen
of Vi Va Feam
of a0
V F
1 oF
V:
a0
0V
VQ, Funm aQ aQu 05..
V Vq, F
a²
oF
aQu
aQ
S.,,, Fuan Upper
s. Fem = Upper
aF
85.
Here. Upper returns the upper triangular matrix of the input matrix and oF OF oF oF OF oF give Euclidean gradients. completeness, the closed form expression the gradients is,
aU 5V. UQ_ UQ. U5_ US,
oF aU Cx.VQ..S aF av
aF U" Cxy VQ.S... UQ
aF QU''CCr
as
aF V CrxUÚQ.S...
aF STQ Sxx.v'"C.v
Cxr
aQ.
Os.
(6)
Gradient update canonical directions
V = V Fga + Fy
VpF= = V + Ve
Vx F Vx Fun where, X is generic entity: E (Qu, Qu, Su, Sc):
Batch update of canonical directions
Exp VV where, generic entity:
{Ü, V. Q.,, Qu, S.)
end
UQ.S and VQ....:
Runtime
800
1-0000
N- N-10500
700
N-20000
Time one paos
800
500
400
800
200
2000
4000
8000
8000
10000
Data dimension d
Figure 3: Runtime of RSG+ under different data dimensions and size of datasets.
16

---

Under review as a conference paper at ICLR 2021
Table 5: Wallclock runtime of one pass through the data of our RSG+ and MSG on MNIST, Mediamill and CIFAR (average of 5 runs).
MNIST
Mediamill
CIFAR k
Time (s)
k=
k
2
k
k
h ==
RSG+ (Ours)
4.16
4.24
4.71
1,89
1.60
1.44
14.80
17.22
22.10
MSG
35.32
42.09
49.17
11.59
14.21
17.34
80.21
100.80
106.55
1 stands for Male and 0 stands for Female) for classification label {0, 1}. We use the difference of false negative rate (conditioned on y = 1) across two groups identified by protected attribute the error metric, and we denote it DEO.
Demographic Parity (DP) Yao Huang (017): A classifier h satisfies DP the likelihodd of making misclassification among the positive predictions of the classifier is independent of the protected attribute 8. We denote the difference of demographic parity between two groups identified by the protected attribute as DDP.
A.7 IMPLEMENTATION DETAILS FAIRNESS EXPERIMENTS
Implementation details. The network is trained for 20 epochs with learning rate 0.01 and batch size 256. We follow Donini et al. (2018) to use NVP (novel validation procedure) to evaluate our result: first search for hyperparameters that achieves the highest classification score and then report the performance the model which gets minimum fairness error metrics with accuracy within the highest 90% accuracies. When apply RSG+ certain layers, first randomized projection to project the feature into 1 dimension, and then extract top-10 canonical components for training. Similar to our previous experiments DeepCCA, the batch method does scale to lk dimension.
A.8 RESNET- 18 ARCHITECTURE AND POSITION OF CoNv-0, IN TABLE 3
The Resnet- 18 contains first convolutional layer followed by normalization, nonlinear activation, and max pooling. Then it has four residual blocks, followed by average polling and fully connected layer. We denote the position after the first convolutional layer as convO, the position after the first residual block as convl and the position after the second residual block conv2. We choose early layers since late layers close to the final fully connected layer will have feature that is more directly relevant the classification variable (attracfiveness in this case).
17


Reviews:
[{"review_id": "yMX_CbHRtO7", "reviewer": "AnonReviewer1", "summary": "", "questions": "", "limitations": "", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "main contribution\n\nThis work offers a stochastic CCA algorithm based on Riemannian optimization approach\n\nStrength\n\n- The paper offers a theory-backed algorithm for CCA under the assumption that the two views are sub-Gaussian. The complexity-accuracy tradeoff of the algorithm seems to be appealing according to the experiments.\n\nWeakness\n\nOverall, the biggest concern is readability. The paper is very packed and many treatments seem to be unbalanced. Important details are missing, and proofs seem to be hastily. The paper may need some re-packaging and re-organization before its core technical contents could be easily followed.\n\n- readability. The biggest concern of the work is that it is very hard to read. This creates a lot of barriers in understanding key aspects of the paper, e.g., contributions, formulations, novelty of the proof, just to name a few. The key formulation and the definitions of the manifolds were not clearly defined. The equivalence of (1) and (3) was not clearly shown. The theorems were presented in a bit abrupt way. Even the algorithm does not appear in the maintext but the appendix (which is also hard to read). The color code is used in a way that is a bit confusing (does ICLR allow writing in different colors?).\n\n- Clarify about the contribution. It is hard to clearly see how much is the contribution. The proofs seem to be short and most of the proofs are presented using “propositions” cited from existing papers. If the authors think the contributions lie in reformulating the CCA problem as a PCA problem, then the reformulating part is perhaps the contribution. But from the current writing it is hard to follow how the reformulation comes through and how the reformulation enables using these existing propositions to prove convergence of the proposed algorithm.\n\n- It is unclear why the reformulation in (3) is a good approximation for CCA. Some attempts for justifying this were offered in Theorem 1, but it was based on the assumption X and Y are sub-Gaussian, which is at best a special case, even if the proof is correct (which, due to the current organization of the paper, is hard to read and fully understand).\n\n- “SO(n): group/manifold of nxn special orthogonal matrices.” what is the definition of “special” here? are they the commonly understood orthogonal matrices?\n\n- The paper has this upper triangular structure of the S matrices but this point seems to have no detailed explanation. If one understands U = \\tilde{U}S as the QR decomposition, then indeed S is upper triangular, but why is there a Q in the middle? Why is this Q useful?\n\n- the “high-level” description of the algorithm says the idea is to connect CCA with PCA using this reformulation, but this point was not clearly explained.\n\n- Complexity. From the current writing, it is very hard to see how the complexity of the algorithm is calculated. The gradients needed are tabulated in the supplementary materials, but it is very hard for a reader to directly see why the algorithm saves computational complexity.\n\n- It is also unclear how the convergence and convergence rate analyses come together. These parts may need to be elaborated.\n", "extra_scores": {}, "extra_text": {}}, {"review_id": "wfYblg_7hRr", "reviewer": "AnonReviewer3", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "This paper aims to reduce the computational complexity of canonical correlation analysis.By decomposing the CCA projection matrices into a product of several structured matrices, a stochastic gradient-based optimization on a Riemannian manifold is provided reducing the computational complexity from $d^3$ to $d^2k$.\n\n\nStrength:\nCCA is a classic and still important method, especially in combination with deep neural networks (e.g., multi-view learnings). \nThe proposed method enables the applications of CCA to high-dimensional vectors with small memory.\nThis makes it easier to use CCA as an objective function of deep neural networks, which is trained on GPUs.\nExperiments show the benefits of their proposed method, in particular, the computational speed is 5-10 times faster than the existing method (MSG).\n\nWeakness:\nAlthough the authors claim that their proposed method captures more correlation than MSG, two of the three datasets in which their method is superior (MNIST and CIFAR) are not realistic setting (i.e., correlation between left/right half of images).\n\nQuestion:\nIs it possible to make a (numerical) comparison with Yger+2012, which reformulates CCA as an optimization on the generalized Stiefel manifold?", "extra_scores": {}, "extra_text": {}}, {"review_id": "bUrYiJNWEff", "reviewer": "AnonReviewer2", "summary": "", "questions": "", "limitations": "", "rating": 7, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "The paper presents an approach to find canonical directions in a streaming fashion, i.e. without direct calculation of covariance matrices (which becomes hard when the number of examples is large). This solution to that task is not obvious, because the objective function of CCA, together with whitening constraints, does not allow simple additive decomposition.\n\nFirst, the optimization task is reformulated as a task over certain Riemannian manifolds, and a natural initialization is suggested. It is shown that under a certain assumption, this initialization is already a solution of good quality. Then, a natural minimization algorithm is presented, which is based on stochastic gradient descent on a Riemannian manifold. The key aspect of the algorithm is a combination of 2 types of gradient, the gradient for top-k principal vectors and standard gradient.\n\nThe experimental part shows that the algorithm successfully solves CCA in a streaming fashion. Also, it can be effectively combined with deep feature learning (Andrew, 2013) to find common features for multi-view representation learning tasks. Experiments look convincing.", "extra_scores": {}, "extra_text": {}}, {"review_id": "D_YudmtLrJr", "reviewer": "AnonReviewer4", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "1. Paper summary:\n\nThis paper proposes a method for solving linear CCA on high dimensional data. Linear CCA has a closed form solution. The solution requires a whitening step that costs O(d^3). This makes it not applicable to data in high dimensional spaces, e.g. representations learnt by deep networks. \n\nTo resolve the issue, the authors propose a reformulation of linear CCA which decomposes the transformation matrix U (V) into a product of three matrices. Those three matrices have the following properties:\n\n- Their initial values can be obtained by efficient streaming PCA on original view matrix X (Y). Streaming PCA costs O(d^2 * k) only where k is the top k eigenvectors.\n\n- They allow for Riemannian stochastic gradient descent which ensures their updated values lie on the same manifold.\n\n2. Strong points of the paper:\n\nThe new linear CCA formulation justifies the rationale of batch CCA training.\n\nUnder Gaussian distribution assumption:\n- The absolute difference between correlation found by original linear CCA and stochastic one is bounded.\n- The convergence of the training process is proven.\n\nThe experiments are performed on different aspects:\n- Recovering groundtruth transformations on MNIST, CIFAR and Mediamill data sets.\n- Learning deep features on MNIST data set.\n- Improving fairness in deep learning by adding CCA term to the loss function.\n\n3. Weak points of the paper:\n\nThe theoretical results are obtained under a strong assumption that X and Y both have Gaussian distribution.\n\nMost propositions are from other papers.", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: This paper gives a new algorithm for the CCA problem. The main idea of the new algorithm is to reformulate the matrices in the CCA problem as a product of three matrices: one orthonormal matrix, one rotation and one upper-diagonal matrix. The algorithm then performs remannian gradient descent to these components. The per-iteration complexity of the algorithm is O(d^2k) while the (local) convergence rate is O(1/t). Overall the reformulation is interesting and the algorithm seems effective in practice. On the other hand the convergence rate proof relies on local strong convexity and it's not clear why the algorithm converges globally (or even what is the radius of convergence locally).
Decision: Reject