Anastasia Podosinnikova

Francis Bach

Simon Lacoste-Julien

Beyond CCA: Moment Matching for Multi-View Models

We introduce three novel semi-parametric extensions of probabilistic canonical correlation analysis with identifiability guarantees. We consider moment matching techniques for estimation in these models. For that, by drawing explicit links between the new models and a discrete version of independent component analysis (DICA), we first extend the DICA cumulant tensors to the new discrete version of CCA. By further using a close connection with independent component analysis, we introduce generalized covariance matrices, which can replace the cumulant tensors in the moment matching framework, and, therefore, improve sample complexity and simplify derivations and algorithms significantly. As the tensor power method or orthogonal joint diagonalization are not applicable in the new setting, we use non-orthogonal joint diagonalization techniques for matching the cumulants. We demonstrate performance of the proposed models and estimation techniques on experiments with both synthetic and real datasets.

Introduction

Canonical correlation analysis (CCA), originally introduced by [START_REF] Hotelling | Relations between two sets of variates[END_REF], is a common statistical tool for the analysis of multi-view data. Examples of such data include, for instance, representation of some text in two languages (e.g., Vinokourov et al., 2002) or images aligned with text data (e.g., [START_REF] Hardoon | Canonical correlation analysis: An overview with application to learning methods[END_REF][START_REF] Gong | A multi-view embedding space for modeling internet images, tags, and their semantics[END_REF]. Given two multidimensional variables (or datasets), CCA finds two linear transformations (factor loading matrices) that mutually maximize the correlations between the transformed variables (or datasets). Together with its kernelized version (see, e.g., [START_REF] Cristianini | An Introduction to Support Vector Machines and Other Kernel-based Learning Methods[END_REF][START_REF] Bach | Kernel independent component analysis[END_REF], CCA has a wide range of applications (see, e.g., [START_REF] Hardoon | Canonical correlation analysis: An overview with application to learning methods[END_REF] for an overview). [START_REF] Bach | A probabilistic interpretation of canonical correlation analysis[END_REF] provide a probabilistic interpretation of CCA: they show that the maximum likelihood estimators of a particular Gaussian graphical model, which we refer to as Gaussian CCA, is equivalent to the classical CCA by [START_REF] Hotelling | Relations between two sets of variates[END_REF]. The key idea of Gaussian CCA is to allow some of the covariance in the two observed variables to be explained by a linear transformation of common independent sources, while the rest of the covariance of each view is explained by their own (unstructured) noises. Importantly, the dimension of the common sources is often significantly smaller than the dimensions of the observations and, potentially, than the dimensions of the noise. Examples of applications and extensions of Gaussian CCA are the works by [START_REF] Socher | Connecting modalities: Semi-supervised segmentation and annotation of images using unaligned text corpora[END_REF], for mapping visual and textual features to the same latent space, and [START_REF] Haghighi | Learning bilingual lexicons from monolingual corpora[END_REF], for machine translation applications.

Gaussian CCA is subject to some well-known unidentifiability issues, in the same way as the closely related factor analysis model (FA; [START_REF] Bartholomew | Latent Variable Models and Factor Analysis[END_REF][START_REF] Basilevsky | Statistical Factor Analysis and Related Methods: Theory and Applications[END_REF] and its special case, the probabilistic principal component analysis model (PPCA; [START_REF] Tipping | Probabilistic principal component analysis[END_REF][START_REF] Roweis | EM algorithms for PCA and SPCA[END_REF]. Indeed, as FA and PPCA are identifiable only up to multiplication by any orthogonal rotation matrix, Gaussian CCA is only identifiable up to multiplication by any invertible matrix. Although this unidentifiability does not affect the predictive performance of the model, it does affect the factor loading matrices and hence the interpretability of the latent factors. In FA and PPCA, one can enforce additional constraints to recover unique factor loading matrices (see, e.g., [START_REF] Murphy | Machine Learning: A Probabilistic Perspective[END_REF]. A notable identifiable version of FA is independent component analysis (ICA; [START_REF] Jutten | Calcul neuromimétique et traitement du signal: Analyse en composantes indépendantes[END_REF][START_REF] Jutten | Blind separation of sources, part I: An adaptive algorithm based on neuromimetric architecture[END_REF][START_REF] Comon | Handbook of Blind Source Separation: Independent Component Analysis and Applications[END_REF]. One of our goals is to introduce identifiable versions of CCA.

The main contributions of this paper are as follows. We first introduce for the first time, to the best of our knowledge, three new formulations of CCA: discrete, non-Gaussian, and mixed (see Section 2.1). We then provide identifiability guarantees for the new models (see Section 2.2). Then, in order to use a moment matching framework for estimation, we first derive a new set of cumulant tensors for the discrete version of CCA (Section 3.1). We further replace these tensors with their approximations by generalized covariance matrices for all three new models (Section 3.2). Finally, as opposed to standard approaches, we use a particular type of non-orthogonal joint diagonalization algorithms for extracting the model parameters from the cumulant tensors or their approximations (Section 4).

Models. The new CCA models are adapted to applications where one or both of the data-views are either counts, like in the bag-of-words representation for text, or continuous data, for instance, any continuous representation of images. A key feature of CCA compared to joint PCA is the focus on modeling the common variations of the two views, as opposed to modeling all variations (including joint and marginal ones).

Moment matching. Regarding parameter estimation, we use the method of moments, also known as "spectral methods". It recently regained popularity as an alternative to other estimation methods for graphical models, such as approximate variational inference or MCMC sampling. Estimation of a wide range of models is possible within the moment matching framework: ICA (e.g., [START_REF] Cardoso | Independent component analysis, A survey of some algebraic methods[END_REF][START_REF] Comon | Handbook of Blind Source Separation: Independent Component Analysis and Applications[END_REF], mixtures of Gaussians (e.g., [START_REF] Arora | Learning mixtures of separated nonspherical Gaussians[END_REF]Hsu & Kakade, 2013), latent Dirichlet allocation and topic models [START_REF] Arora | Learning topic models -Going beyond SVD[END_REF][START_REF] Arora | A practical algorithm for topic modeling with provable guarantees[END_REF][START_REF] Anandkumar | A spectral algorithm for latent Dirichlet allocation[END_REF][START_REF] Podosinnikova | Rethinking LDA: Moment matching for discrete ICA[END_REF], supervised topic models [START_REF] Wang | Spectral methods for supervised topic models[END_REF], Indian buffet process inference [START_REF] Tung | Spectral methods for Indian buffet process inference[END_REF], stochastic languages [START_REF] Balle | Method of moments for learning stochastic languages: Unified presentation and empirical comparison[END_REF], mixture of hidden Markov models [START_REF] Sübakan | Spectral learning of mixture of hidden Markov models[END_REF], neural networks (see, e.g., [START_REF] Anandkumar | Learning mixed membership community models in social tagging networks through tensor methods[END_REF][START_REF] Janzamin | Beating the perils of non-convexity: Guaranteed training of neural networks using tensor methods[END_REF], and other models (see, e.g., Anandkumar et al., 2014, and references therein).

Moment matching algorithms for estimation in graphical models mostly consist of two main steps: (a) construction of moments or cumulants with a particular diagonal structure and (b) joint diagonalization of the sample estimates of the moments or cumulants in order to estimate the parameters.

Cumulants and generalized covariance matrices. By using the close connection between ICA and CCA, we first derive in Section 3.1 the cumulant tensors for the discrete version of CCA from the cumulant tensors of a discrete version of ICA (DICA) proposed by [START_REF] Podosinnikova | Rethinking LDA: Moment matching for discrete ICA[END_REF]. Extending the ideas from the ICA literature [START_REF] Yeredor | Blind source separation via the second characteristic function[END_REF]Todros & Hero, 2013), we further generalize in Section 3.2 cumulants as the derivatives of the cumulant generating function. This allows us to replace cumulant tensors with "generalized covariance matrices", while preserving the rest of the framework. As a consequence of working with the second order information only, the derivations and algorithms get significantly simplified and the sample complexity potentially improves.

Non-orthogonal joint diagonalization. When estimating model parameters, both CCA cumulant tensors and generalized covariance matrices for CCA lead to non-symmetric approximate joint diagonalization problems. Therefore, the workhorses of the method of moments in similar context -orthogonal diagonalization algorithms, such as the tensor power method [START_REF] Anandkumar | Tensor decompositions for learning latent variable models[END_REF], and orthogonal joint diagonalization [START_REF] Bunse-Gerstner | Numerical methods for simultaneous diagonalization[END_REF][START_REF] Cardoso | Jacobi angles for simultaneous diagonalization[END_REF] -are not applicable. As an alternative, we use a particular type of non-orthogonal Jacobi-like joint diagonalization algorithms (see Section 4). Importantly, the joint diagonalization problem we deal with in this paper is conceptually different from the one considered, e.g., by [START_REF] Kuleshov | Tensor factorization via matrix factorization[END_REF] (and references therein) and, therefore, the respective algorithms are not applicable here.

2 Multi-view models 2.1 Extensions of Gaussian CCA Gaussian CCA. Classical CCA [START_REF] Hotelling | Relations between two sets of variates[END_REF] aims to find projections D 1 ∈ R M1×K and D 2 ∈ R M2×K , of two observation vectors x 1 ∈ R M1 and x 2 ∈ R M2 , each representing a data-view, such that the projected data, D 1 x 1 and D 2 x 2 , are maximally correlated. Similarly to classical PCA, the solution boils down to solving a generalized SVD problem. The following probabilistic interpretation of CCA is well known [START_REF] Browne | The maximum-likelihood solution in inter-battery factor analysis[END_REF][START_REF] Bach | A probabilistic interpretation of canonical correlation analysis[END_REF][START_REF] Klami | Bayesian canonical correlation analysis[END_REF]. Given that K sources are i.i.d. standard normal random variables, α ∼ N (0, I K ), the Gaussian CCA model is given by

x 1 | α, µ 1 , Ψ 1 ∼ N (D 1 α + µ 1 , Ψ 1 ), x 2 | α, µ 2 , Ψ 2 ∼ N (D 2 α + µ 2 , Ψ 2 ), (1)
where the matrices Ψ 1 ∈ R M1×M1 and Ψ 2 ∈ R M2×M2 are positive semi-definite. Then, the maximum likelihood solution of (1) coincides (up to permutation, scaling, and multiplication by any invertible matrix) with the classical CCA solution. The model ( 1) is equivalent to

x 1 = D 1 α + ε 1 , x 2 = D 2 α + ε 2 , (2)
α k ε 1 ε 2 x 1m x 2m D 1 D 2 K M 1 M 2 N
Figure 1: Graphical models for non-Gaussian (4), discrete (5), and mixed (6) CCA.

where the noise vectors are normal random variables, i.e. ε 1 ∼ N (µ 1 , Ψ 1 ) and ε 2 ∼ N (µ 2 , Ψ 2 ), and the following independence assumptions are made:

α 1 , . . . , α K are mutually independent, α ⊥ ⊥ ε 1 , ε 2 and ε 1 ⊥ ⊥ ε 2 .
(3)

The following three models are our novel semi-parametric extensions of Gaussian CCA (1)-( 2).

Multi-view models. The first new model follows by dropping the Gaussianity assumption on α, ε 1 , and ε 2 .

In particular, the non-Gaussian CCA model is defined as

x 1 = D 1 α + ε 1 , x 2 = D 2 α + ε 2 , (4)
where, as opposed to (2), no assumptions are made on the sources α and the noise ε 1 and ε 2 except for the independence assumption (3).

By analogy with [START_REF] Podosinnikova | Rethinking LDA: Moment matching for discrete ICA[END_REF], we can further "discretize" non-Gaussian CCA (4) by applying the Poisson distribution to each view (independently on each variable):

x 1 | α, ε 1 ∼ Poisson(D 1 α + ε 1 ), x 2 | α, ε 2 ∼ Poisson(D 2 α + ε 2 ). (5)
We obtain the (non-Gaussian) discrete CCA (DCCA) model, which is adapted to count data (e.g., such as word counts in the bag-of-words model of text). In this case, the sources α, the noise ε 1 and ε 2 , and the matrices D 1 and D 2 have non-negative components.

Finally, by combining non-Gaussian and discrete CCA, we also introduce the mixed CCA (MCCA) model:

x 1 = D 1 α + ε 1 , x 2 | α, ε 2 ∼ Poisson(D 2 α + ε 2 ), (6)
which is adapted to a combination of discrete and continuous data (e.g., such as images represented as continuous vectors aligned with text represented as counts). Note that no assumptions are made on distributions of the sources α except for independence (3).

The plate diagram for the models (4)-( 6) is presented in Fig. 1. Depending on the context, the matrices D 1 and D 2 are called differently: topic matrices1 in the topic learning context, factor loading or projection matrices in the FA and/or PPCA context, mixing matrices in the ICA context, or dictionaries in the dictionary learning context. In this paper, we will use the name factor loading matrices to refer to D 1 and D 2 .

Relation between PCA and CCA. The important difference between Gaussian CCA and the closely related FA/PPCA models is that the noise in each view of Gaussian CCA is not assumed to be isotropic unlike in FA/PPCA. In other words, the components of the noise are not assumed to be independent or, equivalently, the noise covariance matrix does not have to be diagonal and may exhibit a strong structure. In this paper, we never assume any diagonal structure of the covariance matrices of the noises of the models (4)-(6).

The following example illustrates the mentioned relation. Assuming a linear structure for the noise, (non-) Gaussian CCA (NCCA) takes the form

x 1 = D 1 α + F 1 β 1 , x 2 = D 2 α + F 2 β 2 , (7)
where

ε 1 = F 1 β 1 with β 1 ∈ R K1 and ε 2 = F 2 β 2 with β 2 ∈ R K2
. By stacking the vectors on the top of each other

x = x 1 x 2 , D = D 1 F 1 0 D 2 0 F 2 , z =   α β 1 β 2   , (8)
one can rewrite the model as x = Dz. If the noise sources β 1 and β 2 are assumed to have mutually independent components, ICA is recovered. If the sources z are further assumed to be Gaussian, x = Dz corresponds to PPCA. However, we do not assume the noise in Gaussian CCA (as well as in the models (4)-( 6)) to have a very specific low dimensional structure.

Related work. Some extensions of Gaussian CCA were proposed in the literature: exponential family CCA [START_REF] Virtanen | Bayesian exponential family projections[END_REF][START_REF] Klami | Bayesian exponential family projections for coupled data sources[END_REF] and Bayesian CCA (see, e.g., Klami et al., 2013, and references therein).

Although exponential family CCA can also be discretized, it assumes in practice that the prior of the sources is a specific combination of Gaussians. Bayesian CCA models the factor loading matrices and the covariance matrix of Gaussian CCA. Sampling or approximate variational inference are used for estimation and inference in both models. Both models, however, lack our identifiability guarantees and are quite different from the models (4)-( 6). [START_REF] Song | Nonparametric estimation of multi-view latent variable models[END_REF] consider a multi-view framework to deal with non-parametric mixture components, while our approach is semi-parametric with an explicit linear structure (our loading matrices) and makes the explicit link with CCA.

Identifiability

In this section, the identifiability of the factor loading matrices D 1 and D 2 is discussed. In general, for the type of models considered, the unidentifiability to permutation and scaling cannot be avoided. In practice, this unidentifiability is however easy to handle and, in the following, we only consider identifiability up to permutation and scaling.

ICA can be seen as an identifiable analog of FA/PPCA. Indeed, it is known that the mixing matrix D of ICA is identifiable if at most one source is Gaussian [START_REF] Comon | Independent component analysis, A new concept?[END_REF]. The factor loading matrix of FA/PPCA is unidentifiable since it is defined only up to multiplication by any orthogonal rotation matrix.

Similarly, the factor loading matrices of Gaussian CCA (1), which can be seen as a multi-view extension of PPCA, are identifiable only up to multiplication by any invertible matrix [START_REF] Bach | A probabilistic interpretation of canonical correlation analysis[END_REF]. We show the identifiability results for the new models (4)-( 6): The factor loading matrices of these models are identifiable if at most one source is Gaussian (see Appendix A for a proof).

Theorem 1. Assume that matrices D 1 ∈ R M1×K and D 2 ∈ R M2×K , where K ≤ min(M 1 , M 2 ), have full rank. If the covariance matrices cov(x 1 ) and cov(x 2 ) exist and if at most one source α k , for k = 1, . . . , K, is Gaussian and none of the sources are deterministic, then the models (4)-( 6) are identifiable (up to scaling and joint permutation).

Importantly, the permutation unidentifiability does not destroy the alignment in the factor loading matrices, that is, for some permutation matrix P , if D 1 P is the factor loading matrix of the first view, than D 2 P must be the factor loading matrix of the second view. This property is important for the interpretability of the factor loading matrices and, in particular, is used in the experimental Section 5.

The cumulants and generalized covariances

In this section, we first derive the cumulant tensors for the discrete CCA model (Section 3.1) and then generalized covariance matrices (Section 3.2) for the models (4)-( 6). We show that both cumulants and generalized covariances have a special diagonal form and, therefore, can be efficiently used within the moment matching framework (Section 4).

From discrete ICA to discrete CCA

In this section, we derive the DCCA cumulants as an extension of the cumulants of discrete independent component analysis (DICA; [START_REF] Podosinnikova | Rethinking LDA: Moment matching for discrete ICA[END_REF].

Discrete ICA. [START_REF] Podosinnikova | Rethinking LDA: Moment matching for discrete ICA[END_REF] consider the discrete ICA model ( 9), where x ∈ R M has conditionally independent Poisson components with mean Dα and α ∈ R K has independent non-negative components:

x | α ∼ Poisson(Dα). (9)
For estimating the factor loading matrix D, [START_REF] Podosinnikova | Rethinking LDA: Moment matching for discrete ICA[END_REF] propose an algorithm based on the moment matching method with the cumulants of the DICA model. In particular, they define the DICA Scovariance matrix and T-cumulant tensor as

S := cov(x) -diag [Ex] , [T ] m1m2m3 := cum(x) m1m2m3 + [τ ] m1m2m3 , (10)
where the indices m 1 , m 2 , and m 3 take the values in 1, . . . , M , and

[τ ] m1m2m3 = 2δ m1m2m3 Ex m1 -δ m2m3 cov(x) m1m2 -δ m1m3 cov(x) m1m2 -δ m1m2 cov(x) m1m3
with δ being the Kronecker delta. For completeness, we outline the derivation by [START_REF] Podosinnikova | Rethinking LDA: Moment matching for discrete ICA[END_REF] below. Denoting y := Dα, one obtains by the law of total expectation that E(x) = E(x|y) = E(y) and by the law of total covariance:

cov(x) = E[cov(x|y)] + cov[E(x|y), E(x|y)] = diag[E(y)] + cov(y),
since all the cumulants of a Poisson random variable with parameter y are equal to y. Therefore, S = cov(y).

Similarly, from the law of total cumulance, T = cum(y). Then, by the multilinearity property for cumulants, one obtains

S = D cov(α)D , T = cum(α) × 1 D × 2 D × 3 D , (11)
which is called the diagonal form since the covariance cov(α) and cumulant cum(α) of the independent sources are diagonal. Note that × i denotes the i-mode tensor-matrix product (see, e.g., [START_REF] Kolda | Tensor decompositions and applications[END_REF]. This diagonal form is further used for estimation of D (see Section 4).

Noisy discrete ICA. The following noisy version ( 12) of the DICA model reveals the connection between DICA and DCCA. Noisy discrete ICA is obtained by adding non-negative noise ε, such that α ⊥ ⊥ ε, to discrete ICA (9):

x | α, ε ∼ Poisson (Dα + ε) . (12)
Let y := Dα + ε and S and T are defined as in (10). Then a simple extension of the derivations from above gives S = cov(y) and T = cum(y). Since the covariance matrix (cumulant tensor) of the sum of two independent multivariate random variables, Dα and ε, is equal to the sum of the covariance matrices (cumulant tensors) of these variables, the "perturbed" version of the diagonal form (11) follows

S = Dcov(α)D + cov(ε), T = cum(α) × 1 D × 2 D × 3 D + cum(ε). (13)
DCCA cumulants. By analogy with (8), stacking the observations x = [x 1 ; x 2 ], the factor loading matrices D = [D 1 ; D 2 ], and the noise vectors ε = [ε 1 ; ε 2 ] of discrete CCA (5) gives a noisy version of discrete ICA with a particular form of the covariance matrix of the noise:

cov(ε) = cov(ε 1 ) 0 0 cov(ε 2 ) , (14)
which is due to the independence ε 1 ⊥ ⊥ ε 2 . Similarly, the cumulant cum(ε) of the noise has only two diagonal blocks which are non-zero. Therefore, considering only those parts of the S-covariance matrix and T-cumulant tensor of noisy DICA that correspond to zero blocks of the covariance cov(ε) and cumulant cum(ε) gives immediately a matrix and tensor with a diagonal structure similar to the one in (11). Those blocks are the cross-covariance and cross-cumulants of x 1 and x 2 .

We define the S-covariance matrix of discrete CCA2 as the cross-covariance matrix of x 1 and x 2 :

S 12 := cov(x 1 , x 2 ). (15)
From ( 13) and ( 14), the matrix S 12 has the following diagonal form

S 12 = D 1 cov(α)D 2 . ( 16
)
Similarly, we define the T-cumulant tensors of discrete CCA ( T 121 ∈ R M1×M2×M1 and T 122 ∈ R M1×M2×M2 ) through the cross-cumulants of x 1 and x 2 , for j = 1, 2:

[T 12j ] m1m2 mj := [cum(x 1 , x 2 , x j )] m1m2 mj -δ mj mj [cov(x 1 , x 2 )] m1m2 , (17)
where the indices m 1 , m 2 , and mj take the values m 1 ∈ 1, . . . , M 1 , m 2 ∈ 1, . . . , M 2 , and mj ∈ 1, . . . , M j . From ( 11) and the mentioned block structure (14) of cov(ε), the DCCA T-cumulants have the diagonal form:

T 121 = cum(α) × 1 D 1 × 2 D 2 × 3 D 1 , T 122 = cum(α) × 1 D 1 × 2 D 2 × 3 D 2 . (18)
In Section 4, we show how to estimate the factor loading matrices D 1 and D 2 using the diagonal form ( 16) and ( 18). Before that, in Section 3.2, we first derive the generalized covariance matrices of discrete ICA and the CCA models ( 4)-( 6) as an extension of the ideas by [START_REF] Yeredor | Blind source separation via the second characteristic function[END_REF]; Todros & Hero (2013).

Generalized covariance matrices

In this section, we introduce the generalization of the S-covariance matrix for both DICA and the CCA models ( 4)-( 6), which are obtained through the Hessian of the cumulant generating function. We show that (a) the generalized covariance matrices can be used for approximation of the T-cumulant tensors using generalized derivatives and (b) in the DICA case, these generalized covariance matrices have the diagonal form analogous to (11), and, in the CCA case, they have the diagonal form analogous to ( 16). Therefore, generalized covariance matrices can be seen as a substitute for the T-cumulant tensors in the moment matching framework. This (a) significantly simplifies derivations and the final expressions used for implementation of resulting algorithms and (b) potentially improves the sample complexity, since only the second order information is used.

Generalized covariance matrices. The idea of generalized covariance matrices3 is inspired by the similar extension of the ICA cumulants by [START_REF] Yeredor | Blind source separation via the second characteristic function[END_REF].

The cumulant generating function (CGF) of a multivariate random variable x ∈ R M is defined as

K x (t) = log E(e t x ), (19)
for t ∈ R M . The cumulants κ s (x), for s = 1, 2, 3, . . . , are the coefficients of the Taylor series expansion of the CGF evaluated at zero. Therefore, the cumulants are the derivatives of the CGF evaluated at zero:

κ s (x) = ∇ s K x (0), s = 1, 2, 3, . . . , where ∇ s K x (t)
is the s-th order derivative of K x (t) with respect to t. Thus, the expectation of x is the gradient E(x) = ∇K x (0) and the covariance of x is the Hessian cov(x) = ∇ 2 K x (0) of the CGF evaluated at zero.

The extension of cumulants then follows immediately: for t ∈ R M , we refer to the derivatives ∇ s K x (t) of the CGF as the generalized cumulants. The respective parameter t is called a processing point. In particular, the gradient, ∇K x (t), and Hessian, ∇ 2 K x (t), of the CGF are referred to as the generalized expectation and generalized covariance matrix, respectively:

E x (t) := ∇K x (t) = E(xe t x ) E(e t x ) , (20)
C x (t) := ∇ 2 K x (t) = E(xx e t x ) E(e t x ) -E x (t)E x (t) . ( 21
)
Some properties of these statistics and their natural finite sample estimators are analyzed by Slapak & Yeredor (2012b).

We now outline the key ideas of this section. When a multivariate random variable α ∈ R K has independent components, its CGF K α (h) = log E(e h α ), for some h ∈ R K , is equal to a sum of decoupled terms:

K α (h) = k log E(e h k α k ). Therefore, the Hessian ∇ 2 K α (h) of the CGF K α (h) is diagonal (see Appendix B.1
). Like covariance matrices, these Hessians (a.k.a. generalized covariance matrices) are subject to the multilinearity property for linear transformations of a vector, hence the resulting diagonal structure of the form (11). This is essentially the previous ICA work [START_REF] Yeredor | Blind source separation via the second characteristic function[END_REF]Todros & Hero, 2013). Below we generalize these ideas first to the discrete ICA case and then to the CCA models ( 4)-( 6).

Discrete ICA generalized covariance matrices. Like covariance matrices, generalized covariance matrices of a vector with independent components are diagonal: they satisfy the multilinearity property C Dα (h) = D C α (h)D , and are equal to covariance matrices when h = 0. Therefore, we can expect that the derivations of the diagonal form (11) of the S-covariance matrices extends to the generalized covariance matrices case. By analogy with (10), we define the generalized S-covariance matrix of DICA:

S(t) := C x (t) -diag[E x (t)]. (22)
To derive the analog of the diagonal form (11) for S(t), we have to compute all the expectations in ( 20) and ( 21) for a Poisson random variable x with the parameter y = Dα. Just to provide some intuition, we compute here one of these expectations (see Appendix B.2 for further derivations):

E(xx e t x ) = E[E(xx e t x | y)] = diag[e t ] E(yy e y (e t -1) )diag[e t ] = diag[e t ]D E(αα e α h(t) ) diag[e t ]D ,
where h(t) = D (e t -1) and e t denotes an M -vector with the m-th component equal to e tm . This gives

S(t) = diag[e t ]D C α (h(t)) diag[e t ]D , (23)
which is a diagonal form similar (and equivalent for t = 0) to ( 11) since the generalized covariance matrix C α (h) of independent sources is diagonal (see ( 42) in Appendix B.1). Therefore, the generalized S-covariance matrices, estimated at different processing points t, can be used as a substitute of the T-cumulant tensors in the moment matching framework. Interestingly enough, the T-cumulant tensor (10) can be approximated by the generalized covariance matrix via its directional derivative (see Appendix B.5).

CCA generalized covariance matrices. For the CCA models ( 4)-( 6), straightforward generalizations of the ideas from Section 3.1 leads to the following definition of the generalized CCA S-covariance matrix:

S 12 (t) := E(x 1 x 2 e t x ) E(e t x ) - E(x 1 e t x ) E(e t x ) E(x 2 e t x ) E(e t x ) , (24)
where the vectors x and t are obtained by vertically stacking x 1 & x 2 and t 1 & t 2 as in (8). In the discrete CCA case, S 12 (t) is essentially the upper right block of the generalized S-covariance matrix S(t) of DICA and has the form

S 12 (t) = diag[e t1 ]D 1 C α (h(t)) diag[e t2 ]D 2 , ( 25
)
where h(t) = D (e t -1) and the matrix D is obtained by vertically stacking D 1 & D 2 by analogy with (8).

For non-Gaussian CCA, the diagonal form is

S 12 (t) = D 1 C α (h(t)) D 2 , (26)
where h(t) = D 1 t 1 + D 2 t 2 . Finally, for mixed CCA,

S 12 (t) = D 1 C α (h(t)) diag[e t2 ]D 2 , (27)
where h(t) = D 1 t 1 + D 2 (e t2 -1). Since the generalized covariance matrix of the sources C α (•) is diagonal, expressions ( 25)-( 27) have the desired diagonal form (see Appendix B.4 for detailed derivations).

Joint diagonalization algorithms

The standard algorithms such as TPM or orthogonal joint diagonalization cannot be used for the estimation of D 1 and D 2 . Indeed, even after whitening, the matrices appearing in the diagonal form ( 16)&( 18) or ( 25)-( 27) are not orthogonal. As an alternative, we use Jacobi-like non-orthogonal diagonalization algorithms [START_REF] Fu | Simultaneous diagonalization with similarity transformation for non-defective matrices[END_REF][START_REF] Iferroudjene | A new Jacobi-like method for joint diagonalization of arbitrary non-defective matrices[END_REF][START_REF] Luciani | Joint eigenvalue decomposition using polar matrix factorization[END_REF]. These algorithms are discussed in this section and in Appendix E.

The estimation of the factor loading matrices D 1 and D 2 of the CCA models ( 4)-( 6) via non-orthogonal joint diagonalization algorithms consists of the following steps: (a) construction of a set of matrices to be jointly diagonalized (using finite sample estimators), (b) a whitening step, (c) a non-orthogonal joint diagonalization step, and (d) the final estimation of the factor loading matrices (Appendix D.4).

Matrices for diagonalization. There are two ways to construct matrices for subsequent joint diagonalization: either with the CCA S-matrices (15) and T-cumulants (17) (only DCCA) or the generalized covariance matrices (24) (D/N/MCCA). Given a dataset, these matrices are estimated using natural finite sample estimators (see Appendices C.1 and C.2).

When dealing with S-and T-cumulants, the matrices are obtained via tensor projections. We define a projection

T (v) ∈ R M1×M2 of a third-order tensor T ∈ R M1×M2×M3 onto a vector v ∈ R M3 as [T (v)] m1m2 := M3 m3=1 [T ] m1m2m3 v m3 . (28)
Note that the projection T (v) is a matrix. Therefore, given 2P vectors {v 11 , v 21 , v 12 , v 22 , . . . , v 1P , v 2P }, one can construct 2P + 1 matrices

{S 12 , T 121 (v 1p ), T 122 (v 2p ), for p = 1, . . . , P }, (29)
which have the diagonal form ( 16) and ( 18). Importantly, the tensors are never constructed (see [START_REF] Anandkumar | A spectral algorithm for latent Dirichlet allocation[END_REF][START_REF] Anandkumar | Tensor decompositions for learning latent variable models[END_REF]; [START_REF] Podosinnikova | Rethinking LDA: Moment matching for discrete ICA[END_REF] and Appendix C.2).

Alternatively to (29), the set of matrices can be constructed by estimating the generalized S-covariance matrices at P + 1 processing points 0, t 1 , . . . , t P ∈ R M1+M2 :

{S 12 = S 12 (0), S 12 (t 1 ), . . . , S 12 (t P )}, (30)
which also have the diagonal form ( 25)-( 27). It is interesting to mention the connection between the Tcumulants and the generalized S-covariance matrices. The T-cumulant can be approximated via the directional derivative of the generalized covariance matrix (see Appendix B.5). However, in general, e.g., S 12 (t) with t = [t 1 ; 0] is not exactly the same as T 121 (t 1 ) and the former can be non-zero even when the latter is zero. This is important since order-4 and higher statistics are used with the method of moments when there is a risk that an order-3 statistic is zero. In general, the use of higher-order statistics increases the sample complexity and makes the resulting expressions quite complicated. Therefore, replacing the T-cumulants with the generalized S-covariance matrices is potentially beneficial.

Whitening. The matrices

W 1 ∈ R K×M1 and W 2 ∈ R K×M2 are called whitening matrices of S 12 if W 1 S 12 W 2 = I K , (31)
where I K is the K-dimensional identity matrix. Such matrices W 1 and W 2 are only defined up to multiplication by any invertible matrix Q ∈ R K×K , since any pair of matrices W 1 = QW 1 and W 2 = Q -W 2 also satisfy (31). In fact, using higher order information (i.e. the T-cumulants or the generalized covariances for t = 0) allows to solve this ambiguity.

The whitening matrices can be computed via SVD of S 12 (see Appendix D.1). When M 1 and M 2 are too large, one can use a randomized SVD algorithm (see, e.g., [START_REF] Halko | Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions[END_REF] to avoid the construction of the large matrix S 12 and to decrease the computational time.

Non-orthogonal joint diagonalization (NOJD). For simplicity, let us consider joint diagonalization of the generalized covariance matrices (30) (the same procedure holds for the S-and T-cumulants (29); see Appendix D.2). Given the whitening matrices W 1 and W 2 , the transformation of the generalized covariance matrices (30) gives P + 1 matrices

{W 1 S 12 W 2 , W 1 S 12 (t p )W 2 , p = 1, . . . , P }, ( 32
)
where each matrix is in R K×K and has reduced dimension since K < M 1 , M 2 . In practice, finite sample estimators are used to construct (30) (see Appendices C.1 and C.2).

Due to the diagonal form ( 16) and ( 25)-( 27), each matrix in (30) has the form4 (W

1 D 1 ) diag(•) (W 2 D 2 ) .
Both D 1 and D 2 are (full) K-rank matrices and W 1 and W 2 are K-rank by construction. Therefore, the square matrices 16) and (31), we get

V 1 = W 1 D 1 and V 2 = W 2 D 2 are invertible. From (
V 1 cov(α)V 2 = I and hence V 2 = diag[var(α) -1 ]V -1 1
(the covariance matrix of the sources is diagonal and we assume they are non-deterministic, i.e. var(α) = 0). Substituting this into W 1 S 12 (t)W 2 and using the diagonal form ( 25)-( 27), we obtain that the matrices in (30

) have the form V 1 diag(•)V -1
1 . Hence, we deal with the problem of the following type: Given P non-defective (a.k.a. diagonalizable) matrices B = {B 1 , . . . , B P }, where each matrix B p ∈ R K×K , find and invertible matrix Q ∈ R K×K such that

QBQ -1 = {QB 1 Q -1 , . . . , QB P Q -1 } (33)
are (jointly) as diagonal as possible. This can be seen as a joint non-symmetric eigenvalue problem. This problem should not be confused with classical joint diagonalization problem by congruence (JDC), where Q -1 is replaced by Q , except when Q is an orthogonal matrix [START_REF] Luciani | Joint eigenvalue decomposition using polar matrix factorization[END_REF]. JDC is often used for ICA algorithms or moment matching based algorithms for graphical models when a whitening step is not desirable (see, e.g., [START_REF] Kuleshov | Tensor factorization via matrix factorization[END_REF] and references therein). However, neither JDC nor the orthogonal diagonalization-type algorithms (such as, e.g., the tensor power method [START_REF] Anandkumar | Tensor decompositions for learning latent variable models[END_REF] are applicable for the problem (33).

To solve the problem (33), we use the Jacobi-like non-orthogonal joint diagonalization (NOJD) algorithms (e.g., [START_REF] Fu | Simultaneous diagonalization with similarity transformation for non-defective matrices[END_REF][START_REF] Iferroudjene | A new Jacobi-like method for joint diagonalization of arbitrary non-defective matrices[END_REF][START_REF] Luciani | Joint eigenvalue decomposition using polar matrix factorization[END_REF]. These algorithms are an extension of the orthogonal joint diagonalization algorithms based on Jacobi (=Givens) rotations [START_REF] Golub | Matrix Computations[END_REF][START_REF] Bunse-Gerstner | Numerical methods for simultaneous diagonalization[END_REF][START_REF] Cardoso | Jacobi angles for simultaneous diagonalization[END_REF]. Due to the space constraint, the description of the NOJD algorithms is moved to Appendix E. Although these algorithms are quite stable in practice, we are not aware of any theoretical guarantees about their convergence or stability to perturbation.

Spectral algorithm. By analogy with the orthogonal case [START_REF] Cardoso | Source separation using higher order moments[END_REF][START_REF] Anandkumar | A spectral algorithm for latent Dirichlet allocation[END_REF], we can easily extend the idea of the spectral algorithm to the non-orthogonal one. Indeed, it amounts to performing whitening as before and constructing only one matrix with the diagonal structure, e.g., B = W 1 S 12 (t)W 2 for some t. Then, the matrix Q is obtained as the matrix of the eigenvectors of B. The vector t can be, e.g., chosen as t = W u, where W = [W 1 ; W 2 ] and u ∈ R K is a vector sampled uniformly at random.

This spectral algorithm and the NOJD algorithms are closely connected. In particular, when B has real eigenvectors, the spectral algorithm is equivalent to NOJD of B. Indeed, in such case, NOJD boils down to an algorithm for a non-symmetric eigenproblem [START_REF] Eberlein | A Jacobi-like method for the automatic computation of eigenvalues and eigenvectors of an arbitrary matrix[END_REF][START_REF] Ruhe | On the quadratic convergene of a generalization of the Jacobi method to arbitrary matrices[END_REF]. In practice, however, due to the presence of noise and finite sample errors, B may have complex eigenvectors. In such case, the spectral algorithm is different from NOJD. Importantly, the joint diagonalization type algorithms are known to be more stable in practice (see, e.g., [START_REF] Bach | Kernel independent component analysis[END_REF][START_REF] Podosinnikova | Rethinking LDA: Moment matching for discrete ICA[END_REF].

While deriving precise theoretical guarantees is beyond the scope of this paper, the techniques outlined by [START_REF] Anandkumar | A spectral algorithm for latent Dirichlet allocation[END_REF] for the spectral algorithm for latent Dirichlet Allocation can potentially be extended. The main difference is obtaining the analogue of the SVD accuracy (Lemma C.3, [START_REF] Anandkumar | A spectral algorithm for latent Dirichlet allocation[END_REF] for the eigen decomposition. This kind of analysis can potentially be extended with the techniques outlined in (Chapter 4, [START_REF] Stewart | Matrix Perturbation Theory[END_REF]). Nevertheless, with appropriate parametric assumptions on the sources, we expect that the above described extension of the spectral algorithm should lead to similar guarantee as the spectral algorithm of [START_REF] Anandkumar | A spectral algorithm for latent Dirichlet allocation[END_REF].

Some important implementation details, including the choice of the processing points, are discussed in Appendix D.

Experiments

Synthetic data. We sample synthetic data to have ground truth information for comparison. We sample from linear DCCA which extends linear CCA (7) such that each view is x j ∼ Poisson(D j α + F j β j ). The sources α ∼ Gamma(c, b) and the noise sources β j ∼ Gamma(c j , b j ), for j = 1, 2, are sampled from the gamma distribution (where b is the rate parameter). Let s j ∼ Poisson(D j α) be the part of the sample due to the sources and n j ∼ Poisson(F j β j ) be the part of the sample due to the noise (i.e., x j = s j + n j ). Then we define the expected sample length due to the sources and noise, respectively, as L js := E[ m s jm ] and 2015)). For the larger dimensional example (Fig. 2, right), each column of the matrices D j and F j , for j = 1, 2, is sampled from the symmetric Dirichlet distribution with the concentration parameter equal to 0.5. For the smaller 2D example (Fig. 2, left), they are fixed:

L jn := E[ m n jm ].
D 1 = D 2 with [D 1 ] 1 = [D 1 ] 2 = 0.5 and F 1 = F 2 with [F 1 ] 11 = [F 1 ] 22 = 0.9 and [F 1 ] 12 = [F 1 ] 21 = 0.1.
For each experiment, D j and F j , for j = 1, 2, are sampled once and, then, the x-observations are sampled for different sample sizes N = {500, 1, 000, 2, 000, 5, 000, 10, 000}, 5 times for each N .

Metric. The evaluation is performed on a matrix D obtained by stacking D 1 and D 2 vertically (see also the comment after Thm. 1). As in [START_REF] Podosinnikova | Rethinking LDA: Moment matching for discrete ICA[END_REF], we use as evaluation metric the normalized 1 -error between a recovered matrix D and the true matrix D with the best permutation of columns err 1 ( D, D) := min π∈PERM

1 2K k d π k -d k 1 ∈ [0, 1].
The minimization is over the possible permutations π ∈ PERM of the columns of D and can be efficiently obtained with the Hungarian algorithm for bipartite matching. The (normalized) 1 -error takes the values in [0, 1] and smaller values of this error indicate better performance of an algorithm.

M 1 = M 2 = K 1 = K 2 = 2, K = 1, c = c 1 = c 2 = 0.
1, and L s = L n = 100; middle (2D data): the x 1 -observations and factor loading matrices for the 2D example (F 1j denotes the j-th column of the noise factor matrix F 1 ); right (20D example): andc 1 = c 2 = 0 Algorithms. We compare DCCA (implementation with the S-and T-cumulants) and DCCAg (implementation with the generalized S-covariance matrices and the processing points initialized as described in Appendix D.3) to DICA and the non-negative matrix factorization (NMF) algorithm with multiplicative updates for divergence [START_REF] Lee | Algorithms for non-negative matrix factorization[END_REF]. To run DICA or NMF, we use the stacking trick (8). DCCA is set to estimate K components. DICA is set to estimate either

M 1 = M 2 = K 1 = K 2 = 20, K = 10, L s = L n = 1, 000, c = 0.3,
K 0 = K + K 1 + K 2 or M = M 1 + M 2 compo-
nents (whichever is the smallest, since DICA cannot work in the over-complete case). NMF is always set to estimate K 0 components. For the evaluation of DICA/NMF, the K columns with the smallest 1 -error are chosen. NMF • stands for NMF initialized with a matrix D of the form (8) with induced zeros; otherwise NMF is initialized with (uniformly) random non-negative matrices.

Synthetic experiment. We first perform an experiment with discrete synthetic data in 2D (Fig. 2) and then repeat the same experiment when the size of the problem is 10 times larger. In practice, we observed that for K 0 < M all models work approximately equally well, except for NMF which breaks down in high dimensions. In the over-complete case as in Fig. 2, DCCA works better. A continuous analogue of this experiment is presented in Appendix F.1.

Real data (translation). Following Vinokourov et al. (2002), we illustrate the performance of DCCA by extracting bilingual topics from the Hansard collection [START_REF] Vinokourov | A probabilistic framework for the hierarchic organisation and classification of document collections[END_REF] with aligned English and French proceedings of the 36-th Canadian Parliament. We first pre-process the dataset with the NLTK toolbox by [START_REF] Bird | Natural Language Processing with Python[END_REF] (see details in Appendix F.3) to obtain N = 12, 932 aligned documents with M 1 = M 2 = 5, 000 English and French words. We then run DCCA with K = 20 to extract the aligned topics. Some of the extracted topics are presented in Table 1 and all the topics as well as the detailed description of this experiment are presented in Appendix F.3.

Running time. For the real experiment above, the runtime of DCCA algorithm is 24 seconds including 22 seconds for SVD at the whitening step. In general, the computational complexity of the D/N/MCCA algorithms is bounded by the time of SVD plus O(RN K) + O(N K 2 ), where R is the largest number of non-zero components in the stacked vector x = [x 1 ; x 2 ], plus the time of NOJD for the matrices of size K.

In practice, DCCAg is faster than DCCA.

Conclusion

We have proposed the first identifiable versions of CCA, together with moment matching algorithms which allow the identification of the loading matrices in a semi-parametric framework, where no assumptions are made regarding the distribution of the source or the noise. We also introduce a new sets of moments (our generalized covariance matrices), which could prove useful in other settings.

Acknowledgments. This work was partially supported by the MSR-Inria Joint Center.

Appendix

The appendix is organized as follows.

-In Appendix A, we present the proof of Theorem 1 stating the identifiability of the CCA models (4)-( 6).

-In Appendix B, we provide some details for the generalized covariance matrices: the form of the generalized covariance matrices of independent variables (Appendix B.1), the derivations of the diagonal form of the generalized covariance matrices of discrete ICA (Appendix B.3), the derivations of the diagonal form of the generalized covariance matrices of the CCA models ( 4)-( 6) (Appendix B.4), and approximation of the T-cumulants with the generalized covariance matrix (Appendix B.5).

-In Appendix C, we provide expressions for natural finite sample estimators of the generalized covariance matrices and the T-cumulant tensors for the considered CCA models.

-In Appendix D, we discuss some rather technical implementation details: computation of whitening matrices (Appendix D.1), selection of the projection vectors for the T-cumulants and the processing points for the generalized covariance matrices (Appendix D.3), and the final estimation of the factor loading matrices (Appendix D.4).

-In Appendix E, we describe the non-orthogonal joint diagonalization algorithms used in this paper.

-In Appendix F, we present some supplementary experiments: a continuous analog of the synthetic experiment from Section 5 (Appendix F.1), an experiment to analyze the sensitivity of the DCCA algorithm with the generalized S-covariance matrices to the choice of the processing points (Appendix F.2), and a detailed description of the experiment with the real data from Section 5 (Appendices F.3 and F.4).

A Identifiability

In this section, we prove that the factor loading matrices D 1 and D 2 of the non-Gaussian CCA (4), discrete CCA ( 5), and mixed CCA ( 6) models are identifiable up to permutation and scaling if at most one source α k is Gaussian. We provide a complete proof for the non-Gaussian CCA case and show that the other two cases can be proved by analogy.

A.1 Identifiability of non-Gaussian CCA (4)

The proof uses the notion of the second characteristic function (SCF) of a random variable x ∈ R M :

φ x (t) = log E(e it x ),
for all t ∈ R M . The SCF completely defines the probability distribution of x (see, e.g., Jacod & Protter, 2004). Important difference between the SCF and the cumulant generating function ( 19) is that the former always exists.

The following property of the SCF is of central importance for the proof: if two random variables, z 1 and z 2 , are independent, then φ A1z1+A2z2 (t) = φ z1 (A 1 t) + φ z2 (A 2 t), where A 1 and A 2 are any matrices of compatible sizes.

We can now use our CCA model to derive an expression of φ x (t). Indeed, defining a vector x by stacking the vectors x 1 and x 2 , the SCF of x for any t = [t 1 ; t 2 ], takes the form

φ x (t) = log E(e it 1 x1+it 2 x2 ) (a) = log E(e iα (D 1 t1+D 2 t2)+iε 1 t1+iε 2 t2 ) (b) = log E(e iα (D 1 t1+D 2 t2) ) + log E(e iε 1 t1 ) + log E(e iε 2 t2 ) = φ α (D 1 t 1 + D 2 t 2 ) + φ ε1 (t 1 ) + φ ε2 (t 2 ),
where in (a) we substituted the definition (4) of x 1 and x 2 and in (b) we used the independence α ⊥ ⊥ ε 1 ⊥ ⊥ ε 2 . Therefore, the blockwise mixed derivatives of φ x are equal to

∂ 1 ∂ 2 φ x (t) = D 1 φ α (D 1 t 1 + D 2 t 2 )D 2 , ( 34
)
where

∂ 1 ∂ 2 φ x (t) := ∇ t1 ∇ t2 φ x (h(t 1 , t 2 )) ∈ R M1×M2 and φ α (u) := ∇ 2 u φ α (u)
, does not depend on the noise vectors ε 1 and ε 2 .

For simplicity, we first prove the identifiability result when all components of the common sources are non-Gaussian. The high level idea of the proof is as follows. We assume two different representations of x 1 and x 2 and using (34) and the independence of the components of α and the noises, we first show that the two potential dictionaries are related by an orthogonal matrix (and not any invertible matrix), and then show that this implies that the two potential sets of independent components are (orthogonal) linear combinations of each other, which, for non-Gaussian components which are not reduced to point masses, imposes that this orthogonal transformation is the combination of a permutation matrix and marginal scaling-a standard result from the ICA literature (Comon, 1994, Theorem 11).

Let us then assume that two equivalent representations of non-Gaussian CCA exist:

x 1 = D 1 α + ε 1 = E 1 β + η 1 , x 2 = D 2 α + ε 2 = E 2 β + η 2 , ( 35
)
where the other sources β = (β 1 , . . . , β K ) are also assumed mutually independent and non-degenerate. As a standard practice in the ICA literature and without loss of generality as the sources have non-degenerate components, one can assume that the sources have unit variances, i.e. cov(α, α) = I and cov(β, β) = I, by respectively rescaling the columns of the factor loading matrices. Under this assumption, the two expressions of the cross-covariance matrix are

cov(x 1 , x 2 ) = D 1 D 2 = E 1 E 2 , (36)
which, given that D 1 , D 2 have full rank, implies that 5

E 1 = D 1 Q, E 2 = D 2 Q -, (37)
where Q ∈ R K×K is some invertible matrix. Substituting the representations (35) into the blockwise mixed derivatives of the SCF (34) and using the expressions (37) give

D 1 φ α (D 1 t 1 + D 2 t 2 )D 2 = D 1 Qφ β (Q D 1 t 1 + Q -1 D 2 t 2 )Q -1 D 2 ,
for all t 1 ∈ R M1 and t 2 ∈ R M2 . Since the matrices D 1 and D 2 have full rank, this can be rewritten as

φ α (D 1 t 1 + D 2 t 2 ) = Qφ β (Q D 1 t 1 + Q -1 D 2 t 2 )Q -1 ,
which holds for all t 1 ∈ R M1 and t 2 ∈ R M2 . Moreover, still since D 1 and D 2 have full rank, we have, for any u 1 , u 2 ∈ R K the existence of t 1 ∈ R M1 and t 2 ∈ R M2 , such that u 1 = D 1 t 1 and u 2 = D 2 t 2 , that is,

φ α (u 1 + u 2 ) = Qφ β (Q u 1 + Q -1 u 2 )Q -1 , (38)
for all u 1 , u 2 ∈ R K .

We will now prove two facts:

(F1) For any vector v ∈ R K , then φ β ((Q Q -I)v) = -I, which will imply that QQ = I because of the non-Gaussian assumptions.

(F2) If QQ = I, then φ α (u) = φ Qβ (u) for any u ∈ R K , which will imply that Q is the composition of a permutation and a scaling. This will end the proof.

Proof of fact (F1). By letting u 1 = Qv and u 2 = -Qv, we get:

φ α (0) = Qφ β ((Q Q -I)v)Q -1 , (39)
Since 6 φ α (0) = -cov(α) = -I, one gets

φ β ((Q Q -I)v) = -I,
for any v ∈ R K .

5 The fact that D 1 , D 2 have full rank and that E 1 , E 2 have K columns, combined with (36), implies that E 1 , E 2 have also full rank.

6 Note that ∇ 2 u φα(u) = -

E(αα e iu α )
E(e iu α )

+ Eα(u)Eα(u) , where Eα(u) =

E(αe iu α ) E(e iu α )
.

Using the property that φ A β (v) = A φ β (Av)A for any matrix A, and in particular with

A = Q Q -I, we have that φ A β (v) = -A A, i.e. is constant.
If the second derivative of a function is constant, the function is quadratic. Therefore, φ A β (•) is a quadratic function. Since the SCF completely defines the distribution of its variable (see,e.g., Jacod & Protter ( 2004)), A β must be Gaussian (the SCF of a Gaussian random variable is a quadratic function). Given Lemma 9 from Comon (1994) (i.e., Cramer's lemma: a linear combination of non-Gaussian random variables cannot be Gaussian unless the coefficients are all zero), this implies that A = 0, and hence Q Q = I, i.e., Q is an orthogonal matrix.

Proof of fact (F2). Plugging Q = Q -1 into (38), with u 1 = 0 and u 2 = u, gives

φ α (u) = Qφ β (Q u)Q = φ Qβ (u), (40)
for any u ∈ R K . By integrating both sides of (40) and using φ α (0) = φ Qβ (0) = 0, we get that φ α (u) = φ Qβ (u) + iγ u for all u ∈ R K for some constant vector γ. Using again that the SCF completely defines the distribution, it follows that α -γ and Qβ have the same distribution. Since both α and β have independent components, this is only possible when Q = ΛP , where P is a permutation matrix and Λ is some diagonal matrix (Comon, 1994, Theorem 11).

A.2 Case of a single Gaussian source

Without loss of generality, we assume that the potential Gaussian source is the first one for α and β. The first change is in the proof of fact (F1). We use the same argument up to the point where we conclude that A β is a Gaussian vector. As only β 1 can be Gaussian, Cramer's lemma implies that only the first row of A can have non-zero components, that is

A = Q Q -I = e 1 f
, where e 1 is the first basis vector and f any vector. Since Q Q is symmetric, we must have

Q Q = I + ae 1 e 1 ,
where a is a constant scalar different than -1 as Q Q is invertible. This implies that Q Q is an invertible diagonal matrix Λ, and hence QΛ -1/2 is an orthogonal matrix, which in turn implies that

Q -1 = Λ -1 Q .
Plugging this into (38) gives, for any u 1 and u 2 :

φ α (u 1 + u 2 ) = Qφ β (Q u 1 + Λ -1 Q u 2 )Λ -1 Q .
Given that diagonal matrices commute and that φ β is diagonal for independent sources (see Appendix B.1), this leads to

φ α (u 1 + u 2 ) = QΛ -1/2 φ β (Q u 1 + Λ -1 Q u 2 )Λ -1/2 Q .
For any given v ∈ R K , we are looking for u 1 and u 2 such that Q u

1 + Λ -1 Q u 2 = Λ -1/2 Q v and u 1 + u 2 = v, which is always possible by setting Q u 2 = (Λ -1/2 + I) -1 Q v and Q u 1 = Q v -Q u 2
by using the special structure of Λ. Thus, for any v,

φ α (v) = QΛ -1/2 φ β (Λ -1/2 Q v)Λ -1/2 Q = φ QΛ -1/2 β (v).
Integrating as previously, this implies that the characteristic function of α and QΛ -1/2 β differ only by a linear function iγ v, and thus, that α -γ and QΛ -1/2 β have the same distribution. This in turn, from Comon (1994, Theorem 11), implies that QΛ -1/2 is a product of a scaling and a permutation, which ends the proof.

A.3 Identifiability of discrete CCA (5) and mixed CCA (6)

Given the discrete CCA model, the SCF φ x (t) takes the form

φ x (t) = φ α (D 1 (e it1 -1) + D 2 (e it2 -1)) + φ ε1 (e it1 -1) + φ ε2 (e it2 -1),
where e itj , for j = 1, 2, denotes a vector with the m-th element equal to e i[tj ]m , and we used the arguments analogous with the non-Gaussian case. The rest of the proof extends with a correction that sometimes one has to replace D j with diag[e itj ]D j and that u j = D j (e itj -1) for j = 1, 2. For the mixed CCA case, only the part related to x 2 and D 2 changes in the same way as for the discrete CCA case.

B The generalized expectation and covariance matrix B.1 The generalized expectation and covariance matrix of the sources

The sources α = (α 1 , . . . , α K ) are mutually independent. Therefore, for some h ∈ R K , their CGF (19) K α (h) = log E(e α h ) takes the form

K α (h) = k log E(e α k h k ) .
Therefore, the k-th element of the generalized expectation (20

) of α is (separable in α k ) [E α (h)] k = E(α k e α k h k ) E(e α k h k ) (41)
and the generalized covariance (21) of α is diagonal due to the separability and its k-th diagonal element is

[C α (h)] kk = E(α 2 k e α k h k ) E(e α k h k ) -[E α (h)] 2 k . (42)
B.2 Some expectations of a Poisson random variable Let x ∈ R M be a multivariate Poisson random variable with mean y ∈ R M + . Then, for some t ∈ R M , E(e t x ) = e y (e t -1) , E(x m e t x ) = y m e tm e y (e t -1) , E(x 2 m e t x ) = y m e tm + 1 y m e tm e y (e t -1) , E(x m x m e t x ) = y m e tm y m e t m e y (e t -1) , m = m , where e t denotes an M -vector with the m-th element equal to e tm .

B.3 The generalized expectation and covariance matrix of discrete ICA

In this section, we use the expectations of a Poisson random variable presented in Appendix B.2.

Given the discrete ICA model ( 9), the generalized expectation (20) of x ∈ R M takes the form

E x (t) = E(xe t x ) E(e t x ) = E E(xe t x |α) E E(e t x |α) = diag[e t ]D E(αe α h(t) ) E(e α h(t) ) = diag[e t ]DE α (h(t)),
where t ∈ R M is a parameter, h(t) = D (e t -1), and e t denotes an M -vector with the m-th element equal to e tm . Note that in the last equation we used the definition (20) of the generalized expectation E α (•).

Further, the generalized covariance (21) of x takes the form

C x (t) = E(xx e t x ) E(e t x ) -E x (t)E x (t) = E E(xx e t x |α) E E(e t x |α) -E x (t)E x (t) .
Plugging into this expression the expression for E x (t) and

E(xx e t x |α) = diag[e t ]DE(αα e α h(t) )D diag[e t ] + diag[e t ]diag DE(αe α h(t) ) we get C x (t) = diag[E x (t)] + diag[e t ]DC α (h(t))D diag[e t ],
where we used the definition (21) of the generalized covariance of α.

B.4 The generalized CCA S-covariance matrix

In this section we sketch the derivation of the diagonal form ( 27) of the generalized S-covariance matrix of mixed CCA (6). Expressions ( 25) and ( 26) can be obtained in a similar way.

Denoting x = [x 1 ; x 2 ] and t = [t 1 ; t 2 ] (i.e. stacking the vectors as in ( 8)), the CGF (19) of mixed CCA (6) can be written as

K x (t) = log E(e t 1 x1+t 2 x2 ) = log E E(e t 1 x1+t 2 x2 | α, ε 1 , ε 2 ) (a) = log E E(e t 1 x1 | α, ε 1 )E(e t 2 x2 | α, ε 2 ) (b)
= log E e t 1 (D1α+ε1) e (D2α+ε2) (e t 2 -1)

(c)
= log E e α h(t) + log E e ε 2 (e t 2 -1) + log E(e t 1 ε1 ),

where h(t) = (D 1 t 1 + D 2 (e t2 -1), in (a) we used the conditional independence of x 1 and x 2 , in (b) we used the first expression from Appendix B.2, and in (c) we used the independence assumption (3).

The generalized CCA S-covariance matrix is defined as

S 12 (t) := ∇ t2 ∇ t1 K x (t).
Its gradient with respect to t 1 is

∇ t1 K x (t) = D 1 E(αe α h(t) ) E(e α h(t) ) + E(ε 1 e t 1 ε1 ) E(e t 1 ε1 )
,

where the last term does not depend on t 2 . Computing the gradient of this expression with respect to

t 2 gives S 12 (t) = D 1 C α (h(t)) diag[e t2 ]D 2 ,
where we substituted expression (42) for the generalized covariance of the independent sources.

B.5 Approximation of the T-cumulants with the generalized covariance matrix

Let f mm (t) = [C x (t)] mm be a function R → R M corresponding to the (m, m )-th element of the generalized covariance matrix. Then the following holds for its directional derivative at t 0 along the direction t:

∇f mm (t 0 ), t = lim δ→0 f mm (t 0 + δt) -f mm (t 0 ) δ ,
where •, • stands for the inner product. Therefore, when using the fact that ∇f (t 0 ) = ∇C x (t) is the generalized cumulant of x at t 0 and the definition of a projection of a tensor onto a vector (28), one obtains for t 0 = 0 the approximation of the cumulant cum(x) with the generalized covariance matrix C x (t).

Let us define v 1 = W 1 u 1 and v 1 = W 2 u 2 for some u 1 , u 2 ∈ R K . Then, approximations for the Tcumulants (17) of discrete CCA take the following form: W 1 T 121 (v 1 )W 2 is approximated by the generalized S-covariances (24) S 12 (t) via the following expression

W 1 T 121 (v 1 )W 2 ≈ W 1 S 12 (δt 1 )W 2 -W 1 S 12 (0)W 2 δ -W 1 diag(v 1 )S 12 W 2 ,
where t 1 = v 1 0 and W 1 T 122 (v 2 )W 2 is approximated by the generalized S-covariances S 12 (t) via

W 1 T 122 (v 2 )W 2 ≈ W 1 S 12 (δt 2 )W 2 -W 1 S 12 (0)W 2 δ -W 1 S 12 diag(v 2 )W 2 ,
where t 2 = 0 v 2 and δ are chosen to be small.

C Finite sample estimators

C.1 Finite sample estimators of the generalized expectation and covariance matrix Following [START_REF] Yeredor | Blind source separation via the second characteristic function[END_REF]; Slapak & Yeredor (2012b), we use the most direct way of defining the finite sample estimators of the generalized expectation (20) and covariance matrix (21).

Given a finite sample X = {x 1 , x 2 , . . . , x N }, an estimator of the generalized expectation is

E x (t) = n x n w n n w n
where weights w n = e t xn and an estimator of the generalized covariance is

C x (t) = n x n x n w n n w n -E x (t) E x (t) .
Similarly, an estimator of the generalized S-covariance matrix is then

C x1,x2 (t) = n x 1n x 2n w n n w n - n x 1n w n n w n n x 2n w n n w n , where x = [x 1 ; x 2 ] and t = [t 1 ; t 2 ] for some t 1 ∈ R M1 and t 2 ∈ R M2 .
Some properties of these estimators are analyzed by Slapak & Yeredor (2012b).

C.2 Finite sample estimators of the DCCA cumulants

In this section, we sketch the derivation of unbiased finite sample estimators for the CCA cumulants S 12 , T 121 , and T 122 . Since the derivation is nearly identical to the derivation of the estimators for the DICA cumulants (see Appendix F.2 of [START_REF] Podosinnikova | Rethinking LDA: Moment matching for discrete ICA[END_REF]), all details are omitted.

Given a finite sample X 1 = {x 11 , x 12 , . . . , x 1N } and X 2 = {x 21 , x 22 , . . . , x 2N }, the finite sample estimator of the discrete CCA S-covariance (15), i.e., S 12 := cum(x 1 , x 2 ), takes the form

S 12 = η 1 X 1 X 2 -N E(x 1 ) E(x 2 ) , (43)
where E(x 1 ) = N -1 n x 1n , E(x 2 ) = N -1 n x 2n , and η 1 = 1/(N -1). Substitution of the finite sample estimators of the 2nd and 3rd cumulants (see, e.g., Appendix C.4 of [START_REF] Podosinnikova | Rethinking LDA: Moment matching for discrete ICA[END_REF]) into the definition of the DCCA T-cumulants (17) leads to the following expressions

W 1 T 12j (v j ) W 2 = η 2 [( W 1 X 1 )diag(X j v j )] ⊗ ( W 2 X 2 ) + η 2 v j , E(x j ) 2N [ W 1 E(x 1 )] ⊗ [ W 2 E(x 2 )] -η 2 v j , E(x j ) ( W 1 X 1 ) ⊗ ( W 2 X 2 ) -η 2 [( W 1 X 1 )(X j v j )] ⊗ [ W 2 E(x 2 )] -η 2 [ W 1 E(x 1 )] ⊗ [( W 2 X 2 )(X j v j )] -η 1 ( W (j) 1 X 1 ) ⊗ ( W (j) 2 X 2 ) + η 1 N [ W (j) 1 E(x 1 )] ⊗ [ W (j) 2 E(x 2 )],
where η 2 = N/((N -1)(N -2)) and W

(1) 1

= W 1 diag(v 1 ), W (1) 2 = W 2 , W (2) 1 = W 1 , and W (2) 2 = W 2 diag(v 2 ).
In the expressions above, W 1 and W 2 denote whitening matrices of S 12 , i.e. such that W 1 S 12 W 2 = I.

D Implementation details D.1 Computation of whitening matrices

One can compute such whitening matrices (31) via the singular value decomposition (SVD) of S 12 . Let S 12 = U ΣV be the SVD of S 12 , then one can define W 1 = U 1:K Λ and W 2 = V 1:K Λ, where U 1:K and V 1:K are the first K left-and right-singular vectors and Λ = diag(σ

-1/2 1 , . . . , σ -1/2 K
) and σ 1 , . . . , σ K are the K largest singular values.

Although SVD is computed only once, the size of the matrix S 12 can be significant even for storage. To avoid construction of this large matrix and speed up SVD, one can use randomized SVD techniques [START_REF] Halko | Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions[END_REF]. Indeed, since the sample estimator S 12 has the form (43), one can reduce this matrix by sampling two Gaussian random matrices Ω 1 ∈ R K×M1 and Ω 2 ∈ R K×M2 , where K is slightly larger than K. Now, if U and V are the K largest singular vectors of the reduced matrix Ω 1 S 12 Ω 2 , then Ω † 1 U and Ω † 2 V are approximately (and up to permutation and scaling of the columns) the K largest singular vectors of S 12 .

D.2 Applying whitening transform to DCCA T-cumulants

Transformation of the T-cumulants (29) with whitening matrices W 1 and W 2 gives new tensors T 12j ∈ R K×K×K :

T 12j := T 12j × 1 W 1 × 2 W 2 × 3 W j , (44)
where j = 1, 2. Combining this transformation with the projection (28), one obtains 2P + 1 matrices

W 1 S 12 W 2 , W 1 T 12j (W j u jp )W 2 , (45)
where p = 1, . . . , P and j = 1, 2 and we used v jp = W j u jp to take into account whitening along the third direction. By choosing u jp ∈ R K to be the canonical vectors of the R K , the number of tensor projections is reduced from M = M 1 + M 2 to 2K.

D.3 Choice of projection vectors or processing points

For the T-cumulants (29), we choose the K projection vectors as v 1p = W 1 e p and v 2p = W 2 e p , where e p is one of the columns of the K-identity matrix (i.e., a canonical vector). For the generalized S-covariances (30), we choose the processing points as t 1p = δ 1 v 1p and t 2p = δ 2 v 2p , where δ j , for j = 1, 2 are set to a small value such as 0.1 divided by m E(|x jm |)/M j , for j = 1, 2.

When projecting a tensor T 12j onto a vector, part of the information contained in this tensor gets lost. To preserve all information, one could project a tensor T 12j onto the canonical basis of R Mj to obtain M j matrices. However, this would be an expensive operation in terms of both memory and computational time.

In practice, we use the fact, that the tensor T 12j , for J = 1, 2, is transformed with whitening matrices (44). Hence, the projection vector has to include multiplication by the whitening matrices. Since they reduce the dimension to K, choosing the canonical basis in R K becomes sufficient. Hence, the choice v 1p = W 1 e p and v 2p = W 2 e p , where e p is one of the columns of the K-identity matrix.

Importantly, in practice, the tensors are never constructed (see Appendix C.2).

The choice of the processing points of the generalized covariance matrices has to be done carefully. Indeed, if the values of t 1 or t 2 are too large, the exponents blow up. Hence, it is reasonable to maintain the values of the processing points very small. Therefore, for j = 1, 2, we set t jp = δ j v jp where δ j is proportional to a parameter δ which is set to a small value (δ = 0.1 by default), and the scale is determined by the inverse of the empirical average of the component of x j , i.e.: The non-orthogonal joint diagonalization algorithm outputs an invertible matrix Q. If the estimated factor loading matrices are not supposed to be non-negative (continuous case of NCCA (4)), then

δ j := δ N M j N n=1 Mj m=1 [|X j |] mn , (46)
D 1 = W † 1 Q, D 2 = W † 2 Q -1 , (47)
where † stands for the pseudo-inverse. For the spectral algorithm, where Q are eigenvectors of a nonsymmetric matrix and are not guaranteed to be real, only real parts are kept after evaluating matrices D 1 and D 2 in accordance with (47).

If the matrices D 1 and/or D 2 have to be non-negative (the discrete case of DCCA ( 5) and MCCA ( 6)), they have to be further mapped. For that, we select the sign of each column such that the vector (column) has less negative than positive components, which is measured by the sum of squares of the components of each sign, (this is necessary since the scaling unidentifiability includes the scaling by -1) and then truncate all negative values at 0.

In practice, due to the scaling unidentifiability, each column of the obtained matrices D 1 and D 2 can be further normalized to have the unit 1 -norm. This is applicable in all cases (D/M/NCCA).

E Jacobi-like joint diagonalization of non-symmetric matrices

Given N non-defective (a.k.a. diagonalizable) not necessary normal7 matrices

A = {A 1 , A 2 , . . . , A N } , where each matrix A n ∈ R M ×M , find such matrix Q ∈ R M ×M that matrices Q -1 AQ = Q -1 A 1 Q, Q -1 A 2 Q, . . . , Q -1 A N Q
are (jointly) as diagonal as possible. We refer to this problem as a non-orthogonal joint diagonalization (NOJD) problem.8

Algorithm 1 Non-orthogonal joint diagonalization (NOJD)

1: Initialize: A (0) ← A and Q (0) ← I M and iterations = 0 2: for sweeps k = 1, 2, . . . do 3:

for p = 1, . . . , M -1 do 4:

for q = p + 1, . . . , M do 5:

Increase = + 1 6:

Find the (approx.) shear parameter y * defined in (54) 7:

Find the Jacobi angle θ * defined in (53) The data is the discrete synthetic data as described in Section 5 with the parameters set as in Fig. 2 (right).

8: Update Q ( ) ← Q ( -1) S ( ) * U ( ) * 9: Update A ( ) ← U
The three NOJD algorithms can have slightly different convergence properties, however, for the purposes of this paper their performance can hardly be distinguished. That is, the difference in the performance of the algorithms in terms of the 1 -error of the factor loading matrices is hardly noticeable. For the experiments, we use the JDTM algorithm, the other two algorithms could be equally used. To the best of our knowledge, no theoretical analysis of the NOJD algorithms is available, except for the single matrix case when they boil down to the (non-symmetric) eigenproblem [START_REF] Eberlein | A Jacobi-like method for the automatic computation of eigenvalues and eigenvectors of an arbitrary matrix[END_REF][START_REF] Ruhe | On the quadratic convergene of a generalization of the Jacobi method to arbitrary matrices[END_REF].

The following intuitively explains why the normality measure, i.e. the sum of the squared Frobenius norms, has to be minimized at the shear transform. As [START_REF] Ruhe | On the quadratic convergene of a generalization of the Jacobi method to arbitrary matrices[END_REF]) mention, for every matrix A and non-singular

Q: inf Q Q -1 AQ 2 F = Λ 2 F ,
where Λ is the diagonal matrix containing the eigenvalues of A. Therefore, a diagonalized version of the matrix A must have the smallest Frobenius norm. Since the unitary transform does not change the Frobenius norm, it can only be minimized with the shear transform. Further, if a matrix is normal, i.e. A A = AA with a symmetric matrix as a particular case, the upper triangular matrix in its Schur decomposition is zero (Golub & Van Loan, 1996, Chapter 7) and then the Schur vectors correspond to the (orthogonal in this case) eigenvectors of this matrix. Therefore, a normal non-defective matrix can be diagonalized by an orthogonal matrix, which preserves the Frobenius norm. Hence, the shear transform by minimizing the normality measure decreases the deviation from normality and then the unitary transform by minimizing the diagonality measure decreases the deviation from diagonality.

F Supplementary experiments F.1 Continuous synthetic data

This experiment is essentially a continuous analogue to the synthetic experiment with the discrete data from Section 5.

Synthetic data. We sample synthetic data from the linear non-Gaussian CCA (NCCA) model ( 7) with each view x j = D j α + F j β j . The (non-Gaussian) sources are α ∼ z α Gamma(c, b), where z α is a Rademacher random variable (i.e., takes the values -1 or 1 with the equal probabilities). The noise sources are β j ∼ z βj Gamma(c j , b j ), for j = 1, 2, where again z βj is a Rademacher random variable. Parameters of the gamma distribution are initialized by analogy with the discrete case (see Section 5). The elements of the matrices D j and F j , for j = 1, 2, are sampled i.i.d. for the uniform distribution in [-1, 1]. Each column of D j and F j , for j = 1, 2, is normalized to have the unit 1 -norm.

Algorithms. We compare gNCCA (the implementation of NCCA with the generalized S-covariance matrices with the default values of the parameters δ 1 and δ 2 as described in Appendix D.3) the spectral algorithm for NCCA (also with the generalized S-covariance matrices) to the JADE algorithm11 [START_REF] Cardoso | Blind beamforming for non Gaussian signals[END_REF] for independent component analysis (ICA) and to classical CCA.

Synthetic experiment. In Fig. 3 (left and middle), the results of the experiment for the different number of topics are presented. The error of the classical CCA is high due to the mentioned unidentifiability issues.

F.2 Sensitivity of the generalized covariance matrices to the choice of the processing points

In this section, we experimentally analyze the performance of the DCCAg algorithm based on the generalized S-covariance matrices vs. the parameters δ 1 and δ 2 . We use the experimental setup of the synthetic discrete data from Section 5 with K 1 = K 2 = K = 10. The results are presented in Fig. 3 (right).

F.3 Real data experiment -translation topics

For the real data experiment, we estimate the factor loading matrices (topics, in the following) D 1 and D 2 of aligned proceedings of the 36-th Canadian Parliament in English and French languages. This Hansard collection can be found at http://www.isi.edu/natural-language/download/hansard/.

Although going into details of natural language processing (NLP) related problems is not the goal of this paper, we do minor pre-processing (see Appendix F.4) of this text data to improve the presentation of the estimated bilingual topics D 1 and D 2 .

The 20 topics obtained with DCCA are presented in Tables 23456. For each topic, we display the 20 most frequent words (ordered from top to bottom in the decreasing order). Most of the topic have quite clear interpretation. Moreover, we can often observe the pairs of words which are each others translations in the topics. Take, e.g.,

-the topic 10: the phrase "pension plan" can be translated as "régime de retraite", the word "benefits" as "prestations", and abbreviations "CPP" and "RPC" stand for "Canada Pension Plan" and "Régime de pensions du Canada", respectively;

-the topic 3: "OTAN" is the French abbreviation for "NATO", the word "war" is translated as "guerre", and the word "peace" as "paix"; -the topic 9: "Nisga" is the name of an Indigenous (or "aboriginal") people in British Columbia, the word "aboriginal" translates to French as "autochtontes", and, e.g., the word "right" can be translated as "droit".

Note also that, e.g., in topics 10, although the French words "ans" and "années" are present in the French topic, their English translation "year" is not, since it was removed as one of the 15 most frequent words in English (see Appendix F.4).

F.4 Data preprocessing

For the experiment, we use House Debate Training Set of the Hansard collection. To process this text data, we perform case conversion, stemming, and removal of some stop words. For stemming, the SnowballStemmer of the NLTK toolbox by [START_REF] Bird | Natural Language Processing with Python[END_REF] was used for both English and French languages. Although this stemmer has particular problems (such as mapping several different forms of a word to a single stem in one language but not in the other), they are left beyond our consideration. Moreover, in addition to the standard stop words of the NLTK toolbox, we also removed the following words that we consider to be stop words for our task12 (and their possible forms):

-from English: ask, become, believe, can, could, come, cost, cut, do, done, follow, get, give, go, know, let, like, listen, live, look, lost, make, may, met, move, must, need, put, say, see, show, take, think, talk, use, want, will, also, another, back, day, certain, certainly, even, final, finally, first, future, general, good, high, just, last, long, major, many, new, next, now, one, point, since, thing, time, today, way, well, without; -from French (translations in brackets): demander (ask), doit (must), devenir (become), dit (speak, talk), devoir (have to), donner (give), ila (he has), met (put), parler (speak, talk), penser (think), pourrait (could), pouvoir (can), prendre (take), savoir (know), aller (go), voir (see), vouloir (want), actuellement, après (after), aujourd'hui (today), autres (other), bien (good), beaucoup (a lot), besoin (need), cas (case), cause, cela (it), certain, chose (thing), déjà (already), dernier (last), égal (equal), entre (between), fac ¸on (way), grand (big), jour (day), lorsque (when), neuf (new), passé (past), plus, point, présent, prêts (ready), prochain (next), quelque (some), suivant (next), unique.

  For sampling, the target values L s = L 1s = L 2s and L n = L 1n = L 2n are fixed and the parameters b and b j are accordingly set to ensure these values: b = Kc/L s and b j = K j c j /L n (see Appendix B.2 of Podosinnikova et al. (

Figure 2 :

 2 Figure 2: Synthetic experiment with discrete data. Left (2D example):M 1 = M 2 = K 1 = K 2 = 2, K = 1, c = c 1 = c 2 = 0.1, and L s = L n = 100; middle (2D data): the x 1 -observations and factor loading matrices for the 2D example (F 1j denotes the j-th column of the noise factor matrix F 1 ); right (20D example): M 1 = M 2 = K 1 = K 2 = 20, K = 10, L s = L n = 1, 000, c = 0.3, and c 1 = c 2 = 0.1.

  for j = 1, 2. See Appendix F.2 for an experimental comparison of different values of δ (the default value used in other experiments is δ = 0.1). D.4 Finalizing estimation of D 1 and D 2

Figure 3 :

 3 Figure 3: (left and middle) The continuous synthetic data experiment from Appendix F.1 with M 1 = M 2 = 20, c = c 1 = c 2 = 0.1 and L n = L s = 1000. The number of factors:(left) K 1 = K 2 = K = 1 and (middle) K 1 = K 2 = K = 10. (right):An experimental analysis of the performance of DCCAg with generalized covariance matrices using different parameters δ j for the processing points. The numbers in the legend correspond to the values of δ defining δ j via (46) in Appendix D.3. The default value (def) is δ = 0.1. The data is the discrete synthetic data as described in Section 5 with the parameters set as in Fig.2(right).

Table 1 :

 1 .1. Factor loadings (a.k.a. topics) extracted from the Hansard collection for K = 20 with DCCA.

	nato	otan	work	travail	board	commission	nisga	nisga
	kosovo	kosovo	workers	négociations	wheat	blé	treaty	autochtones
	forces	militaires	strike	travailleurs	farmers	agriculteurs	aboriginal	traité
	military	guerre	legislation	grève	grain	administration	agreement	accord
	war	international	union	emploi	producers	producteurs	right	droit
	troops	pays	agreement	droit	amendment	grain	land	nations
	country	réfugiés	labour	syndicat	market	conseil	reserve	britannique
	world	situation	right	services	directors	ouest	national	indiennes
	national	paix	services	accord	western	amendement	british	terre
	peace	yougoslavie	negotiations	voix	election	comité	columbia	colombie

Note that[START_REF] Podosinnikova | Rethinking LDA: Moment matching for discrete ICA[END_REF] show that DICA is closely connected (and under some conditions is equivalent) to latent Dirichlet allocation(Blei et al.,

2003). Due to the close relation of DCCA and DICA, the former is thus closely related to the multi-view topic models (see, e.g.,Blei & Jordan, 2003).

Note that S 21 := cov(x 2 , x 1 ) is just the transpose of S 12 .

We find the name "generalized covariance matrix" to be more meaningful than "charrelation" matrix as was proposed by previous authors (see, e.g.Slapak & Yeredor, 2012a,b).

Note that when the diagonal form has terms diag[e t ], we simply multiply the expression by diag[e -t ].

A real matrix A is normal if A A = AA .

An orthogonal joint diagonalization problem corresponds to the case where the matrices A 1 , A 2 , . . . , A N are normal and, hence, diagonalizable by an orthogonal matrix Q.

In the JUST algorithm[START_REF] Iferroudjene | A new Jacobi-like method for joint diagonalization of arbitrary non-defective matrices[END_REF], this objective function is also considered for the (shear transformed) matrix A ( ) .

The JUST algorithm is an exception here, since it minimizes the diagonality measure Off[A ( ) ] of the (shear transformed) matrices A ( ) with respect to y.

The code is available at: http://perso.telecom-paristech.fr/ cardoso/Algo/Jade/jadeR.m

This list of words was obtained by looking at words that appear in the top-20 words of a large number of topics in a first experiment. Removing these words did not change much the content of the topics, but made them much more interpretable.

end for 12: end for 13: Output: Q ( ) Algorithm. Non-orthogonal Jacobi-like joint diagonalization algorithms have the high level structure which is outlined in Alg. 1.

The algorithm iteratively constructs the sequence of matrices A ( ) = A = U ( ) (θ * ) for the chosen in accordance with some rules (see below) optimal shear parameter y * and optimal Jacobi (=Givens) angle θ * .

For the theoretical analysis purposes, the two transforms are considered separately:

A ( ) n = S ( )-1 (y)A ( -1) n S ( ) (y),

Each such iteration is a combination of the iteration k and the pivots p and q (see Alg. 1). The iteration k is referred to as a sweep. Within each sweep k, M (M -1)/2 pivots p < q are chosen in accordance with the lexicographical rule. The rule for the choice of pivots can affect convergence as was analyzed for the single matrix case (see, e.g., [START_REF] Ruhe | On the quadratic convergene of a generalization of the Jacobi method to arbitrary matrices[END_REF][START_REF] Eberlein | A Jacobi-like method for the automatic computation of eigenvalues and eigenvectors of an arbitrary matrix[END_REF], where more sophisticated rules were proposed for the algorithm to have a quadratic convergence phase. However, up to our best knowledge, no such analysis was done for the several matrices case. We assume the simple lexicographical rule all over the paper.

The shear transform is defined by the hyperbolic rotation matrix S ( ) = S ( ) (y) which is equal to the identity matrix except for the following entries

where the shear parameter y ∈ R. The unitary transform is defined by the Jacobi (=Givens) rotation matrix U ( ) = U ( ) (θ) which is equal to the identity matrix except for the following entries

where the Jacobi (=Givens) angle θ ∈ -π 4 , π 4 . The following two objective functions are of the central importance for this type of algorithms: (a) the sum of squares of all the off-diagonal elements of the matrices 9 A ( ) which are the transformed with the unitary transform U ( ) matrices A ( ) :

and (b) the sum of the squared Frobenius norms of the matrices A ( ) which are the transformed with the share transform S ( ) matrices A ( -1) :

We refer to (51) as the diagonality measure and to (52) as the normality measure.

All the considered algorithms find the optimal Jacobi angle θ * as the minimizer of the diagonality measure of the (unitary transformed) matrices A ( ) (48):

which admits a unique closed form solution [START_REF] Cardoso | Jacobi angles for simultaneous diagonalization[END_REF]. The optimal shear parameter y * is found 10 as a minimizer of the normality measure of the (shear transformed) matrices A ( ) (48):

All the considered algorithms [START_REF] Fu | Simultaneous diagonalization with similarity transformation for non-defective matrices[END_REF][START_REF] Iferroudjene | A new Jacobi-like method for joint diagonalization of arbitrary non-defective matrices[END_REF][START_REF] Luciani | Joint eigenvalue decomposition using polar matrix factorization[END_REF] solve this step only approximately. In particular, the sh-rt algorithm [START_REF] Fu | Simultaneous diagonalization with similarity transformation for non-defective matrices[END_REF] approximates the equation for finding the nulls of the gradient of the objective; the JUST algorithm [START_REF] Iferroudjene | A new Jacobi-like method for joint diagonalization of arbitrary non-defective matrices[END_REF] replaces the normality measure with the diagonality measure and provides a closed form solution for the resulting problem; and the JDTM algorithm [START_REF] Luciani | Joint eigenvalue decomposition using polar matrix factorization[END_REF] replaces the normality measure with the sum of only two squared elements A n,pq and A n,qp and provides a closed form solution for the resulting problem. After stemming and removing stop words, several files had different number of documents in each language and had to be removed too. The numbers of these files are: 16, 36, 49 55, 88, 103, 110, 114, 123, 155, 159, 204, 229, 240, 2-17, 2-35. We also removed the 15 most frequent words from each language. These include:

-in English: Mr, govern, member, speaker, minist(er), Hon, Canadian, Canada, bill, hous(e), peopl(e), year, act, motion, question;

-in French: gouvern(er), président, loi, déput(é), ministr(e), canadien, Canada, projet, Monsieur, question, part(y), chambr(e), premi(er), motion, Hon.

Removing these words is not necessary, but improves the presentation of the learned topics significantly. Indeed, the most frequent words tend to appear in nearly every topic (often in pairs in both languages as translations of each other, e.g., "member" and "député" or "Canada" in both languages, which confirms one more time the correctness of our algorithm).

Finally, we select M 1 = M 2 = 5, 000 words for each language to form matrices X 1 and X 2 each containing N = 11, 969 documents in columns. As stemming removes the words endings, we map the stemmed words to the respective most frequent original words when showing off the topics in Tables 23456.

Supplementary References


T Detroux
email: tdetroux@ulg.ac.be

L Renson
email: l.renson@ulg.ac.be

L Masset
email: luc.masset@ulg.ac.be

G Kerschen
email: g.kerschen@ulg.ac.be

Presqu'île de Giens (Var) Bifurcation Analysis of Large-Scale Dynamical Systems Using the Harmonic Balance Method

Keywords: continuation of periodic solutions, bifurcation tracking, harmonic balance method, quasiperiodic oscillations, detached resonance curves

The harmonic balance (HB) method has been widely used in the past few years, as a numerical tool for the study of nonlinear industrial models. However, in its classical formulation the HB method is limited to the approximation of periodic solutions. The present paper proposes to extend the method to the detection and tracking of codimension-1 bifurcations in the system parameters space. As an application, the frequency response of a spacecraft and its bifurcations are studied, together with two nonlinear phenomena, namely quasiperiodic oscillations and detached resonance curves.

Introduction

Because nowadays engineering structures are designed to be lighter and operate in more severe conditions, nonlinear phenomena such as amplitude jumps, modal interactions, limit cycle oscillations and quasiperiodic (QP) oscillations are expected to occur [START_REF] Nayfeh | Nonlinear oscillations[END_REF]. For most of these nonlinear systems, bifurcations of periodic solutions play a key role in the response dynamics; for example, fold bifurcations indicate a stability change for the solutions, while QP oscillations are encountered in the vicinity of Neimark-Sacker (NS) bifurcations. In that regard, it seems relevant to include stability analysis and bifurcation monitoring while performing parametric studies of the structures.

Time-domain methods such as shooting technique and orthogonal collocation, which deal with the resolution of a boundary value problem, usually prove accurate to analyze periodic solutions and bifurcations of low-dimensional structures. When applied to larger systems however, their computational burden become substantial. As an efficient alternative, most engineers apply a frequency-domain method, the so-called harmonic balance (HB) method (see, e.g., [START_REF] Jaumouillé | An adaptive harmonic balance method for predicting the nonlinear dynamic responses of mechanical systems-Application to bolted structures[END_REF]). Indeed, through the approximation of the solutions with truncated Fourier series, the HB method involves algebraic equations with usually less unknowns than time-domain techniques. Nevertheless, in spite of its performance and accuracy, the HB method has never been extended to track bifurcations of mechanical structures. This is why the present paper proposes an extension of the methodology to study bifurcations in the system parameters space.

The first part of this paper is devoted to the HB theory and its formulation in the framework of a continuation algorithm for bifurcation tracking. In section 3, the methodology is validated with the analysis of the nonlinear dynamics of an industrial spacecraft.

2 Harmonic balance method and bifurcation tracking

Periodic solutions

The method will be applied to general non-autonomous nonlinear dynamical systems with n degreesof-freedom (DOFs) whose equations of motion are

Mẍ + Cẋ + Kx = f ext (ω,t) -f nl (x, ẋ) = f(x, ẋ, ω,t) (1)
where M,C and K are the mass, damping and stiffness matrices respectively, x represents the displacements, the dots refer to the derivatives with respect to time t, f nl represents the nonlinear forces and f ext stands for the periodic external forces (harmonic excitation, for example) with frequency ω. The term f gathers both the external and nonlinear forces.

The HB procedure consists in approximating x and f with their Fourier coefficients z and b, respectively, and applying a Galerkin projection. Equations (1) then become

h(z, ω) ≡ A(ω)z -b(z, ω) = 0 (2)
where A is the matrix describing the linear dynamics of the system (see [START_REF] Jaumouillé | An adaptive harmonic balance method for predicting the nonlinear dynamic responses of mechanical systems-Application to bolted structures[END_REF] for further details). When embedded in a continuation procedure, solutions of ( 2) can be tracked with respect to the frequency ω, which gives a frequency response of the system.

Stability analysis and detection of bifurcations

In the frequency domain, the stability of a periodic solution is usually assessed with Hill's method, which provides approximation of the Floquet exponents. In [START_REF] Detroux | The harmonic balance method for advanced analysis and design of nonlinear mechanical systems[END_REF], it is shown that these exponents, here denoted λ, can be found among the eigenvalues of B, a matrix whose components are by-products of the HB method, i.e., depend only on M, C, K and J z , the jacobian of (2) with respect to z. If at least one of the Floquet exponents has a positive real part, then the solution is unstable, otherwise it is asymptotically stable.

In this work, the detection of fold (F) and NS bifurcations along the frequency response is sought. This is achieved through the evaluation of test functions at each iteration of the branch, which have the property to pass through zero at a bifurcation. From the Floquet theory, it is known that a fold bifurcation is detected when one of the Floquet exponents λ crosses the imaginary axis through the real axis. As a consequence, the following test function has a root at a fold bifurcation:

ϕ F = B (3)
The NS bifurcations are characterized by a pair of Floquet exponents crossing the imaginary axis as complex conjugates. A dedicated test function is thus

ϕ NS = B⊙ ( 4
)
where B⊙ stands for the bialternate matrix product of B, which has the property to be singular when B has two complex-conjugate eigenvalues with no real part.

Tracking of bifurcations

In order to track bifurcations in a codimension-2 parameter space, one has to append one equation to [START_REF] Govaerts | Numerical methods for bifurcations of dynamical equilibria[END_REF], describing the bifurcation of interest:

{ h = 0 ϕ = 0 (5)
Instead of considering ϕ = ϕ F and ϕ = ϕ NS , which can lead to scaling problems for large matrices, in this paper the use of the so-called bordering technique [START_REF] Govaerts | Numerical methods for bifurcations of dynamical equilibria[END_REF] is proposed. The idea behind this technique is to replace the evaluation of the determinant of a given matrix G by the evaluation of a scalar function, herein denoted g, which vanishes as regular zero for the same system state and parameters as the determinant. It can be shown that a candidate for g is obtained by solving the bordered system

[ G p q * 0 ] [ w g ] = [ 0 1 ] (6)
In this system, * denotes a conjugate transpose, and vectors p and q are chosen to ensure the nonsingularity of the bordered matrix. From (3-4), one then proposes the following equations for the extended system (5):

Fold bifurcations: ϕ = g F with G F = B NS bifurcations: ϕ = g NS with G NS = B⊙ (7) 2
3 Validation of the method on the study of an industrial, complex model with strong nonlinearities: the SmallSat

The example studied is referred to as the SmallSat, a 64-kg structure represented in Figure 1(a) and which was conceived by EADS-Astrium as a platform for small satellites. As depicted in Figure 1(b), a support bracket connects to one of the eight walls the so-called wheel elastomer mounting system (WEMS) device which is loaded with an 8-kg dummy inertia wheel. The WEMS device is a mechanical filter which mitigates disturbances coming from the inertia wheel through the presence of a soft elastomeric interface between its mobile part, i.e. the inertia wheel and a supporting metallic cross, and its fixed part, i.e. the bracket and by extension the spacecraft. Moreover, eight mechanical stops limit the axial and lateral motions of the WEMS mobile part during launch, which gives rise to strongly nonlinear dynamical phenomena.

A reduced finite element model of the SmallSat with 37 DOFs is developed to conduct numerical experiments. Proportional damping is considered and the high dissipation in the elastomer components of the WEMS is described using lumped dashpots with coefficients c ax and c lat for axial (vertical) and lateral directions, respectively. Bilinear and trilinear springs are finally introduced to model the nonlinearities of the connections between the WEMS and the rest of the SmallSat. The DOF studied in this paper, referred to as NC1-Z, represents the vertical displacements of one of the nodes involved in the nonlinear connections.

Figure 1 :

 1 Figure 1: SmallSat spacecraft. (a) Photograph. (b) Schematic of the WEMS.

Figure 2 (Figure 2 :

 22 Figure 2: Frequency responses (black lines) and bifurcation curves (orange and blue lines) of the NC1-Z node of the SmallSat, for harmonic excitations applied vertically to the inertia wheel. (a) Tracking of fold bifurcations in the forcing amplitude F -frequency ω parameters space, for a fixed value of axial damping c ax = 63 Ns/m. (b) Tracking of NS bifurcations in the axial damping c ax -frequency ω parameters space, for a fixed value of forcing amplitude F = 155 N.