# An Effective Iterative Solution for Independent Vector Analysis with Convergence Guarantees

Clément Cosserat, Ben Gabrielson, Emilie Chouzenoux, Jean-Christophe

Pesquet, Tülay Adali

► To cite this version:

Clément Cosserat, Ben Gabrielson, Emilie Chouzenoux, Jean-Christophe Pesquet, Tülay Adali. An Effective Iterative Solution for Independent Vector Analysis with Convergence Guarantees. IEEE Transactions on Signal Processing, 2025, 73, pp.4718 - 4733. 10.1109/TSP.2025.3620539. hal-04785828v2

HAL Id: hal-04785828

https://hal.science/hal-04785828v2

Submitted on 13 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# An Effective Iterative Solution for Independent Vector Analysis with Convergence Guarantees

Clément Cosserat, Student Member IEEE, Ben Gabrielson, Student Member, IEEE, Emilie Chouzenoux, Senior Member, IEEE, Jean-Christophe Pesquet, Fellow, IEEE, and Tülay Adali, Fellow IEEE

**Abstract**—Independent vector analysis (IVA) is an attractive solution to address the problem of joint blind source separation (JBSS), that is, the simultaneous extraction of latent sources from several datasets implicitly sharing some information. Among IVA approaches, we focus here on the celebrated IVA-G model, that describes observed data through the mixing of independent Gaussian source vectors across the datasets. IVA-G algorithms usually seek the values of demixing matrices that maximize the joint likelihood of the datasets, estimating the sources using these demixing matrices. Instead, we write the likelihood of the data with respect to both the demixing matrices and the precision matrices of the source estimates. This allows us to formulate a cost function whose mathematical properties enable the use of a proximal alternating algorithm based on closed form operators with provable convergence to a critical point. After establishing the convergence properties of the new algorithm, we illustrate its desirable performance in separating sources with covariance structures that represent varying degrees of difficulty for JBSS.

**Index Terms**—IVA, PALM Algorithm, Maximum-Likelihood, Blind Source Separation, Proximal methods

## I. INTRODUCTION

THE *blind source separation* problem (BSS) aims at factorizing a data matrix as a product of a mixing matrix and a source data matrix. BSS is thus a data-driven method to extract latent features from a dataset, which have a broad variety of uses. It offers a wide range of applications in signal processing and engineering including neuroimaging data analysis [1], communications [2], [3], remote sensing [4], to name a few [5]. The latent sources are interpreted as physical quantities of interest that cannot be measured directly, e.g. brain connectivity networks [1], [6], [7], or can be used as features for further tasks such as classification [8]. The BSS problem generalizes to the *joint blind source separation* problem (JBSS) when multiple datasets are analyzed jointly to benefit from their shared information. JBSS is necessary to fully analyze datasets that share similarities. This case presents itself when the datasets contain measures of a same phenomenon for different subjects [6], [9], different measurement methods [10], or more generally for various modalities [11]. Allowing the interaction between the datasets leads JBSS to achieve more accurate separation than multiple separate BSS in general [6], [12].

One way to address BSS is to model the rows of the source dataset as samples of mutually independent random variables or random processes called sources and the datasets as mixtures of independent sources [13]. This approach is known as *independent component analysis* (ICA) [14], which has been one of the most popular ways to achieve BSS due to its uniqueness guarantees under very general conditions. A couple of decades ago, [15] and [16] introduced statistical links between the datasets by considering *source component*

vectors (SCVs) whose components account respectively for each dataset. This addition to the ICA model was motivated by the will to exploit dependencies between the frequencies in the frequency domain BSS problem [15], [16], [17]. Then, instead of recovering the independent components for each frequency bin, the authors modeled the data as mixtures of independent SCVs, thus generalizing ICA to *independent vector analysis* (IVA) [15], [18]. Within each SCV, the statistical dependence (correlation for IVA-G) across the datasets is taken into account. IVA methods gather a family of algorithms that aim at recovering the SCVs and mixing coefficients that generated the observed data. The JBSS problem can be addressed with many other methods such as groupICA [7], [19], or joint ICA [20], but IVA is more powerful as it offers a greater flexibility and helps preserve the individual variability represented by each dataset, e.g., in multisubject analyses [21], [22]. Moreover, it enables a common factorization of the datasets without needing to realign the sources *a posteriori*, which can be costly, thus alleviating the inherent permutation ambiguity across the datasets [16].

IVA is typically formulated using a *maximum likelihood* (ML) estimator, which means that the estimation is done by solving an optimization problem where the cost-function is derived from the log-likelihood of the data [18]. In [12], [15], [16] we see IVA presented as a minimization of the mutual information of the SCVs to maximize their independence, those two formulations are equivalent when the number of samples tends to infinity [6].

In this paper, we focus on the case of non-degenerate, centered Gaussian SCVs, placing ourselves in the so-called IVA-G framework, as presented in [12]. This model is convenient, as the SCVs are fully described by their covariance matrix (or equivalently by their precision matrix). Besides, since IVA-G algorithms only make use of *second order statistics* (SOS) of the data, this makes them the simplest and most efficient among the IVA algorithms. Previous works have additionally shown that even when using a multivariate pdf model that takes into account higher-order statistics (HOS), i.e., higher than second order, accounting for SOS remains crucial by estimating a scatter matrix [23], [18]. Since SOS can be estimated more efficiently than HOS, and linear dependence is a powerful statistic, this can help provide better alignment of source estimates across the datasets compared with the case when second-order statistics are not accounted for in the SCV pdf model, as shown in [23], in the context of an application to fMRI analysis.

In the Gaussian case, having the demixing matrices gives a natural estimate of the covariance matrices of the SCVs through their empirical covariance, this is why the cost function in [12] only depends on the demixing matrices. In this work however, we choose to write the ML such that it
explicitly takes the SCVs probability density function (pdf) as an input variable, to benefit from analytical properties of the resulting cost.

So far, most IVA-G algorithms are based on standard optimization methods, like Newton's method or gradient descent, without demonstration of explicit convergence guarantees. The idea of applying majorization-minimization (MM) methods to IVA has been explored in [24], [25], [26], [27] or [28], which suggest algorithms with the valuable property of ensuring a provable decrease of the cost function. In this article, we show that our proposed cost function to jointly search for the Gaussian SCVs and demixing matrices generalizes the cost function in [12], in the sense that replacing the precision matrices of the SCVs by the inverse of their empirical covariance matrices for the given demixing matrices, we recover a cost function that only depends on those demixing matrices and that is equal to the one in [12] up to a constant. The choice of a cost function that jointly acts on the sought demixing and covariance matrices offers an attractive structure enabling the explicit incorporation of priors and constraints, and the use of mathematically sound MM algorithms. Hence, we design a Proximal Alternating Linearized Minimization algorithm (PALM) [29], dedicated to the IVA-G problem, and show that it converges to a critical point under some mild assumptions. The resulting scheme, denoted PALM-IVA-G, is also shown to be fast, and to establish a reliable estimation performance in practice, and at least as good as the state-of-the-art methods for IVA-G.

Our contributions can be summarized as follows. First, we provide a novel variational formulation for IVA-G, introducing a new cost function and establishing its mathematical properties. Second, building upon these properties, we design the PALM-IVA-G algorithm, to solve the resulting minimization problem, and show its convergence under mild assumptions. Third, we numerically illustrate the desirable performance of our approach, in comparison with state-of-the-art approaches, on a number of scenarios that cover various degrees of dependence across the datasets.

The paper is organized as follows. Section II introduces the JBSS problem, the IVA-G framework to address it, and presents our cost function based on maximum likelihood estimation. Section III then focuses on our proposition of an original iterative algorithm to solve the IVA-G problem. Convergence results are presented in Section IV. Numerical experiments assessing the validity and good performance of our approach are presented in Section V. Section VI concludes the paper.

## II. THE IVA-G PROBLEM

### A. Notation

Throughout the paper, we use bold upper case symbols for matrices, bold lower case symbols for column or row vectors, calligraphic upper case symbols for tensors of order three or more, and regular lower case symbols for scalars. Italic font is used for deterministic quantities, while the regular one is used for random quantities. For instance, let $\mathbf{u} = [u_1, \dots, u_N]^\top$ be an $\mathbb{R}^N$-valued random vector, for which we draw $\mathcal{V}$ realizations that we stack, columnwise, in a matrix $\mathbf{U} \in \mathbb{R}^{N \times V}$. The latter can also be written as $\mathbf{U} = [\mathbf{u}_1^\top, \dots, \mathbf{u}_N^\top]^\top$ with, for every $n \in \{1, \dots, N\}$, $\mathbf{u}_n \in \mathbb{R}^V$ a column vector of $\mathcal{V}$ realizations of the scalar

random variable $u_n$.

The subset of non-singular (resp. symmetric) matrices of $\mathbb{R}^{N \times N}$ is denoted $\mathcal{GL}_N(\mathbb{R})$ (resp. $\mathcal{S}_N$). $\mathbf{I}_N$ is the $N \times N$ identity matrix. We also denote $\mathbf{S}_N^+ = \text{the set of positive semi-definite matrices of size } N$, and $\mathbf{S}_N^{++} = \mathbf{S}_N^+ \cap \mathcal{GL}_N(\mathbb{R})$ the set of positive definite (PD) matrices. $\| \cdot \|$ denotes the Frobenius norm for any vector, matrix, and tensor of order three or more.

For any matrix $\mathbf{A}$, we denote by $\|\mathbf{A}\|_S$ its spectral norm, equal to the largest singular value of $\mathbf{A}$. For every $n \in \{1, \dots, N\}$, $\mathbf{a}_n^\top$ denotes the $n$-th row of $\mathbf{A}$. If $\mathbf{A}$ is a square matrix, we note $\sigma_{\mathbf{A},l} = (\sigma_{\mathbf{A},l})_{1 \le l \le N} \in \mathbb{R}^N$ the vector of its singular values, tr($\mathbf{A}$) its trace, and det($\mathbf{A}$) its determinant. diag($\mathbf{A}$) $\in \mathbb{R}^N$ is the vector whose entries are the diagonal coefficients of $\mathbf{A}$, whereas for $\mathbf{a} \in \mathbb{R}^N$, Diag($\mathbf{a}$) (with an upper case) is the diagonal matrix whose coefficients are the components of $\mathbf{a}$. We also note Diag($\mathbf{A}$) the diagonal matrix whose diagonal coefficients are the same as those of $\mathbf{A}$. Finally, for any $(N, M) \in (\mathbb{N}\setminus\{0\})^2$, we consider the matricial scalar product defined as $(\forall (\mathbf{A}, \mathbf{B}) \in (\mathbb{R}^{N \times M})^2) \langle \mathbf{A} | \mathbf{B} \rangle = \text{tr}(\mathbf{A}^\top \mathbf{B})$.

### B. JBSS problem

Let $K, N, V$ be positive integers, we consider $K$ datasets denoted $(\mathbf{X}^{[k]})_{1 \le k \le K}$ where $\forall k \in \{1, \dots, K\}$, $\mathbf{X}^{[k]} \in \mathbb{R}^{N \times V}$ is a real-valued matrix. For instance, in an fMRI analysis problem, for $n \in \{1, \dots, N\}$, and $k \in \{1, \dots, K\}$, the $n$-th row of $\mathbf{X}^{[k]}$, $(\mathbf{x}_n^{[k]})^\top$ could model the blood-oxygen-level-dependent (BOLD) contrast in the $V$ voxels of a 3D brain model, measured at the $n$-th acquisition time for the $k$-th subject, within a cohort of $K$ patients [1], [7]. In the context of speech separation, it would typically be the Fourier coefficient corresponding to the $k$-th frequency bin of the signal recorded by the $n$-th microphone over $V$ time frames spread over the duration of the signal [16], [17]. We assume that the observed datasets are obtained by a linear mixing of latent source datasets, i.e., $(\forall k \in \{1, \dots, K\})$:

$$ \mathbf{X}^{[k]} = \mathbf{A}^{[k]} \mathbf{S}^{[k]} \in \mathbb{R}^{N \times V}, \quad (1) $$

where $\mathbf{A}^{[k]} = (a_{m,n}^{[k]})_{1 \le m \le N, 1 \le n \le N} \in \mathbb{R}^{N \times N}$ is a square non-singular mixing matrix and $\mathbf{S}^{[k]} \in \mathbb{R}^{N \times V}$ is a latent matrix whose coefficients are typically interpreted as hidden features of the phenomenon described by $\mathbf{X}^{[k]}$. The JBSS problem consists in simultaneously estimating for all $k \in \{1, \dots, K\}$, both $\mathbf{A}^{[k]}$ and $\mathbf{S}^{[k]}$ from $\mathbf{X}^{[k]}$. We estimate the inverse of the mixing matrices by the so-called demixing matrices $(\mathbf{W}^{[k]})_{1 \le k \le K}$, and deduce estimates $(\mathbf{Y}^{[k]})_{1 \le k \le K}$ for the source datasets $(\mathbf{S}^{[k]})_{1 \le k \le K}$ by calculating

$$ (\forall k \in \{1, \dots, K\}) \quad \mathbf{Y}^{[k]} = \mathbf{W}^{[k]} \mathbf{X}^{[k]}. \quad (2) $$

In a nutshell, using tensor notations, JBSS amounts to providing an estimate $\mathcal{Y} = [\mathbf{Y}^{[1]}, \dots, \mathbf{Y}^{[K]}]$ of the source tensor $\mathcal{S} = [\mathbf{S}^{[1]}, \dots, \mathbf{S}^{[K]}]$ via the demixing tensor $\mathcal{W} = [\mathbf{W}^{[1]}, \dots, \mathbf{W}^{[K]}]$, which is an estimate of the slice-wise inverse of the mixing tensor $\mathcal{A} = [\mathbf{A}^{[1]}, \dots, \mathbf{A}^{[K]}]$, given the datasets $\mathcal{X} = [\mathbf{X}^{[1]}, \dots, \mathbf{X}^{[K]}]$. Here, $\mathcal{S}, \mathcal{X}$ and $\mathcal{Y} \in \mathbb{R}^{N \times V \times K}$ and $\mathcal{A}$ and $\mathcal{W} \in \mathbb{R}^{N \times N \times K}$.
### C. IVA model

Assuming that the datasets share some information that can be leveraged to separate the sources more accurately than if dealt with individually, IVA [30] models this interdependence of the $X^{[k]}$ through statistical links between the latent datasets. More precisely, the rows of the latent datasets with the same index are assumed to be correlated, while being independent from rows with different indices. To formalize this, IVA models the columns of the $X^{[k]}$ (resp. $S^{[k]}$) as independent samples from $\mathbb{R}^N$-valued random vectors $x^{[k]} = [x_1^{[k]}, \dots, x_N^{[k]}]^{\top}$ (resp. $s^{[k]} = [s_1^{[k]}, \dots, s_N^{[k]}]^{\top}$). We can thus rewrite the model in (1) using random vector notations: $(\forall k \in \{1, \dots, K\}), x^{[k]} = A^{[k]}s^{[k]}$. Regrouping the components of the $s^{[k]}$ with corresponding indices, we obtain $N \mathbb{R}^K$-valued random vectors, $s_n = [s_n^{[1]}, \dots, s_n^{[K]}]^{\top}$ for $n \in \{1, \dots, N\}$ called source component vectors, where each entry of an SCV accounts for the corresponding dataset. The goal of IVA is to make the SCVs as independent as possible.

IVA aims at recovering the sources by building demixing matrices $(\mathbf{W}^{[k]})_{1 \le k \le K}$. At the same time, it builds estimated SCVs $(\mathbf{y}_n)_{1 \le n \le N}$ whose distributions have probability density functions respectively denoted $(p_{\mathbf{y}_n})_{1 \le n \le N}$ and that we suppose mutually independent. Similarly as for SCVs, we denote $(\forall n \in \{1, \dots, N\}), \mathbf{y}_n = [y_n^{[1]}, \dots, y_n^{[K]}]^{\top}$ the estimated SCVs, and we reorganize the components to define $(\forall k \in \{1, \dots, K\}), \mathbf{y}^{[k]} = [y_1^{[k]}, \dots, y_N^{[k]}]^{\top}$. Then, we see with (2) that for all $k \in \{1, \dots, K\}, \mathbf{x}^{[k]}$ is estimated by $(\mathbf{W}^{[k]})^{-1}\mathbf{y}^{[k]}$, whose probability distribution is entirely determined by the demixing matrices and the $(p_{\mathbf{y}_n})_{1 \le n \le N}$. Said otherwise, $\mathcal{W}$ and $(p_{\mathbf{y}_n})_{1 \le n \le N}$ yield an estimated generative model for our observed datasets.

### D. IVA-G cost function

In the following, for every $n \in \{1, \dots, N\}$, we denote by $\mathbf{Y}_n \in \mathbb{R}^{K \times V}$ the matrix obtained by stacking vertically the $n$-th rows of $\mathbf{Y}^{[k]}$ for $k \in \{1, \dots, K\}$. Therefore, using (2), we have

$$ (\forall n \in \{1, \dots, N\}) \quad \mathbf{Y}_n = \mathbf{W}_n \mathbf{X}, \quad (3) $$

where

$$ \mathbf{W}_n = \begin{pmatrix} \mathbf{w}_n^{[1]\top} & 0 & \dots & 0 \\ 0 & \mathbf{w}_n^{[2]\top} & \dots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \dots & \mathbf{w}_n^{[K]\top} \end{pmatrix} \in \mathbb{R}^{K \times KN} \quad (4) $$

and

$$ \mathbf{X} = \begin{pmatrix} \mathbf{X}^{[1]} \\ \vdots \\ \mathbf{X}^{[K]} \end{pmatrix} \in \mathbb{R}^{KN \times V}. \quad (5) $$

The objective of IVA is to determine $\mathcal{W}$ and $(p_{\mathbf{y}_n})_{1 \le n \le N}$ that maximize the log-likelihood of $\mathbf{X}$ in our estimated generative model, given by [18]

$$ \sum_{n=1}^{N} \sum_{v=1}^{V} \log p_{\mathbf{y}_n}(\mathbf{y}_{n,v}) + V \sum_{k=1}^{K} \log |\det \mathbf{W}^{[k]}| \quad (6) $$

where $\mathbf{y}_{n,v}$ is the $v$-th column of $\mathbf{Y}_n$. In the Gaussian case, considered here, we model $(\mathbf{y}_n)_{1 \le n \le N}$ as centered

non-degenerate Gaussian vectors, whose pdf is thus entirely determined by their covariance matrices, or equivalently, their precision matrices that we denote $(C_n)_{1 \le n \le N}$. As done previously, to simplify our notation, we introduce the tensor $\mathcal{C} = [C_1, \dots, C_N] \in (S_K^{++})^N$ that gathers the estimated precision matrices of all the SCVs. Moreover, we assume sample independence, so the pdf of $\mathbf{Y}_n$ is the product of the pdfs of its columns. Under this Gaussian modeling, we have,

$$ \begin{gathered} (\forall \mathbf{y} \in \mathbb{R}^K) \\ \log p_{\mathbf{y}_n}(\mathbf{y}) = \frac{1}{2} \log \det \mathcal{C}_n - \frac{K}{2} \log 2\pi - \frac{1}{2} \mathbf{y}^{\top} \mathcal{C}_n \mathbf{y}. \quad (7) \end{gathered} $$

Hence, the problem becomes equivalent to minimizing, with respect to $(\mathcal{W}, \mathcal{C})$, the cost function $J_{\text{IVA-G}}(\mathcal{W}, \mathcal{C})$ defined as

$$ \begin{align} & J_{\text{IVA-G}}(\mathcal{W}, \mathcal{C}) + \frac{1}{2} \sum_{n=1}^{N} \log \det \mathcal{C}_n + \sum_{k=1}^{K} \log |\det \mathbf{W}^{[k]}| \nonumber \\ &= \frac{1}{2} \sum_{n=1}^{N} \frac{1}{V} \operatorname{tr}(\mathcal{C}_n \mathbf{Y}_n \mathbf{Y}_n^{\top}) \nonumber \\ &= \frac{1}{2} \sum_{n=1}^{N} \operatorname{tr}(\mathcal{C}_n \mathbf{W}_n \hat{\mathbf{R}}_x \mathbf{W}_n^{\top}), \tag{8} \end{align} $$

where $\hat{\mathbf{R}}_x = \frac{1}{V} \mathbf{X} \mathbf{X}^{\top}$ is the empirical covariance matrix of $\hat{\mathbf{x}}$.

Note that the domain of the function (8) can be extended to $\mathbb{R}^{N \times N \times K} \times \mathbb{R}^{K \times K \times N}$ by setting $J_{\text{IVA-G}}(\mathcal{W}, \mathcal{C}) = +\infty$ if $\mathbf{W}^{[k]}$ is singular for some $k \in \{1, \dots, K\}$ or if $\mathcal{C}_n$ is not symmetric and positive definite for some $n \in \{1, \dots, N\}$.

**Remark 1** As we will show in Section IV, the non-singularity of $\hat{\mathbf{R}}_x$ is a sufficient condition for the lower-boundedness of (8), and as such, for the well-posedness of our optimization problem. In practice, for a large number of samples — that is $V > KN$ — drawn from a continuous probability distribution, the empirical covariance matrix $\hat{\mathbf{R}}_x$ is non-singular almost surely. Hence, we will suppose that this condition holds in the remainder of this article. This condition is usually verified for fMRI or EEG data [19], [31], [32] where we typically have $V \approx 50000$, $K \approx 100$ and $N \approx 50$. In other applications such as speech separation, it may not be the case, see [16] for instance. A typical workaround is to replace $\hat{\mathbf{R}}_x$ by $\hat{\mathbf{R}}'_x = \hat{\mathbf{R}}_x + a\mathbf{I}_{NK}$, for some $a > 0$. While this has a negligible effect when $\hat{\mathbf{R}}_x$ is already non-singular as we show in Appendix G-B, this ensures $\hat{\mathbf{R}}_x$ is non-singular and the optimization problem is well-posed. However, IVA-G's performance will still be negatively affected if $V$ is too low, as we show in Appendix G-A. In these cases, specific application-dependent regularization approaches may be appropriate.

### E. IVA-G-V and IVA-G-N approaches

In [12], the authors proposed two algorithms, called IVA-G-V and IVA-G-N, to solve the IVA-G estimation problem. To do so, they reformulate the problem as the minimization of the following cost function [12]:

$$ \tilde{J}_{\text{IVA-G}}(\mathcal{W}) = \frac{1}{2} \sum_{n=1}^{N} \log \det(\mathbf{W}_n \hat{\mathbf{R}}_x \mathbf{W}_n^{\top}) - \sum_{k=1}^{K} \log |\det \mathbf{W}^{[k]}|. \quad (9) $$
These approaches implicitly estimate the covariance matrices of the sources by $C_n^{-1} = W_n \hat{R}_x W_n^\top$, that is, the empirical covariance matrices of the $(y_n)_{1 \le n \le N}$. The minimization of (9) is usually performed using either a gradient-based, or a Newton-based solver, leading to IVA-G-V and IVA-G-N schemes, respectively. To the best of our knowledge, these algorithms do not benefit from strong guarantees in terms of convergence of the iterates.

It is easy to show that finding $(\mathcal{W}, \mathcal{C})$ that minimizes the proposed objective function (8) is actually equivalent to finding $\mathcal{W}$ that minimizes the cost (9) used in state-of-the-art IVA-G approaches, as stated in the following known result.

**Theorem 1 (Equivalence between $J_{IVA-G}$ and $\tilde{J}_{IVA-G}$)**

For a given $\mathcal{W} \in \mathbb{R}^{N \times N \times K}$, if all the $\mathcal{W}^{[k]}$ are non-singular, then

(i) $J_{IVA-G}(\mathcal{W}, \cdot)$ is minimized over $\mathbb{R}^{K \times K \times N}$ at $\hat{\mathcal{C}}(\mathcal{W}) = [\hat{C}_1(W_1), \dots, \hat{C}_N(W_N)]$ where

$$ (\forall n \in \{1, \dots, N\}) \quad \hat{C}_n(W_n) = (W_n \hat{R}_x W_n^\top)^{-1}. \quad (10) $$

(ii) We have, for every $\mathcal{W} \in \mathbb{R}^{N \times N \times K}$:

$$ \min_{C \in (S_K^{++})^N} J_{IVA-G}(\mathcal{W}, C) = \tilde{J}_{IVA-G}(\mathcal{W}) + \frac{KN}{2}. \quad (11) $$

Theorem 1 shows that $\tilde{\mathcal{W}}$ is a minimizer for $\tilde{J}_{IVA-G}$ if and only if $(\tilde{\mathcal{W}}, \hat{\mathcal{C}}(\tilde{\mathcal{W}}))$ is a minimizer for $J_{IVA-G}$, hence the equivalence. The advantage of relying on the proposed cost function that takes $\mathcal{C}$ as an input block variable is that it offers a more structured form, allowing an efficient use of a block alternating minimization scheme, and benefiting from sounded convergence, as we show in Section III.

### F. Ambiguities in IVA-G model

There are two ambiguities in the IVA-G generative model which correspond to information on the parameters that cannot be deduced from the observed data. The first is *permutation ambiguity*. New mixing matrices can be obtained by renumbering the SCVs, that is, by defining $\tilde{\mathcal{A}}^{[k]} = \mathcal{A}^{[k]} \mathcal{P}$ with $\mathcal{P}$ a permutation matrix. This change, however, does not modify the value of the cost function (8). The second is *scaling ambiguity*. For every $k \in \{1, \dots, K\}$ and $n \in \{1, \dots, N\}$, and for any $\alpha_n^{[k]} \in \mathbb{R} \setminus \{0\}$, we can replace $\mathbf{s}_n^{[k]}$ with $\hat{\mathbf{s}}_n^{[k]} = \alpha_n^{[k]} \mathbf{s}_n^{[k]}$ and $a_{m,n}^{[k]}$ with $\hat{a}_{m,n}^{[k]} = (\alpha_n^{[k]})^{-1} a_{m,n}^{[k]}$ for every $m \in \{1, \dots, N\}$. Those transformations leave the random vectors $\mathbf{x}_n^{[k]}$ unchanged, and consequently, do not affect the likelihood expression. Hence, the demixing matrices can only estimate the inverse of the ground truth mixing matrices, up to the permutation and scaling ambiguity. The latter ambiguity, moreover, raises the problem that a minimizing sequence of $J_{IVA-G}$ is not necessarily bounded. This motivates our proposition for a regularized version of the cost.

### G. Proposed regularized cost function

Let us note, for all $\mathcal{C} \in (S_K^{++})^N$ and $n \in \{1, \dots, N\}$, $(c_{n,k,k})_{1 \le k \le K} = \text{diag}(\mathcal{C}_n)$. Due to the scaling ambiguity raised above, for any $(\mathcal{W}, \mathcal{C})$, it is possible to rescale the coefficients to define $(\mathcal{W}', \mathcal{C}')$ such that $J_{IVA-G}(\mathcal{W}', \mathcal{C}') = J_{IVA-G}(\mathcal{W}, \mathcal{C})$ and that

$$ (\forall n \in \{1, \dots, N\})(\forall k \in \{1, \dots, K\}) \quad c'_{n,k,k} = 1. \quad (12) $$

To do this, one can set, for each $(k, n)$, $\alpha_n^{[k]} = \frac{1}{\sqrt{c_{n,k,k}}}$. As a consequence, if a minimizer of $J_{IVA-G}$ exists, then there exists at least another minimizer, satisfying $c_{n,k,k} = 1$ for all $(k, n)$. To mitigate this ambiguity, we thus propose to add a quadratic penalty to the cost function to control the distance to 1 of the diagonal coefficients of the precision matrices.<fnref n="1" />

**Remark 2** Let us notice that a minimizer of $J_{IVA-G}$, such that $c_{n,k,k} = 1$ has no reason to be qualitatively better than any other minimizer of $J_{IVA-G}$. The proposed regularization aims at reducing the number of distinct minima, and at ensuring some mathematical properties of the cost we will leverage to prove the convergence of our optimization algorithms. It is still possible, once convergence is reached, to rescale the sources a posteriori.

In addition to the scaling penalty term, we also propose to add an extra term to the cost function, to constrain the singular values of the recovered precision matrices to be positive by a minimum margin. This term aims at avoiding numerical issues that could arise at the boundary of the logarithm determinant definition domain, without damaging the quality of the results. In practice, the constraint is simply imposed, by adding an indicator function $\iota_{[\epsilon, +\infty)^K}$, equal to 0 for entries greater than $\epsilon$, $+\infty$ otherwise. Our objective is to impose that the components of the vector of singular values $\sigma_{C_n} \in \mathbb{R}^K$ of matrix $C_n$ are above a certain $\epsilon > 0$, for every $n$. We thus obtain the final form for our proposed (regularized) cost function, denoted by $J_{IVA-G}^{\text{Reg}}$:

$$ (\forall (\mathcal{W}, \mathcal{C}) \in (\mathcal{GL}_N(\mathbb{R})^K \times (S_K^{++})^N)) \quad J_{IVA-G}^{\text{Reg}}(\mathcal{W}, \mathcal{C}) = J_{IVA-G}(\mathcal{W}, \mathcal{C}) + \frac{\alpha}{2} \sum_{n=1}^{N} \| \text{diag}(C_n) - \mathbf{1}_K \|^2 \\ + \sum_{n=1}^{N} \iota_{[\epsilon, +\infty)^K} (\sigma_{C_n}), \quad (13) $$

$\alpha > 0$ is a hyper-parameter that controls the strength of the introduced regularization.

Our IVA-G method then aims to minimize (13). This is a challenging non-convex and non-smooth problem. The next section is dedicated to discuss further the properties of (13), and to design an efficient optimization algorithm to find a critical point of it.

## III. PROPOSED MINIMIZATION ALGORITHM

### A. Mathematical properties of $J_{IVA-G}^{\text{Reg}}$

In order to design an appropriate minimization algorithm for (13), let us examine the structure and properties of the cost function $J_{IVA-G}^{\text{Reg}}$. First, let us remark that minimizing (13) on $\mathcal{GL}_N(\mathbb{R})^K \times (S_K^{++})^N$ is equivalent to

$$ \minimize_{(\mathcal{W}, \mathcal{C}) \in \mathbb{R}^{N \times N \times K} \times \mathbb{R}^{K \times K \times N}} h(\mathcal{W}, \mathcal{C}) + f(\mathcal{W}) + g(\mathcal{C}), \quad (14) $$

<footnote n="1">¹In IVA-G-V and IVA-G-N, the scale ambiguity is managed by an ad-hoc renormalizing of the rows of the demixing matrices after each iteration of the minimization solver. Note that this normalization step also speeds up these algorithms, as it allows some simplification in the computations, as explained in [12, Section 4.B].</footnote>
with

$$2h(\mathcal{W}, \mathcal{C}) = \sum_{n=1}^{N} \mathrm{tr}(\mathcal{C}_n \mathbf{W}_n \hat{\mathbf{R}}_x \mathbf{W}_n^\top) + \alpha \|\mathrm{diag}(\mathcal{C}_n) - \mathbf{1}_K\|^2, \quad (15)$$

$$f(\mathcal{W}) = \begin{cases} -\sum_{k=1}^{K} \log |\det \mathbf{W}^{[k]}| & \text{if } (\forall k \in \{1, \dots, K\}) \mathbf{W}^{[k]} \in \mathcal{GL}_N(\mathbb{R}) \\ +\infty & \text{otherwise,} \end{cases} \quad (16)$$

$$g(\mathcal{C}) = \begin{cases} -\frac{1}{2} \sum_{n=1}^{N} \log \det \mathcal{C}_n & \text{if } (\forall n \in \{1, \dots, N\}) \mathcal{C}_n - \epsilon \mathbf{I}_K \in \mathcal{S}_K^+ \\ +\infty & \text{otherwise.} \end{cases} \quad (17)$$

As we already mentioned, $\epsilon > 0$ serves to obtain better regularity for the cost function (with closed definition domain). It is typically taken very small (e.g., $\epsilon = 10^{-12}$ in our experiments) to ensure the problem is essentially equivalent to a search for a maximum of the likelihood. A similar strategy was adopted in [33]. We can deduce from (4) that the presence of $\epsilon$ only plays a role in our iterative algorithm if, for any $n \in \{1, \dots, N\}$, the estimated $\mathcal{C}_n$ has a singular value in $(0, \epsilon)$ at some iteration. In this case, this singular value is set to $\epsilon$ instead. As we show in Appendix F, this case did not happen in practice in our experimental settings as soon as $\epsilon$ is less than or equal to $10^{-1}$. Formulation (14) has the advantage of isolating a differentiable term $h$ acting on both sets of variables $(\mathcal{W}, \mathcal{C})$, and two non-differentiable terms, $f$ and $g$, acting separately on $\mathcal{W}$ or $\mathcal{C}$.

a) Function $h$: Let us first study function $h$ acting jointly on $\mathcal{W}$ and $\mathcal{C}$ variables. We state the following lemmas, regarding the expression and smoothness properties, of its partial gradients, with respect to the first or second entry, the other being fixed.

**Lemma 1** The partial gradient of $h$ with respect to variable $\mathcal{W}$, evaluated at $(\mathcal{W}, \mathcal{C}) \in \mathbb{R}^{N \times N \times K} \times \mathbb{R}^{K \times K \times N}$, reads:

$$\nabla_{\mathcal{W}} h(\mathcal{W}, \mathcal{C}) = (\nabla_{\mathcal{W}^{[k]}} h(\mathcal{W}, \mathcal{C}))_{1 \le k \le K} \in \mathbb{R}^{N \times N \times K}, \quad (18)$$

where, for all $k \in \{1, \dots, K\}$,

$$\nabla_{\mathcal{W}^{[k]}} h(\mathcal{W}, \mathcal{C}) = \left( \frac{\partial h(\mathcal{W}, \mathcal{C})}{\partial w_{n,m}^{[k]}} \right)_{1 \le n, m \le N}$$

$$= \left( \begin{array}{cccc} [C_1 W_1 \hat{R}_x]_{k,(k-1)N+1} & \cdots & [C_1 W_1 \hat{R}_x]_{k,(k-1)N+N} \\ \vdots & \ddots & \vdots \\ [C_N W_N \hat{R}_x]_{k,(k-1)N+1} & \cdots & [C_N W_N \hat{R}_x]_{k,(k-1)N+N} \end{array} \right)$$

Moreover, for every $\mathcal{C} \in \mathbb{R}^{K \times K \times N}$, $\nabla_{\mathcal{W}} h(\mathcal{W}, \mathcal{C})$ is Lipschitz continuous, with modulus

$$L_{\mathcal{W}}(\mathcal{C}) = \rho_{\mathcal{C}} \varrho_{\hat{\mathbf{R}}_{\mathbf{x}}}, \quad (19)$$

where $\hat{\mathbf{R}}_{\mathbf{x}}^{[k]}$ is the matrix obtained by extraction of the columns $kN + 1$ to $(k+1)N$ of $\hat{\mathbf{R}}_{\mathbf{x}}$,

$$\varrho_{\hat{\mathbf{R}}_{\mathbf{x}}} = \max_{1 \le k \le K} \| \hat{\mathbf{R}}_{\mathbf{x}}^{[k]} \|_{\mathcal{S}} \quad (20)$$

and

$$(\forall \mathcal{C} \in \mathbb{R}^{K \times K \times N}) \quad \rho_{\mathcal{C}} = \max_{1 \le n \le N} \| \mathcal{C}_n \|_{\mathcal{S}}. \quad (21)$$

Proof: See Appendix A.

**Lemma 2** The partial gradient of $h$ with respect to $\mathcal{C}$, evaluated at $(\mathcal{W}, \mathcal{C}) \in \mathbb{R}^{N \times \hat{N} \times K} \times \mathbb{R}^{K \times K \times N}$, leads to

$$\nabla_{\mathcal{C}} h(\mathcal{W}, \mathcal{C}) = (\nabla_{\mathcal{C}_n} h(\mathcal{W}, \mathcal{C}))_{1 \le n \le N} \in \mathbb{R}^{K \times K \times N}, \quad (22)$$

where, for all $n \in \{1, \dots, N\}$,

$$\nabla_{\mathcal{C}_n} h(\mathcal{W}, \mathcal{C}) = \frac{1}{2} \mathbf{W}_n \hat{\mathbf{R}}_x \mathbf{W}_n^\top + \alpha (\mathrm{Diag}(\mathcal{C}_n) - \mathbf{I}_K).$$

Moreover, for every $\mathcal{W} \in \mathbb{R}^{N \times N \times K}$, $\nabla_{\mathcal{W}} h(\mathcal{W}, \cdot)$ is Lipschitz continuous with constant modulus

$$L_{\mathcal{C}} = \alpha. \quad (23)$$

Proof: See Appendix B.

According to Lemmas 1 and 2, both partial derivatives of $h$ are well defined and continuous with respect to $\mathcal{W}$ and $\mathcal{C}$, which shows in particular that $h$ is a $C^1$ function.

b) Functions $f$ and $g$: Functions $f$ and $g$ are proper (i.e., finite-valued at least at one point), and lower-semicontinuous. Furthermore, function $g$ is convex (see, for e.g., [34, Example 24.66]), for a proof), while function $f$ is not. Both functions $f$ and $g$ are not differentiable but still, it is possible to manipulate them efficiently, for minimization purpose, through their proximity operator [35].<fnref n="2" /> The following lemmas provide the expression for these operators, that will then be perused in our algorithm. The proofs for these lemmas mainly rely on the fact that both $f$ and $g$ are so-called spectral functions [36], depending only on the spectral values of their inputs. Note that, despite the non-convexity of $f$, its proximity operator is still uniquely defined, as shown in the proof for Lemma 3.

**Lemma 3** Let $\mathcal{W}' \in \mathbb{R}^{N \times N \times K}$, and some $c > 0$. We define the proximity operator of $f$ at $\mathcal{W}'$ as

$$\begin{aligned} \mathrm{prox}_{c,f}(\mathcal{W}') &= \underset{\mathcal{W} \in \mathbb{R}^{N \times N \times K}}{\operatorname{argmin}} \frac{1}{2} \|\mathcal{W} - \mathcal{W}'\|^2 + cf(\mathcal{W}) \\ &= (U_{\mathcal{W}'[k]} \mathrm{Diag}(\sigma_{\mathcal{W}'[k]}) V_{\mathcal{W}'[k]})_{1 \le k \le K}, \quad (24) \end{aligned}$$

where $(\forall k \in \{1, \dots, K\}) \mathcal{W}'[k] = U_{\mathcal{W}'[k]} \mathrm{Diag}(\sigma_{\mathcal{W}'[k]}) V_{\mathcal{W}'[k]}$ is the singular value decomposition of $\mathcal{W}'[k]$, and

$$\sigma_{\mathcal{W}^{[k]}} = \frac{\sigma_{\mathcal{W}^{[k]}} + \sqrt{\sigma_{\mathcal{W}^{[k]}}^2 + 4c}}{2}. \quad (25)$$

Proof: See Appendix C.

**Lemma 4** Let $\mathcal{C}' \in \mathbb{R}^{K \times K \times N}$, and some $c > 0$. The proximity operator of $g$ at $\mathcal{C}'$ is given by

$$\begin{aligned} \mathrm{prox}_{c,g}(\mathcal{C}') &= \underset{\mathcal{W} \in \mathbb{R}^{K \times K \times N}}{\operatorname{argmin}} \frac{1}{2} \|\mathcal{C}' - \mathcal{C}'\|^2 + cg(\mathcal{C}) \\ &= (U_{\mathcal{C}'_n} \mathrm{Diag}(\sigma_{\mathcal{C}'_n}) V_{\mathcal{C}'_n})_{1 \le n \le N} \quad (26) \end{aligned}$$

where $(\forall n \in \{1, \dots, N\})$, $\mathcal{C}'_n = U_{\mathcal{C}'_n} \mathrm{Diag}(\sigma_{\mathcal{C}'_n}) V_{\mathcal{C}'_n}$ is the singular value decomposition of $\mathcal{C}'_n$, and

$$\sigma_{\mathcal{C}_n} = \max \left( \epsilon, \frac{\sigma_{\mathcal{C}'_n} + \sqrt{(\sigma_{\mathcal{C}'_n})^2 + 2c}}{2} \right). \quad (27)$$

Proof: See Appendix D.

<footnote n="2">²See also https://proximity-operator.net/</footnote>
## B. Proposed PALM-IVA-G algorithm

As shown in the previous subsection, the minimization of (13) amounts to solving Problem (14), which has a structure particularly well suited to block alternating minimization. More precisely, we have shown that the cost function includes partially Lipschitz differentiable terms acting on both variables ($W, C$) (Lemmas 1 and 2), as well as two terms acting separately on these variables. Despite being non-differentiable, these terms are proper, lower-semicontinuous and with a tractable proximity operator (Lemmas 3 and 4). These results pave the way for applying a block alternating proximal gradient algorithm, as studied for instance in [37], [38], [39]. Here, we opted for PALM introduced in [29], due to its powerful convergence results. We adapted PALM mechanism to the minimization of the cost (13) and thus designed PALM-IVA-G presented in Alg. 1. In PALM initial study, the convergence was shown in the case of sequential block updates. Here, we instead opted for a more versatile update scheme that follows the so-called essentially cyclic rule [37], allowing each block to be updated more than once, per main iteration. This assumption hence gives more flexibility to our algorithm and it is straightforward to adapt the proof given in [29] to this case, using a technique similar to [37].

At each step $i \in \mathbb{N}$ of the PALM-IVA-G main loop, we update $\mathcal{W}$ (resp. $\mathcal{C}$) a number $n_{\mathcal{W}}(i) \le \bar{n}_{\mathcal{W}}$ (resp. $n_{\mathcal{C}}(i) \le \bar{n}_{\mathcal{C}}$) of times, with $\bar{n}_{\mathcal{W}}$ and $\bar{n}_{\mathcal{C}}$ positive integers. The updates include gradient and proximal steps, as follows. First, a gradient step on $h$ is performed with respect to the active block, $\mathcal{W}$ or $\mathcal{C}$, with a respective positive step size $c_{\mathcal{W}}^{(i)}$ or $c_{\mathcal{C}}$. Then, a proximal step on $f$ or $g$ is run, using the same step size. Inner and outer loops are controlled by a maximum number of iterations and furthermore include early stopping tests, based on comparing the following quantities with the (small) precision parameter $\delta > 0$:

$$ (\forall (\mathcal{W}, \mathcal{W}') \in (\mathbb{R}^{N \times N \times K})^2) \\ \theta_{\mathcal{W}}(\mathcal{W}, \mathcal{W}') = \max_{\substack{1 \le n \le N \\ 1 \le k \le K}} \frac{\|\mathbf{w}_n^{[k]} - \mathbf{w}_n^{[k]}\|^2}{2N}, \quad (28) $$

$$ (\forall (\mathcal{C}, \mathcal{C}') \in (\mathbb{R}^{K \times K \times N})^2) \\ \theta_{\mathcal{C}}(\mathcal{C}, \mathcal{C}') = \max_{\substack{1 \le n \le N \\ 1 \le k \le K}} \frac{\|\mathbf{c}_{n,k}' - \mathbf{c}_{n,k}\|^2}{2K}, \quad (29) $$

where, for $n \in \{1, \dots, N\}$ and $k \in \{1, \dots, K\}$, the notation $c_{n,k}$ refers to the $k$-th row of the matrix $\mathcal{C}_n$.

Let us now turn to Section IV, with the aim of establishing the convergence of PALM-IVA-G iterates $(W^{(i)}, C^{(i)})_{i \in \mathbb{N}}$. Practical settings for PALM-IVA-G hyper-parameters will be discussed in Section V.

## IV. CONVERGENCE RESULT

Let us state our convergence theorem for the proposed PALM-IVA-G algorithm.

**Theorem 2 (Convergence of PALM-IVA-G)** Assuming the setting of Algorithm 1, the infinite sequence $(W^{(i)}, C^{(i)})_{i \in \mathbb{N}}$ converges to a critical point $(W^*, C^*)$ of function $J_{\text{IVA-G}}^{\text{Reg}}$ given in (13).

### Algorithm 1 PALM-IVA-G

**Require:** Empirical covariance $\hat{R}_x$, initial tensors $(W^{(0)}, C^{(0)}) \in \mathcal{G}\mathcal{L}_N(\mathbb{R})^K \times (\epsilon I_K + S_K^+)^N$, penalty weight $\alpha > 0$, step sizes $\gamma_C \in (0, 2)$, $\gamma_W \in (0, 1)$, maximal inner/outer loops $\bar{n}_W \in \mathbb{N} \setminus \{0\}$, $\bar{n}_C \in \mathbb{N} \setminus \{0\}$, $\bar{N} \in \mathbb{N} \setminus \{0\}$, precision $\delta > 0$.

**Ensure:** $(W_{\text{out}}, C_{\text{out}})$
1: $\triangleright \text{ Initialization}$ ▷
2: $\mathcal{W}^{(0,0)} \leftarrow \mathcal{W}^{(0)}$
3: $C^{(0,0)} \leftarrow C^{(0)}$
4: $c_C \leftarrow \frac{\gamma_C}{\alpha}$
5: $i \leftarrow 0$
6: $\theta_{\text{ext}}^{(0)} \leftarrow +\infty$
7: $\triangleright \text{ Start Main Loop}$ ▷
8: **while** $\theta_{\text{ext}}^{(i)} > \delta$ **or** $i < \bar{N}$ **do**

9: $\triangleright \text{ Update } \mathcal{W}$ ▷
10: $c_{\mathcal{W}}^{(i)} \leftarrow \frac{\gamma_{\mathcal{W}}}{L_{\mathcal{W}}(C^{(i)})}$ using (19)
11: $j \leftarrow 0$
12: $\theta_{\text{int}}^{(0)} \leftarrow +\infty$
13: **while** $\theta_{\text{int}}^{(j)} > \delta$ **or** $j < \bar{n}_{\mathcal{W}}$ **do**
14: $\mathcal{W}^{(i,j+1)} \leftarrow \text{prox}_{c_{\mathcal{W}}^{(i)}}(\mathcal{W}^{(i,j)} - c_{\mathcal{W}}^{(i)} \nabla_{\mathcal{W}} h(\mathcal{W}^{(i,j)}, C^{(i)}))$
15: using (18) and (24)
16: $\theta_{\text{int}}^{(j)} \leftarrow \theta_{\mathcal{W}}(\mathcal{W}^{(i,j+1)}, \mathcal{W}^{(i,j)})$
17: $j \leftarrow j + 1$
18: $\mathcal{W}^{(i+1,0)} \leftarrow \mathcal{W}^{(i,j)}$
19: $\mathcal{W}^{(i+1)} \leftarrow \mathcal{W}^{(i+1,0)}$
20: $\triangleright \text{ Update } \mathcal{C}$ ▷
21: $j \leftarrow 0$
22: $\theta_{\text{int}}^{(0)} \leftarrow +\infty$
23: **while** $\theta_{\text{int}}^{(j)} > \delta$ **or** $j < \bar{n}_{\mathcal{C}}$ **do**
24: $C^{(i,j+1)} \leftarrow \text{prox}_{c_{\mathcal{C}}g}(C^{(i,j)} - c_{\mathcal{C}} \nabla_{\mathcal{C}} h(C^{(i,j)}, \mathcal{W}^{(i+1)}))$
25: using (22) and (26)
26: $\theta_{\text{int}}^{(j)} \leftarrow \theta_{\mathcal{C}}(C^{(i,j+1)}, C^{(i,j)})$
27: $j \leftarrow j + 1$
28: $\mathcal{C}^{(i+1,0)} \leftarrow \mathcal{C}^{(i,j)}$
29: $\mathcal{C}^{(i+1)} \leftarrow \mathcal{C}^{(i+1,0)}$
30: $\triangleright \text{ Evaluate Stopping Criteria}$ ▷
31: $\theta_{\mathcal{W}}^{(i+1)} \leftarrow \theta_{\mathcal{W}}(\mathcal{W}^{(i+1)}, \mathcal{W}^{(i)})$ using (28)
32: $\theta_{\mathcal{C}}^{(i+1)} \leftarrow \theta_{\mathcal{C}}(\mathcal{C}^{(i+1)}, \mathcal{C}^{(i)})$ using (29)
33: $\theta_{\text{ext}}^{(i+1)} \leftarrow \max(\theta_{\mathcal{W}}^{(i+1)}, \theta_{\mathcal{C}}^{(i+1)})$
34: $i \leftarrow i + 1$
35: $(W_{\text{out}}, C_{\text{out}}) \leftarrow (\mathcal{W}^{(i)}, C^{(i)})$
return $(W_{\text{out}}, C_{\text{out}})$

Here, a critical point is defined as in [29, Rem.1 (iv)], i.e., $0 \in \partial J_{\text{IVA-G}}^{\text{Reg}}(W^*, C^*)$, where $\partial$ denotes the limiting subdifferential operator. The proof of the above result relies on [29, Theorem 1] in the case of a cyclic update of the blocks. The latter states the convergence of a generic form of **PALM** method under several assumptions regarding the properties of the cost function, and provided the sequence is bounded. It is hence sufficient to show that such conditions hold in our case, to prove Theorem 2. In the previous section, we have already seen that $f$ and $g$ are proper and lower semi-continuous functions such that the proximal operators $\text{prox}_{c_f}$ and $\text{prox}_{c_g}$
are defined for all $c > 0$, at any point $W \in GL_N(\mathbb{R})^K$ and $C \in (S_K^{++})^N$, respectively. We also outlined that $h$ is a $C^1$ function and for every given $W' \in \mathbb{R}^{N \times N \times K}$ (resp. $C' \in \mathbb{R}^{K \times K \times N}$), the function $C \mapsto h(W', C)$ (resp. $W \mapsto h(W, C')$ is $C_{L_C(W')}^{1,1}$, i.e. its partial gradient $C \mapsto \nabla_C h(W', C)$ is globally Lipschitz with modulus $L_C(W')$ (resp. $C_{L_W(C')}^{1,1}$). We note that those first properties were necessary to well-defining the algorithm itself.

Using our notation, to complete the proof, we need to demonstrate each item of the following Assumption 1 is satisfied.

**Assumption 1**

1) $\inf_{\mathbb{R}^{N \times N \times K} \times \mathbb{R}^{K \times K \times N}} J_{\text{IVA-G}}^{\text{Reg}} > -\infty$.

2) $\nabla h$ is Lipschitz continuous on bounded subsets of $\mathbb{R}^{N \times N \times K'} \times \mathbb{R}^{K \times K \times N}$.

3) $(W^{(i)}, C^{(i)})_{i \in \mathbb{N}}$ is bounded.

4) There exists $(\lambda_W^+, \lambda_W^-; \lambda_C^+, \lambda_C^-) > 0$ such that:

$$ (\forall i \in \mathbb{N}) \quad \lambda_W^- \le L_W(C^{(i)}) \le \lambda_W^+,$$ \\ 
and \\ 
$(\forall i \in \mathbb{N}) \quad \lambda_C^- \le L_C(W^{(i)}) \le \lambda_C^+.$

5) $J_{\text{IVA-G}}^{\text{Reg}}$ is a KL function.

Indeed, following the steps of [29], the update scheme of the algorithms makes the sequence of the costs decreasing, and using the first item of the above assumption, this sequence has a finite limit. Then, items 2) to 4) enable to prove that the limit set of $(W^{(i)}, C^{(i)})_{i \in \mathbb{N}}$ is nonempty, compact, and contains only critical points of $J_{\text{IVA-G}}^{\text{Reg}}$. Finally, the last item is used to prove that $(W^{(i)}, C^{(i)})_{i \in \mathbb{N}}$ is a Cauchy sequence, hence convergent. The verification of Assumption 1 is technically involved and is therefore provided in Appendix E.

## V. EXPERIMENTAL RESULTS

We now present a set of experiments to assess the quantitative, qualitative, and computational performance of PALM-IVA-G on tasks of independent Gaussian source separation <fnref n="3" />. A. Experimental protocol

1) **Qualitative evaluation:** In all experiments, the proposed method PALM-IVA-G, as well as the benchmarks, are evaluated by means of the so-called jISI (joint Intersymbol Interference) score, also used in [12]. This score is an extension of the ISI score that was introduced in [40] to assess the performance of ICA methods. The jISI score measures the correspondence between demixing matrices $W = [W^{[1]}, \dots, W^{[K]}]$ and the ground truth mixing matrices $A = [A^{[1]}, \dots, A^{[K]}]$, up to a common permutation and a rescaling of their rows. Let us note

$$ (\forall k \in \{1, \dots, K\}) \quad G^{[k]} = W^{[k]} A^{[k]} \quad (30) $$

and, for every $(m, n) \in \{1, \dots, N\}^2$, the mean $\bar{g}_{n,m} = \sum_{k=1}^K |g_{n,m}^{[k]}|$ of the $(n, m)$-th entry of the tensor $G = [\mathbf{G}^{[1]}, \dots, \mathbf{G}^{[K]}]$. The jISI score of the pair $(W, A)$ is

$$ \text{jISI}(W, A) = \frac{1}{2N(N-1)} \left[ \sum_{n=1}^{N} \left( \sum_{m=1}^{N} \frac{\bar{g}_{n,m}}{\max_p \bar{g}_{n,p}} - 1 \right) + \sum_{m=1}^{N} \left( \sum_{n=1}^{N} \frac{\bar{g}_{n,m}}{\max_p \bar{g}_{p,m}} - 1 \right) \right]. \quad (31) $$

As defined in (31), the jISI score is a real number between 0 and 1, which equals 0 if and only if $\overline{G} = (\bar{g}_{n,m})_{1 \le n, m \le N}$ has exactly one positive coefficient by row and by column, that is, when it is a (possibly permuted) diagonal matrix. This happens when the $(W^{[k]} A^{[k]})_{1 \le k \le K}$ are simultaneously permuted diagonal matrices, which is the best situation one can expect from a source separation step. Hence, the smallest the jISI score, the better quality for the source separation.

2) **Benchmark methods:** Our algorithm PALM-IVA-G is compared against two state of the art algorithms for independent Gaussian source separation, namely IVA-G-V and IVA-G-N, both introduced in [12]. Those two algorithms perform the minimization of the cost function $\tilde{J}_{\text{IVA-G}}$ (9) (which, as we remind, only depends on the variable $W$) using, respectively, a gradient descent and a Newton's method. Both implement a backtracking linesearch and a normalization of the rows of the demixing matrix at each iteration. Hyper-parameters are set as recommended in [12]. Let us emphasize that both of these algorithms are empirical and, to our knowledge, do not benefit from established convergence guarantees, though in practice we did not observe any failure.

For every experiment and each algorithm, we obtain an estimate for which we record the jISI score reached at convergence (i.e., when the stopping criterion is reached) and the computational time, in seconds, denoted $T$, to reach this convergence point. For fair comparison, all algorithms are initialized with the same random $W^{(0)}$ drawn from the standard normal distribution. All algorithms are implemented in Python 3.10.12 and run on a Dell Precision 5820 Workstation with 11th Gen Intel(R) Core(TM) i9-10900X at 3.00GHz, equipped with 32Go Ram.

3) **Synthetic dataset generation:** For each experimental setup, we define a generative model defined by the ground-truth variables $(A, \Sigma)$, and use this model to generate source data $S$ of length $V = 10000$, and then mix it into observed data $\mathcal{X}$. The goal is to recover estimates of $([A^{[1]}]^{-1}, \dots, (A^{[K]})^{-1}]$, $[\Sigma_1^{-1}, \dots, \Sigma_N^{-1}]$) up to the scaling and permutation ambiguities. All trials are repeated over 100 independent runs, and we compute the mean $\mu_{\text{jISI}}$ (resp. $\mu_{\text{T}}$) and standard deviation $\sigma_{\text{jISI}}$ (resp. $\sigma_{\text{T}}$) values for the jISI score (resp. computational time $T$). Most experiments are conducted for various values of the dimensions $(K, N)$, specified in each test case. Depending on the nature of the phenomenon modeled, the covariance and the mixing matrices may have various properties, leading us to define several sets of experiments, detailed hereafter.

We aim at exploring the impact of the overall level of correlation across the datasets (given by the extra-diagonal coefficients of the SCVs covariance matrices) and the variability of those correlations. Our generative model used to simulate SCVs depends on two parameters $\rho = [\rho_1, \dots, \rho_N]^\top \in [0, 1]^N$ and $\lambda \in [0, 1]$. Given those, we compute a ground truth tensor $\Sigma = [\Sigma_1, \dots, \Sigma_N]$ with, for every $n$,

$$ \Sigma_n = \rho_n \mathbf{1}_n^\top + \frac{\lambda}{R} \mathbf{Q} \mathbf{Q}^\top + \eta_n \mathbf{I}_K \quad (32) $$

with

$$ \eta_n = 1 - \rho_n - \lambda \in [0, 1], \quad (33) $$

matrix $Q$ randomly drawn in $\mathbb{R}^{K \times R}$ with elements $q_{i,j} \sim N(0, 1)$ mutually independent, and $R \in \mathbb{N} \setminus \{0\}$ a predefined rank value.

<footnote n="3">³The Python code for PALM-IVA-G and all the software used to run the experiments can be found at https://github.com/Mentell/independent_vector_analysis</footnote>
Then, we have, for every $n$, and every entry $(i, j)$,

$$ \left\{ \begin{array}{l} \mathbb{E}\{(\Sigma_n)_{i,j}\} = \rho_n + (\lambda + \eta_n)\delta_{i,j} \\ \text{Var}\{(\Sigma_n)_{i,j}\} = \frac{\lambda^2}{R}(1 + \delta_{i,j}). \end{array} \right. \quad (34) $$

This means that $\rho$ controls the average correlation across the datasets while $\lambda$ controls the variability between the correlations. The third term ensures that the covariance matrices we use are positive definite. In our experiments, we opt for four scenarios, corresponding to low/high variability, and low/high correlation, as defined below:

* Case A: low correlation, low variability. $\lambda = 0.04$ and $\rho$, regularly sampled in $[0.2, 0.3]^N$,

* Case B: low correlation, high variability. $\lambda = 0.25$ and $\rho$, regularly sampled in $[0.2, 0.3]^N$,

* Case C: high correlation, low variability. $\lambda = 0.04$ and $\rho$, regularly sampled in $[0.6, 0.7]^N$,

* Case D: high correlation, high variability. $\lambda = 0.25$ and $\rho$, regularly sampled in $[0.6, 0.7]^N$,

The sources are expected to be easier to separate in case D, while case A is more challenging, since increasing the correlation or the variability of the sources decreases the Cramer-Rao Lower-Bound of the ML estimator and thus tends to improve separation [18]. We investigate each of the four cases for various sizes $K \in \{5, 10, 20\}$ and $N \in \{10, 20\}$ and set $R = K + 10$. In all experiments, the ground truth tensor $\mathcal{A}$ is simulated by drawing its entries independently from a zero-mean Gaussian distribution with a standard deviation of 1.

For each simulated pair $(\mathcal{A}, \Sigma)$, we build the corresponding input data $\mathcal{X}$, following the mixing model (1). Note that, as recommended in [12], we systematically performed a whitening of the data, before applying the algorithms. This step amounts to multiplying, for every $k \in \{1, \dots, K\}$, the latent sources by full-rank matrices $(\mathbf{B}^{[k]})_{1 \le k \le K}$, and solving the demixing problem on $\mathbf{X}^{[k]} = \mathbf{B}^{[k]}\mathbf{X}^{[k]} = \mathbf{B}^{[k]}\mathbf{A}^{[k]}\mathbf{S}^{[k]}$. The whitening matrices are set to decrease the spectral norm of $\hat{\mathbf{R}}_{\mathcal{X}}$, with the aim of improving the conditioning of the loss function, and hence accelerate the empirical convergence of the methods.

4) Algorithmic settings: The implementation of the proposed PALM-IVA-G algorithm requires the setting of (i) the step sizes, (ii) the regularization weight, and (iii) the stopping conditions for internal and external loops.

We set the step sizes to constant values $(\gamma_C, \gamma_W) = (1.99, 0.99)$, i.e., as large as possible to meet the assumptions of the convergence theorem. The penalty parameter $\alpha$ appears to have little influence on the performance of PALM-IVA-G, in terms of the jISI metric and computational time, as long as it is chosen in a reasonable range. This can be observed empirically in Fig. 1 summarizing the results for Case A and various values of $\alpha$. Indeed, except for $\alpha = 10$ which seems to yield systematically slower computations, we cannot observe that a value of $\alpha$ consistently gives better results than the others. Similar behaviors were obtained for Cases B to D. For the experiments, we thus retain $\alpha = 1$, as it achieves a good compromise in terms of time complexity.

The PALM-IVA-G algorithm and its competitors IVA-G-V and IVA-G-N, are run until a certain stopping criterion is reached, with a maximum number of $\bar{N} = 20000$ iterations. The precision threshold is set to $\delta = 10^{-10}$, and the maximum numbers of iterations within the internal loops of PALM-IVA-G are $\bar{n}_C = 1$ and $\bar{n}_W = 15$. For IVA-G-V (resp. IVA-G-N), we monitor only the value of (28) between consecutive iterates, stopping once lower than $\delta = 10^{-6}$ (resp. $\delta = 10^{-7}$).

<figure><img src="image_2.png" /><figcaption>Fig. 1: Distributions of jISI scores and computation times using PALM-IVA-G for Case A over 20 runs, with $\alpha$ taking values in $\{0.1, 0.5, 1, 5, 10\}$. The jISI score varies little with $\alpha$, the latter impacting mostly the time, with good compromise at $\alpha = 1$.</figcaption></figure>

Note that the values of $(\bar{n}_C, \bar{n}_W, \delta)$ have been empirically set to reach the best trade-off between the jISI score and the computational time, ensuring fair comparisons of the methods.

The validity of our settings can be visually assessed in Fig. 2, showing the evolution of the cost function using our PALM-IVA-G algorithm in a representative example.

## B. Results

We now present the results of the experiments and comment on the strengths and weaknesses of our algorithm in comparison with the benchmarks. The results are illustrated in Figs. 3 and 4, and summarized in Tab. I. In Fig. 3 and 4, we display, for each of the three algorithms, two boxplots representing respectively the distribution of the jISI scores (in a log scale) and the distribution of the computation times in seconds. In Tab. I, we highlight in bold font (resp. italic bold font) the best results in terms of $\mu_{j\text{ISI}}$ (resp. $\mu_T$), considering as having a similar quality the jISI scores with a difference less than $10^{-4}$ (resp. computation times with a difference less than $10^{-2}$).

All algorithms achieve what appears to be optimal separation in easy cases (B and D). In difficult cases, it is generally either IVA-G-V or PALM-IVA-G that have the best jISI, with an advantage for PALM-IVA-G in small dimensions and an advantage for IVA-G-V in larger dimensions, but in all cases the performance remains close. IVA-G-N is generally not as
<figure><img src="image_3.png" /><figcaption>Fig. 2: Top: Empirical convergence of PALM-IVA-G, on a synthetic example from Case B, with $(K, N) = (20, 20)$, $\lambda = 0.25$ and $\rho \in [0.2, 0.3]^T$. Cost function across iterations (left) and time in second (right). Bottom: On the same example, evolution of jISI score along iterations (left) and time in seconds (right), for the compared methods.</figcaption></figure>

good, but manages to keep its performance close to that of the other algorithms, except for Case A, where it gives jISI values 20 to 50% higher than the other algorithms on average.

Computations of PALM-IVA-G are tractable and stay under two minutes. The running time is sensitive to the number of sources and seems to grow linearly with the number of datasets within the tested dimensions. For $N = 10$, our algorithm takes less than 15 seconds to run on average in all cases. Meanwhile, IVA-G-V is less sensitive to $N$ and also manages to separate the sources in two minutes on average at most for the dimensions in the experiment. However, IVA-G-N is the slowest of all three algorithms, it takes several tens of seconds in small dimensions for Cases A and C, and its computational cost becomes prohibitive in larger dimensions, taking an average time of 15 minutes to run in Case A for $K = 20$ and $N = 20$.

In contrast to gradient descent, iterations using Newton's method or proximal gradient are more informed and are expected to find the local minimum more efficiently. On the other hand, these methods are more costly, since they involve Hessian inversions for IVA-G-N, and SVDs for our PALM-IVA-G, and such cost is not necessarily offset by a gain in number of iterations. Besides, the update scheme of PALM-IVA-G implies many more updates of the block $W$ than the block $C$, so most of the matrices for which we compute the SVD are of size $N \times N$ rather than $K \times K$, this is why the computation time of PALM-IVA-G increases faster with respect to $N$ than $K$. Finally, inside an iteration of PALM-IVA-G, the nature of the updates of $W$ and $C$ allows to parallelize their computations for $k \in \{1, \dots, K\}$ and for $n \in \{1, \dots, N\}$ respectively, which enables an efficient implementation.

In conclusion, in addition to its established convergence guarantees, PALM-IVA-G appears to be competitive with the state of the art IVA-G algorithm, consistently achieving good jISI scores and taking reasonable time to run, especially when the number of sources is not too high.

<figure><img src="image_4.png" /><figcaption>Case A: $\rho \in [0.2, 0.3]^N$ and $\lambda = 0.04$</figcaption></figure>

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /><figcaption>Case B: $\rho \in [0.2, 0.3]^N$ and $\lambda = 0.25$</figcaption></figure>

<caption>Fig. 3: Distributions of jISI score and computational time in seconds, for PALM-IVA-G (green), IVA-G-V (red) and IVA-G-N (blue), for Case A and Case B.</caption>
<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 4: Distributions of jISI score and computational time in seconds, for PALM-IVA-G (green), IVA-G-V (red) and IVA-G-N (blue), for Cases C and D.</figcaption></figure>

<table><caption>TABLE I: Averaged jISI scores, $\mu_{jISI}$, and averaged computational times, $\mu_T$, in seconds, for PALM-IVA-G, IVA-G-V and IVA-G-N (from top to bottom). Best (i.e., lowest) jISI results (resp. lowest times) are highlighted in bold (resp. italic bold).</caption><thead><tr><th rowspan="2"></th><th rowspan="2"></th><th colspan="2">K = 5</th><th colspan="2">K = 10</th><th colspan="2">K = 20</th></tr><tr><th>N = 10</th><th>N = 20</th><th>N = 10</th><th>N = 20</th><th>N = 10</th><th>N = 20</th></tr></thead><tbody><tr><td rowspan="4">PALM-IVA-G</td><td rowspan="2">Case A</td><td>μ<sub>jISI</sub></td><td>9.79E-02</td><td>1.36E-01</td><td>5.40E-02</td><td>6.74E-02</td><td>3.88E-02</td><td>4.18E-02</td></tr><tr><td>μ<sub>T</sub></td><td>3.0</td><td>23.8</td><td>4.3</td><td>41.2</td><td>9.2</td><td>84.2</td></tr><tr><td rowspan="2">Case B</td><td>μ<sub>jISI</sub></td><td>2.14E-02</td><td>2.20E-02</td><td>1.40E-02</td><td>1.41E-02</td><td>1.13E-02</td><td>1.13E-02</td></tr><tr><td>μ<sub>T</sub></td><td>0.5</td><td>3.1</td><td>0.6</td><td>3.7</td><td>1.5</td><td>8.6</td></tr><tr><td rowspan="2">Case C</td><td>μ<sub>jISI</sub></td><td>4.63E-02</td><td>5.29E-02</td><td>2.51E-02</td><td>2.63E-02</td><td>1.70E-02</td><td>1.70E-02</td></tr><tr><td>μ<sub>T</sub></td><td>4.9</td><td>39.2</td><td>6.6</td><td>49.7</td><td>13.6</td><td>107.2</td></tr><tr><td rowspan="2">Case D</td><td>μ<sub>jISI</sub></td><td>9.45E-03</td><td>9.46E-03</td><td>6.03E-03</td><td>6.03E-03</td><td>4.97E-03</td><td>4.98E-03</td></tr><tr><td>μ<sub>T</sub></td><td>0.7</td><td>4.6</td><td>0.9</td><td>5.6</td><td>2.4</td><td>14.9</td></tr><tr><td rowspan="4">IVA-G-V</td><td rowspan="2">Case A</td><td>μ<sub>jISI</sub></td><td>1.01E-01</td><td>1.37E-01</td><td>5.61E-02</td><td>6.64E-02</td><td>3.81E-02</td><td>4.11E-02</td></tr><tr><td>μ<sub>T</sub></td><td>1.8</td><td>5.4</td><td>12.2</td><td>30.4</td><td>44.6</td><td>123.7</td></tr><tr><td rowspan="2">Case B</td><td>μ<sub>jISI</sub></td><td>2.15E-02</td><td>2.20E-02</td><td>1.40E-02</td><td>1.41E-02</td><td>1.13E-02</td><td>1.13E-02</td></tr><tr><td>μ<sub>T</sub></td><td>1.1</td><td>2.3</td><td>6.0</td><td>15.2</td><td>28.5</td><td>58.2</td></tr><tr><td rowspan="2">Case C</td><td>μ<sub>jISI</sub></td><td>4.68E-02</td><td>5.03E-02</td><td>2.47E-02</td><td>2.58E-02</td><td>1.63E-02</td><td>1.67E-02</td></tr><tr><td>μ<sub>T</sub></td><td>1.6</td><td>4.1</td><td>10.1</td><td>21.3</td><td>36.0</td><td>86.8</td></tr><tr><td rowspan="2">Case D</td><td>μ<sub>jISI</sub></td><td>9.48E-03</td><td>9.45E-03</td><td>6.08E-03</td><td>6.07E-03</td><td>4.99E-03</td><td>5.00E-03</td></tr><tr><td>μ<sub>T</sub></td><td>1.2</td><td>2.6</td><td>8.3</td><td>15.8</td><td>32.9</td><td>74.7</td></tr><tr><td rowspan="4">IVA-G-N</td><td rowspan="2">Case A</td><td>μ<sub>jISI</sub></td><td>1.36E-01</td><td>1.69E-01</td><td>7.62E-02</td><td>9.79E-02</td><td>4.76E-02</td><td>4.95E-02</td></tr><tr><td>μ<sub>T</sub></td><td>25.0</td><td>82.8</td><td>76.4</td><td>285.7</td><td>234.7</td><td>1023.5</td></tr><tr><td rowspan="2">Case B</td><td>μ<sub>jISI</sub></td><td>2.19E-02</td><td>2.21E-02</td><td>1.41E-02</td><td>1.41E-02</td><td>1.13E-02</td><td>1.13E-02</td></tr><tr><td>μ<sub>T</sub></td><td>5.9</td><td>21.9</td><td>7.7</td><td>24.3</td><td>17.0</td><td>52.4</td></tr><tr><td rowspan="2">Case C</td><td>μ<sub>jISI</sub></td><td>5.05E-02</td><td>5.54E-02</td><td>2.52E-02</td><td>2.63E-02</td><td>1.69E-02</td><td>1.67E-02</td></tr><tr><td>μ<sub>T</sub></td><td>19.5</td><td>89.5</td><td>36.5</td><td>138.3</td><td>81.1</td><td>304.3</td></tr><tr><td rowspan="2">Case D</td><td>μ<sub>jISI</sub></td><td>9.51E-03</td><td>9.45E-03</td><td>6.08E-03</td><td>6.06E-03</td><td>4.99E-03</td><td>5.00E-03</td></tr><tr><td>μ<sub>T</sub></td><td>1.8</td><td>6.2</td><td>2.3</td><td>6.8</td><td>6.5</td><td>17.6</td></tr></tbody></table>

## VI. CONCLUSION

In this article, we addressed the problem of joint blind source separation, through the IVA-G formulation. First, we derived a cost function, parameterised by the demixing matrices and the precision matrices of the sources, to solve IVA-G based on the maximum likelihood estimator. We then introduced an additional term to fix the scaling ambiguity, hence forcing the precision matrices of output solution to have ones on its diagonal, completing the design of our cost function.

Then, we studied the different terms of this non-convex and non-smooth cost, in particular, we gave explicit formulas for the partial gradients of the smooth term with respect to both blocks of variable and for the proximity operators of the separable terms. Based on these results, we proposed an algorithm adapted from the PALM optimizer and proved that the sequence of iterates converges globally to a critical point. We compared our method to two state-of-the-art IVA-G algorithms and showed that our method is competitive in terms of jISI score and computation time, especially for a moderate number of sources.

As a future work, those encouraging results on synthetic datasets should be confirmed on real data, for instance from fMRI or MEEG. Depending on the application considered, we could leverage model knowledge in the form of new regularization terms in the cost function. A new formulation of the problem could also be used in order to enhance the robustness of our method to errors on the empirical covariance $\hat{R}_x$, in particular to guarantee its validity when working with noisy data or few samples. The proposed proximal alternating algorithm is versatile and would be easily adaptable to a large class of penalties and constraints.
## APPENDIX A

### PROOF OF LEMMA 1

Let $(W, C) \in \mathbb{R}^{N \times N \times K} \times \mathbb{R}^{K \times K \times N}$. The partial derivative of $h$ in (15), with respect to the $(n, m, k)$-th entry of tensor $W$ reads

$$ \frac{\partial h(W, C)}{\partial w_{n,m}^{[k]}} = \frac{1}{2} \sum_{n'=1}^{N} \frac{\partial \operatorname{tr}(C_{n'} W_{n'} \hat{R}_{x} W_{n'}^{\top})}{\partial w_{n,m}^{[k]}}. \quad (35) $$

For $n' \neq n$, the terms of the above sum are null, hence

$$ \frac{\partial h(W, C)}{\partial w_{n,m}^{[k]}} = \frac{1}{2} \frac{\partial \operatorname{tr}(C_n W_n \hat{R}_x W_n^{\top})}{\partial w_{n,m}^{[k]}}. \quad (36) $$

Then, by applying the formula of the derivative of a matrix scalar product [41],

$$ \begin{aligned} \frac{\partial h(W, C)}{\partial w_{n,m}^{[k]}} &= \frac{1}{2} \frac{\partial \langle C_n W_n \hat{R}_x | W_n \rangle}{\partial w_{n,m}^{[k]}} \\ &= \frac{1}{2} \left( \langle \frac{\partial (C_n W_n \hat{R}_x)}{\partial w_{n,m}^{[k]}}, |W_n\rangle + \langle C_n W_n \hat{R}_x | \frac{\partial W_n}{\partial w_{n,m}^{[k]}} \rangle \right) \\ &= \frac{1}{2} \left( \langle C_n \frac{\partial W_n}{\partial w_{n,m}^{[k]}}, \hat{R}_x | W_n \rangle + \langle C_n W_n \hat{R}_x | \frac{\partial W_n}{\partial w_{n,m}^{[k]}} \rangle \right) \\ &= \operatorname{tr} \left( \frac{\partial W_n}{\partial w_{n,m}^{[k]}}, \hat{R}_x W_n^{\top} C_n \right). \end{aligned} \quad (37) $$

Hereabove, $\frac{\partial W_n}{\partial w_{n,m}^{[k]}}$ is a matrix of dimension $K \times KN$ whose elements are equal to 0, except for one, at row index $k$ and column index $(k-1)N+m$, which is equal to 1. We deduce

$$ \frac{\partial h(W, C)}{\partial w_{n,m}^{[k]}} = [C_n W_n \hat{R}_x]_{k, (k-1)N+m}. \quad (38) $$

The above expression can be reexpressed in a more concise matrix form which gives (18), and proves the first part of the lemma.

Let $C \in \mathbb{R}^{K \times K \times N}$. From (18), we can see that $\nabla_W h(\cdot, C)$ is linear, and thus Lipschitz continuous. Moreover, for every $W \in \mathbb{R}^{N \times N \times K}$,

$$ \begin{aligned} \left\| \frac{\partial h(W, C)}{\partial W} \right\|^2 &= \sum_{k=1}^{K} \sum_{n=1}^{N} \sum_{m=1}^{N} \left( |C_n W_n \hat{R}_x|_{k, (k-1)N+m} \right)^2 \\ &= \sum_{n=1}^{N} \sum_{k=1}^{K} \left( \sum_{m=1}^{N} \left( |C_n W_n \hat{R}_x|_{k, (k-1)N+m} \right)^2 \right) \\ &= \sum_{n=1}^{N} \sum_{k=1}^{K} \left\| C_n W_n \hat{R}_x \right\|_{k, (k-1)N+1, \dots, N} \right\|^2 \\ &= \sum_{n=1}^{N} \sum_{k=1}^{K} \left\| c_{n,k} W_n \hat{R}_x^{[k]} \right\|^2. \end{aligned} \quad (39) $$

Using the identity $\|AB\| \le \|A\|_S \|B\|$, we have, for all $k \in \{1, \dots, K\}, n \in \{1, \dots, N\}$,

$$ \|c_{n,k} W_n \hat{R}_x^{[k]}\|^2 \le \|\hat{R}_x^{[k]}\|_S^2 \|c_{n,k} W_n\|^2 \le \varrho_{\hat{R}_x}^2 \|c_{n,k} W_n\|^2. \quad (40) $$

$C_n W_n$ is a $K \times KN$ matrix whose $k$-th line is equal to $c_{n,k} W_n$, hence,

$$ \sum_{k=1}^{K} \|c_{n,k} W_n\|^2 = \|C_n W_n\|^2 \le \|C_n\|_S^2 \|W_n\|^2. \quad (41) $$

Overall,

$$ \begin{aligned} \left\| \frac{\partial h(W, C)}{\partial W} \right\|^2 &\le \sum_{n=1}^{N} \varrho_{\hat{R}_x}^2 \sum_{k=1}^{K} \|c_{n,k} W_n\|^2 \\ &\le \varrho_{\hat{R}_x}^2 \sum_{n=1}^{N} \rho_C^2 \|W_n\|^2 \\ &= \rho_C^2 \varrho_{\hat{R}_x}^2 \|W\|^2. \end{aligned} \quad (42) $$

This concludes the second part of the Lemma, i.e. $W \to h(W, C)$ is Lipschitz differentiable with constant (19).

## APPENDIX B

### PROOF OF LEMMA 2

For all $W \in \mathbb{R}^{N \times N \times K}$, $h(W, \cdot)$ in (15) is quadratic, which yields the expression (22) in a straightforward manner, and concludes the first part of the proof. Let $W' \in \mathbb{R}^{N \times N \times K}$. According to (22), it can be easily checked that $\nabla_C h(W', \cdot)$ is Lipschitz with moduli $\alpha$, hence $C \to h(W', C)$ is Lipschitz differentiable with constant (23), which ends the proof.

## APPENDIX C

### PROOF OF LEMMA 3

Function $f$ in (16), reads, equivalently, for all $W \in \mathbb{R}^{N \times N \times K}$, $f(W) = \sum_{k=1}^{K} \hat{f}(W^{[k]})$, where, for all $M \in \mathbb{R}^{N \times N}$, with singular values $\sigma_M = (\sigma_{M,\ell})_{1 \le \ell \le N}$, $\hat{f}(M) = -\log |\det M| = \phi_f(\sigma_M)$, with

$$ \phi_f: \mathbb{R}^N \mapsto (-\infty, +\infty] \\ \sigma = (\sigma_1, \dots, \sigma_N) \mapsto \begin{cases} +\infty & \text{if } (\exists \ell \in \{1, \dots, N\}) \sigma_\ell \le 0 \\ -\sum_{l=1}^N \log \sigma_\ell & \text{otherwise.} \end{cases} \quad (43) $$

Let $c > 0$ and $W' \in \mathbb{R}^{N \times N \times K}$. Then,

$$ c f(W) + \frac{1}{2} \|W - W'\|^2 = \sum_{k=1}^{K} c \hat{f}(W^{[k]}) + \frac{1}{2} \|W^{[k]} - W'^{[k]}\|^2. $$

Hence, we can minimize the sum by minimizing all its terms independently. In other words, $W \in \operatorname{prox}_{c f}(W')$, if and only if, for every $k \in \{1, \dots, K\}$,

$$ W^{[k]} \in \operatorname{prox}_{c f}(W'^{[k]}). $$

It remains to derive an explicit form for $\operatorname{prox}_{c f}$. To do so, let us use the following lemma whose proof is given in [42] (see also [36], [34] for similar results on the proximity operators of spectral functions).

**Lemma 5** Let $N \ge 1$ and $\phi: \mathbb{R}^N \mapsto (-\infty, +\infty]$ be a function whose domain is not empty and is contained in $[0, +\infty)^N$, such that $\phi$ is invariant by any permutation of its arguments and that the proximity operator of $\phi$ is defined. Let $\Phi: \mathbb{R}^{N \times N} \mapsto (-\infty, +\infty]$, $M \mapsto \phi(\sigma_M)$ where $\sigma_M$ is a vector containing the $N$ singular values of $M$ in any order. Then, the proximity operator of $\Phi$ is defined, and for any $M' \in \mathbb{R}^{N \times N}$, whose singular value decomposition reads $M' = U_{M'} \operatorname{Diag}(\sigma_{M'}) V_{M'}^\top$, if $\sigma_M \in \operatorname{prox}_\phi(\sigma_{M'})$, then we have $U_{M'} \operatorname{Diag}(\sigma_M) V_{M'}^\top \in \operatorname{prox}_\Phi(M')$.
Function (43) is invariant by permutation of its arguments and verifies $\emptyset \neq \text{dom } \phi_f \subset [0, +\infty)^N$. Moreover, for any $\sigma' = (\sigma'_1, \dots, \sigma'_N) \in \mathbb{R}^N$,

$$ 
\begin{aligned}
\operatorname{prox}_{c\phi_f}(\boldsymbol{\sigma}') &= \underset{\boldsymbol{\sigma} \in \mathbb{R}^N}{\operatorname{argmin}} c f(\boldsymbol{\sigma}) + \frac{1}{2} \|\boldsymbol{\sigma} - \boldsymbol{\sigma}'\|^2 \\
&= \underset{\boldsymbol{\sigma} \in (0, +\infty)^N}{\operatorname{argmin}} \sum_{\ell=1}^N \frac{1}{2}(\sigma_\ell - \sigma'_\ell)^2 - c \log \sigma_\ell. \quad (44)
\end{aligned}
$$

The above prox calculation requires to minimize a sum of functions, each term acting on a different variable, which means we can minimize each term independently. Let $l \in \{1, \dots, N\}$, then

$$ 
\begin{aligned}
\forall \sigma_\ell > 0: & \frac{d\left(\frac{1}{2}(\sigma_\ell - \sigma'_\ell)^2 - c \log \sigma_\ell\right)}{d\sigma_\ell} = \sigma_\ell - \sigma'_\ell - c\sigma_\ell^{-1} = 0 \\
\iff & \sigma_\ell^2 - \sigma'_\ell \sigma_\ell - c = 0 \\
\iff & \sigma_\ell \in \left\{ \frac{\sigma_\ell' - \sqrt{(\sigma_\ell')^2 + 4c}}{2}, \frac{\sigma_\ell' + \sqrt{(\sigma_\ell')^2 + 4c}}{2} \right\}.
\end{aligned}
$$

As $\sigma_\ell \mapsto \frac{1}{2}(\sigma_\ell - \sigma'_\ell)^2 - c \log \sigma_\ell$ diverges toward $+\infty$ in $0^+$ and $+\infty$ and is $C^1$ on $(0, +\infty)$, it must have a minimum where its derivative is equal to 0. The only point where it happens is $\sigma_\ell = \frac{\sigma_\ell' + \sqrt{(\sigma_\ell')^2 + 4c}}{2}$, so we can conclude that this point is a global minimum.

Finally, $\operatorname{prox}_{c\phi_f}$ is uniquely defined for any $\sigma' \in \mathbb{R}^N$ and has the explicit form:

$$ \sigma = \operatorname{prox}_{c\phi_f}(\sigma') = \frac{\sigma' + \sqrt{(\sigma')^2 + 4c}}{2} $$

where the square and square root operations are applied component-wise. Applying Lemma 5, with $M = W^{[k]}$, $M' = W'^{[k]}$, for every $k \in \{1, \dots, K\}$, and $\Phi = \hat{f}$, allows to conclude the proof.

## APPENDIX D

### PROOF OF LEMMA 4

For every $C \in \mathbb{R}^{K \times K \times N}$, we can rewrite function (17), as $g(C) = \sum_{n=1}^N \hat{g}(C_n)$ where, for every $M \in \mathbb{R}^{K \times K}$, with singular values $\sigma_M = (\sigma_{M,\ell})_{1\le\ell\le K}$, $\hat{g}(M) = \phi_{\hat{g}}(\sigma_M)$, with

$$ 
\begin{aligned}
\phi_{\hat{g}} : \mathbb{R}^K &\to (-\infty, +\infty] \\
\sigma = (\sigma_1, \dots, \sigma_K) &\mapsto \left\{ \begin{array}{ll} +\infty & \text{if } (\exists \ell \in \{1, \dots, K\}) \sigma_\ell \le 0 \\ -\sum_{l=1}^K \log \sigma_\ell & \text{otherwise.} \end{array} \right.
\end{aligned}
$$

Moreover, for any $\sigma' \in \mathbb{R}^K$,

$$ 
\begin{aligned}
\operatorname{prox}_{c\phi_{\hat{g}}}(\boldsymbol{\sigma}') &= \underset{\boldsymbol{\sigma} \in [\epsilon, +\infty)}{\operatorname{argmin}} \frac{1}{2} \sum_{\ell=1}^{K} (\sigma_{\ell} - \sigma'_{\ell})^2 - c \log \sigma_{\ell} \\
&= \max \left( \epsilon, \frac{\boldsymbol{\sigma}' + \sqrt{(\boldsymbol{\sigma}')^2 + 2c}}{2} \right).
\end{aligned}
$$

Hence the result, applying Lemma 5 for every $n \in \{1, \dots, N\}$, with $\Phi = \hat{g}$, and $(M, M') = (C_n, C'_n)$.

## APPENDIX E

### CONVERGENCE PROOF

a) Lower-boundedness of the cost function: Using Theorem 1,

$$ 
\begin{aligned}
& (\forall (W, C) \in \mathbb{R}^{N \times N \times K} \times \mathbb{R}^{K \times K \times N}) \\
& \qquad J_{\text{IVA-G}}(W, C) \ge \tilde{J}_{\text{IVA-G}}(W) + \frac{KN}{2}.
\end{aligned}
$$

Let us note $\sigma^-$ the smallest eigenvalue of $\hat{R}_x$, then $\hat{R}_x - \sigma^- I_{KN}$ is a symmetric matrix whose eigenvalue are non-negative, i.e., $\hat{R}_x - \sigma^- I_{KN}$ is symmetric positive. Therefore, it is also the case of $W_n(\hat{R}_x - \sigma^- I_{KN})W_n^\top$ for any $n \in \{1, \dots, N\}$, and

$$ W_n \hat{R}_x W_n^\top \succeq \sigma^- W_n W_n^\top, \quad (45) $$

where $\succeq$ denotes the Loewner order relationship defined on $(S_N)^2$ as $(\forall (A, B) \in (S_N)^2) A \succeq B \iff A - B \in S_N^+$.

Yet, $W_n W_n^\top$ is a diagonal matrix whose $k$-th coefficient is $\|w_n^{[k]}\|^2$. Hence $\det(W_n \hat{R}_x W_n^\top) \ge (\sigma^-)^K \det(W_n W_n^\top)$ and

$$ \frac{1}{2} \log \det(W_n \hat{R}_x W_n^\top) \ge \frac{1}{2} K \log(\sigma^-) + \frac{1}{2} \sum_{k=1}^{K} \log \|w_n^{[k]}\|^2. \quad (46) $$

Besides, by Hadamard inequality applied to the lines of $W_n^{[k]}$ for any $k \in \{1, \dots, K\}$, $|\det W_n^{[k]}|^2 \le \prod_{n=1}^N \|w_n^{[k]}\|^2$, hence,

$$ -\log |\det W^{[k]}| \geq -\frac{1}{2} \sum_{n=1}^{N} \log \|w_n^{[k]}\|^2. \quad (47) $$

By summing (46) for $n \in \{1, \dots, N\}$, and (47) for $k \in \{1, \dots, K\}$, it comes:

$$ (\forall W \in \mathbb{R}^{N \times N \times K}) \quad \tilde{J}_{\text{IVA-G}}(W) \geq \frac{KN}{2}(1 + \log(\sigma^{-})). \quad (48) $$

We can conclude that $J_{\text{IVA-G}}$ is lower bounded on $\mathbb{R}^{N \times N \times K} \times \mathbb{R}^{K \times K \times N}$. Finally, since the quadratic regularization term is positive, $J_{\text{IVA-G}}^{\text{Reg}}$ is lower bounded too, which ends the proof.

b) Lipschitz continuity of the gradient: The expressions of the partial gradients of $h$ given in Lemmas 1 and 2 show clearly that $\nabla h$ is a $C^1$ function on $\mathbb{R}^{N \times N \times K} \times \mathbb{R}^{K \times K \times N}$. Consequently, the mean value theorem can be applied to show that it has Lipschitz continuity on any bounded domain.

c) Boundedness of the sequence:

* Boundedness of $(C^{(i)})_{i \in \mathbb{N}}$

First, we notice that (21) defines a norm on $\mathbb{R}^{K \times K \times N}$. For all $(i, j) \in \mathbb{N}^2$ such that $C^{(i,j)}$ belongs to the sequence generated by PALM-IVA-G, let us note:

$$ 
\begin{aligned}
\bar{C}_n^{(i,j)} &= C_n^{(i,j)} - c_C(\alpha(\operatorname{Diag}(C_n^{(i,j)}) - I_K) \\
&\qquad - \frac{1}{2} W_n^{(i+1)} \hat{R}_x W_n^{(i+1)\top}). \\
\|C_n^{(i,j+1)}\|_S &= \max \left( \epsilon, \frac{\|\bar{C}_n^{(i,j)}\|_S + \sqrt{\|\bar{C}_n^{(i,j)}\|_S^2 + 2c_C}}{2} \right) \\
&\le \|\bar{C}_n^{(i,j)}\|_S + \sqrt{\frac{c_C}{2}},
\end{aligned}
\quad (49) $$
assuming that $\epsilon < \sqrt{\frac{\epsilon c}{2}}$, which is verified in our experimental settings.

As $\mathbf{W}_n^{(i+1)} \hat{\mathbf{R}}_x \mathbf{W}_n^{(i+1)\top}$ is symmetric positive,

$$ \begin{aligned} \|\overline{\mathbf{C}}_n^{(i,j)}\|_{\mathrm{S}} &\le \|\mathbf{C}_n^{(i,j)} - \gamma_c(\mathrm{Diag}(\mathbf{C}_n^{(i,j)}) - \mathbf{I}_K)\|_{\mathrm{S}} \\ &= \|\mathbf{C}_n^{(i,j)} - \gamma_c \mathrm{Diag}(\mathbf{C}_n^{(i,j)})\|_{\mathrm{S}} + \gamma_c. \end{aligned} \quad (50) $$

Let $\mathbf{u} \in \mathbb{R}^K$ be a unit norm vector. Using the inequality for matrices of $S_K^+$, $\forall (k,l) \in \{1, \dots, K\} |c_{n,k,l}^{(i,j)}| \le \sqrt{c_{n,k,k}^{(i,j)} \sqrt{c_{n,l,l}^{(i,j)}}}$, it comes

$$ \begin{aligned} \mathbf{u}^{\top} \mathbf{C}_{n}^{(i,j)} \mathbf{u} &= \sum_{1 \le k,l \le K} c_{n,k,l}^{(i,j)} u_k u_l \le \sum_{1 \le k,l \le K} |c_{n,k,l}^{(i,j)}| |u_k| |u_l| \\ &\le \sum_{1 \le k,l \le K} \sqrt{c_{n,k,k}^{(i,j)}} |u_k| \sqrt{c_{n,l,l}^{(i,j)}} |u_l| \\ &\le \sum_{1 \le k,l \le K} \frac{c_{n,k,k}^{(i,j)} u_k^2}{2} + \frac{c_{n,l,l}^{(i,j)} u_l^2}{2} = K \mathbf{u}^{\top} \mathrm{Diag}(\mathbf{C}_n^{(i,j)}) \mathbf{u}. \end{aligned} $$

It follows that, for all $\mathbf{u} \in \mathbb{R}^K$ such that $\|\mathbf{u}\| = 1$, $\mathbf{u}^{\top}(\mathbf{C}_n^{(i,j)} - \gamma_c \mathrm{Diag}(\mathbf{C}_n^{(i,j)}))\mathbf{u} \le (1 - \frac{\gamma_c}{K})\mathbf{u}^{\top}\mathbf{C}_n^{(i,j)}\mathbf{u} \le (1 - \frac{\gamma_c}{K})\|\mathbf{C}_n^{(i,j)}\|_{\mathrm{S}}$. Combining this result with (49) and (50), we deduce that $\|\mathbf{C}_n^{(i,j+1)}\|_{\mathrm{S}} \le (1 - \frac{\gamma_c}{K})\|\mathbf{C}_n^{(i,j)}\|_{\mathrm{S}} + \gamma_c + \sqrt{\frac{\epsilon c}{2}}$, then we can prove by recurrence that for all $(i, j)$, $\|\mathbf{C}_n^{(i,j)}\|_{\mathrm{S}} \le \max (\|\mathbf{C}_n^{(0)}\|_{\mathrm{S}}, K(1 + \sqrt{\frac{1}{2\alpha\gamma_c}}))$. Consequently, if we define

$$ \bar{\rho} = \max (\|\rho_{C^{(0)}}\|_{\mathrm{S}}, K(1 + \sqrt{\frac{1}{2\alpha\gamma_c}})), \quad (51) $$

$\bar{\rho}$ is an upper bound for $(\rho_{C^{(i)}})_{i \in \mathbb{N}}$, which proves that $(C^{(i)})_{i \in \mathbb{N}}$ is bounded.

• Boundedness of $(\mathcal{W}^{(i)})_{i \in \mathbb{N}}$

Again, we fix $i \in \mathbb{N}$, the proximal operator used in Algorithm 1 ensures that $\mathcal{C}^{(i)} \in (\epsilon \mathbf{I}_K + S_K^+)^N$. Thus, $(\forall n \in \{1, \dots, N\}) \mathrm{tr}(\mathbf{C}_n^{(i)} \mathbf{W}_n^{(i)} \hat{\mathbf{R}}_x \mathbf{W}_n^{(i)\top}) \ge \epsilon \mathrm{tr}(\mathbf{W}_n^{(i)} \hat{\mathbf{R}}_x \mathbf{W}_n^{(i)\top}) \ge \epsilon \sigma^- \|\mathbf{W}_n^{(i)}\|^2$ using (45). By summing these inequalities,

$$ \frac{1}{2} \sum_{n=1}^{N} \mathrm{tr}(\mathbf{C}_n^{(i)} \mathbf{W}_n^{(i)} \hat{\mathbf{R}}_x \mathbf{W}_n^{(i)\top}) \ge \frac{\epsilon \sigma^-}{2} \|\mathbf{W}^{(i)}\|^2. \quad (52) $$

Besides, let us consider the concave inequality on the logarithm function at point $\frac{2}{\epsilon\sigma^-} > 0$,

$$ (\forall x \in (0, +\infty)) \quad \log x \le \log \frac{2}{\epsilon\sigma^-} + \frac{\epsilon\sigma^-}{2}x - 1. \quad (53) $$

Taking $k \in \{1, \dots, K\}$ and summing (53) applied on $\|\mathbf{w}_n^{(i)[k]}\|$ for $n \in \{1, \dots, N\}$, we obtain

$$ \sum_{n=1}^{N} \log \|\mathbf{w}_n^{(i)[k]}\| \le N(\log \frac{2}{\epsilon\sigma^-} - 1) + \frac{\epsilon\sigma^-}{2} \|\mathbf{W}^{(i)[k]}\|^2. \quad (54) $$

And combining this inequality with (46) gives

$$ -\log |\det \mathbf{W}^{(i)[k]}| \ge -\frac{\epsilon\sigma^-}{4} \|\mathbf{W}^{(i)[k]}\|^2 - \frac{N}{2} (\log \frac{2}{\epsilon\sigma^-} - 1). $$

Finally, summing for $k \in \{1, \dots, K\}$,

$$ -\sum_{k=1}^{K} \log |\det \mathbf{W}^{(i)[k]}| \ge -\frac{\epsilon\sigma^-}{4} \|\mathbf{W}^{(i)}\|^2 - \frac{KN}{2} (\log \frac{2}{\epsilon\sigma^-} - 1). \quad (55) $$

Summing (55) and (52) yields

$$ \begin{aligned} & \frac{\epsilon\sigma^-}{4} \|\mathbf{W}^{(i)}\|^2 - \frac{KN}{2}(\log \frac{2}{\epsilon\sigma^-} - 1) \\ &\le \frac{1}{2} \sum_{n=1}^{N} \mathrm{tr}(\mathbf{C}_n^{(i)} \mathbf{W}_n^{(i)} \hat{\mathbf{R}}_x \mathbf{W}_n^{(i)\top}) - \sum_{k=1}^{K} \log |\det \mathbf{W}^{(i)[k]}| \\ &= J_{\mathrm{IVA-G}}(\mathcal{W}^{(i)}, \mathcal{C}^{(i)}) + \frac{1}{2} \sum_{n=1}^{N} \log \det \mathbf{C}_n^{(i)} \\ &\le J_{\mathrm{IVA-G}}(\mathcal{W}^{(i)}, \mathcal{C}^{(i)}) + \frac{KN}{2} \log \bar{\rho}. \end{aligned} $$

Following the proof of [29, Lemma 2], it is sufficient that the first four items of Assumption 1 hold to obtain that $J_{\mathrm{IVA-G}}^{\mathrm{Reg}}(\mathcal{W}^{(i)}, \mathcal{C}^{(i)})$ is a non-increasing sequence. We thus have

$$ \|\mathcal{W}^{(i)}\|^2 \le \frac{4}{\epsilon\sigma^-} \left( J_{\mathrm{IVA-G}}^{\mathrm{Reg}}(\mathcal{W}^{(0)}, \mathcal{C}^{(0)}) + \frac{KN}{2} \log \bar{\rho} \left( \frac{2}{\epsilon\sigma^-} - 1 \right) \right). \quad (56) $$

This shows that $(\mathcal{W}^{(i)})_{i \in \mathbb{N}}$ is bounded too.

d) Bounds for the corrected Lipschitz moduli: Using (19), we can set $\lambda_W^- = \epsilon\rho\hat{\mathbf{R}}_x$ and $\lambda_W^+ = \bar{\rho}\rho\hat{\mathbf{R}}_x$. We can also set $\lambda_C^- = \lambda_C^+ = \alpha$.

e) KL property: KL property is a key tool from functional analysis to demonstrate convergence of iterates in the non-convex setting [29]. It is sufficient here to apply the result from [43, Section 4.3], which states that a function verifies KL on the domain of its subdifferential, if it is lower semi-continuous, proper, and definable in an o-minimal structure.

We rely on the o-minimal structure $\mathcal{S}(\mathbb{R}_{\mathrm{an,exp}})$, defined in [44, Section 2, Example (6)], which contains log, exp and all semi-algebraic functions. The properties in [44, Section 5] can then be used to show that $\mathcal{S}(\mathbb{R}_{\mathrm{an,exp}})$ contains our cost function $J_{\mathrm{IVA-G}}^{\mathrm{Reg}}$, after identifying $\mathbb{R}^{N \times N \times K} \times \mathbb{R}^{K \times K \times N}$ with $\mathbb{R}^{KN(K+N)}$.

## APPENDIX F

### EFFECT OF $\epsilon$

We provide here an empirical analysis of the role of the constraint term $\sum_{n=1}^N \iota_{[\epsilon, +\infty)^K}(\sigma_{C_n})$, in our cost function. We ran 10 simulations in Case B for $(K, N) = (10, 10)$ and $\epsilon \in \{10^{-12}, 10^{-9}, 10^{-6}, 10^{-3}, 10^{-2}, 10^{-1}\}$. The boxplots of the jISI scores and computation times obtained are presented in 5. As one can notice, the results are almost identical for all these values of $\epsilon$, with a slight increase of time for $\epsilon = 10^{-1}$. This proves that taking $\epsilon = 10^{-12}$ as we did is extremely conservative and should ensure that the returned $(\mathcal{W}^*, \mathcal{C}^*)$ is optimal, in the ML estimator.

## APPENDIX G

### ROBUSTNESS OF PALM-IVA-G TO NOISE AND DOWN SAMPLING

Whereas the experiments were done with a relatively large number of samples ($V = 10000$) and with no noise, dealing
<figure><img src="image_9.png" /><figcaption>Fig. 5: Distributions of the jISI scores and computation times for various values of $\epsilon$ on 10 runs, $(K, N) = (10, 10)$.</figcaption></figure>

with real data may be less than ideal, and lead to a discrepancy between the empirical covariance of the data $\hat{\mathbf{R}}_x$ and the real covariance of the random vector $\mathbf{x}$, which in turn can affect the quality of the separation. This section presents two additional comparisons to assess performance in these scenarios: (i) a reduction of the number of samples, and (ii) the presence of a white noise over the observed data, here simulated by adding a term $a\mathbf{I}_{NK}$ to $\hat{\mathbf{R}}_x$.

We ran 100 simulations in Case B for $(K, N) = (10, 10)$, first varying $V \in \{10000, 5000, 1000, 500, 200, 150, 120, 100, 50\}$, and then varying $a \in \{0, 10^{-3}, 10^{-2}, 10^{-1}, 1, 10\}$ with $V = 10000$ samples. The results are presented, for a maximum number of iterations to $\bar{N} = 5000$ for PALM-IVA-G and $\bar{N} = 20000$ for IVA-G-N and IVA-G-V.

### A. Robustness to downsampling

We observe a steady increase of the jISI score as we lower $V$. Down to $V = 500$, the separation is still acceptable, but for 200 samples and less, the algorithm is unable to separate the sources. On the computation times, we also notice an overall increasing difficulty to reach convergence as the number of samples decreases, with a failure to converge when $V \le 100$. For PALM-IVA-G in particular, the limit of 40 seconds corresponding to the number of iterations meeting its upper bound of 5000 is almost systematically reached, which is consistent with the fact that $\hat{\mathbf{R}}_x$ is singular in this case for these dimensions of the problem. This behavior was expected, and is observed for all competing methods based on IVA-G model.

In conclusion, it is recommended to run PALM-IVA-G (and other IVA-G based techniques) with $V \gg KN$ in practice.

### B. Robustness to noise

The jISI score and runtime, for various values of $a$ are displayed in Figure 7. This shows that we can add white noise with a variance up to $a = 1$ to the data before observing a significant degradation in the jISI score and the computation time, which is an excellent marker of the practical robustness of our method.

### ACKNOWLEDGMENT

The work of the authors is funded by the European Research Council Starting Grant MAJORIS ERC-2019-STG-850925 and the U.S. National Science Foundation Grant NSF 2316420.

<figure><img src="image_10.png" /><figcaption>Fig. 6: Distributions of the jISI scores and computation times for various values of $V$ on 100 runs, $(K, N) = (10, 10)$. (top) PALM-IVA-G, (middle) IVA-G-V, (bottom) IVA-G-N.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 7: Distributions of the jISI scores and computation times for various values of $a$ on 100 runs, $(K, N) = (10, 10)$.</figcaption></figure>

### REFERENCES

[1] V. Calhoun and T. Adali, “Unmixing fMRI with independent component analysis,” *IEEE Engineering in Medicine and Biology Magazine*, vol. 25, no. 2, pp. 79–90, Mar. 2006.

[2] Z. Luo, C. Li, and L. Zhu, “A comprehensive survey on blind source separation for wireless adaptive processing: Principles, perspectives, challenges and new research directions,” *IEEE Access*, vol. 6, pp. 66685–66708, 2018.

[3] M. Castella, J.-C. Pesquet, and A. Petropulu, “A family of frequency- and time-domain contrasts for blind separation of convolutive mixtures of temporally dependent signals,” *IEEE Transactions on Signal Processing*, vol. 53, no. 1, pp. 107–120, 2005.

[4] A. Nielsen, “Multiset canonical correlations analysis and multispectral, truly multitemporal remote sensing data,” *IEEE Transactions on Image Processing*, vol. 11, no. 3, pp. 293–305, Mar. 2002, conference Name: IEEE Transactions on Image Processing.

[5] I. Santamaria, “Handbook of blind source separation: independent component analysis and applications (Comon, P. and Jutten, C. ; 2010 [Book
Review])," *IEEE Signal Processing Magazine*, vol. 30, no. 2, pp. 133-134, 2013.

[6] T. Adali, M. Anderson, and G.-S. Fu, "Diversity in independent component and vector analyses: Identifiability, algorithms, and applications in medical imaging," *IEEE Signal Proc. Mag.*, vol. 31, no. 3, pp. 18–33, May 2014.

[7] V. D. Calhoun and T. Adali, "Multisubject independent component analysis of fMRI: A decade of intrinsic networks, default mode, and neurodiagnostic discovery," *IEEE Reviews in Biomedical Engineering*, vol. 5, pp. 60–73, 2012.

[8] R. García-Bermúdez, F. R. Ruiz, J. G. Peñalver, O. V. Cansino, L. V. Pérez, C. Torres, and R. Becerra-García, "Evaluation of electro-oculography data for ataxia sca-2 classification: A blind source separation approach," in *2010 10th International Conference on Intelligent Systems Design and Applications*, 2010, pp. 237–241.

[9] B. Gabrielson, M. A. B. S. Akhonda, S. Bhinge, J. Brooks, Q. Long, and T. Adali, "Joint-IVA for identification of discriminating features in EEG: Application to a driving study," *Biomedical Signal Processing and Control*, vol. 61, p. 101948, Aug. 2020.

[10] D. Lahat, T. Adali, and C. Jutten, "Multimodal data fusion: An overview of methods, challenges, and prospects," *Proceedings of the IEEE*, vol. 103, no. 9, pp. 1449–1477, 2015.

[11] Z. Boukouvalas, M. Puerto, D. C. Elton, P. W. Chung, and M. D. Fuge, "Independent vector analysis for molecular data fusion: Application to property prediction and knowledge discovery of energetic materials," in *2020 28th European Signal Processing Conference (EUSIPCO)*. Amsterdam, Netherlands: IEEE, Jan. 2021, pp. 1030–1034.

[12] M. Anderson, X.-L. Li, and T. Adali, "Joint blind source separation with multivariate Gaussian model: Algorithms and performance analysis," *IEEE Trans. Signal Processing*, vol. 60, no. 4, pp. 2049–2055, April 2012.

[13] J.-C. Pesquet and E. Moreau, "Cumulant-based independence measures for linear mixtures," *IEEE Transactions on Information Theory*, vol. 47, no. 5, pp. 1947–1956, 2001.

[14] P. Comon, "Independent component analysis, a new concept?" *Signal processing*, vol. 36, no. 3, pp. 287–314, 1994.

[15] T. Kim, T. Eltoft, and T.-W. Lee, "Independent vector analysis: An extension of ica to multivariate components," in *International conference on independent component analysis and signal separation*. Springer, 2006, pp. 165–172.

[16] A. Hiroe, "Solution of permutation problem in frequency domain ica, using multivariate probability density functions," in *Independent Component Analysis and Blind Signal Separation: 6th International Conference, ICA 2006, Charleston, SC, USA, March 5-8, 2006. Proceedings 6*. Springer, 2006, pp. 601–608.

[17] T. Kim, H. T. Attias, S.-Y. Lee, and T.-W. Lee, "Blind source separation exploiting higher-order frequency dependencies," *IEEE transactions on audio, speech, and language processing*, vol. 15, no. 1, pp. 70–79, 2006.

[18] M. Anderson, G.-S. Fu, R. Phlypo, and T. Adali, "Independent vector analysis: Identification conditions and performance bounds," *IEEE Transactions on Signal Processing*, vol. 62, no. 17, pp. 4399–4410, 2014.

[19] T. Adali, M. Akhonda, and V. D. Calhoun, "ICA and IVA for data fusion: An overview and a new approach based on disjoint subspaces," *IEEE sensors letters*, vol. 3, no. 1, p. 7100404, Jan. 2019.

[20] V. D. Calhoun and T. Adali, "Feature-based fusion of medical imaging data," *IEEE Transactions on Information Technology in Biomedicine*, vol. 13, no. 5, pp. 711–720, 2009.

[21] A. Michael, M. Anderson, R. Miller, T. Adali, and V. Calhoun, "Preserving subject variability in group fMRI analysis: Performance evaluation of GICA vs. IVA," *Frontiers in systems neuroscience*, vol. 8, p. 106, 06 2014.

[22] J. Laney, K. P. Westlake, S. Ma, E. Woytowicz, V. D. Calhoun, and T. Adali, "Capturing subject variability in fMRI data: A graph-theoretical analysis of GICA vs. IVA," *Journal of neuroscience methods*, vol. 247, pp. 32–40, 2015.

[23] S. Bhinge, R. Mowakea, V. D. Calhoun, and T. Adali, "Extraction of time-varying spatiotemporal networks using parameter-tuned constrained iva," *IEEE transactions on medical imaging*, vol. 38, no. 7, pp. 1715–1725, 2019.

[24] N. Ono, "Stable and fast update rules for independent vector analysis based on auxiliary function technique," in *2011 IEEE Workshop on Applications of Signal Processing to Audio and Acoustics (WASPAA)*. IEEE, 2011, pp. 189–192.

[25] R. Scheibler and N. Ono, "Fast and stable blind source separation with rank-1 updates," in *ICASSP 2020-2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)*. IEEE, 2020, pp. 236–240.

[26] —, "Mm algorithms for joint independent subspace analysis with application to blind single and multi-source extraction," *arXiv preprint arXiv:2004.03926*, 2020.

[27] R. Scheibler, "Independent vector analysis via log-quadratically penalized quadratic minimization," *IEEE Transactions on Signal Processing*, vol. 69, pp. 2509–2524, 2021.

[28] K. Yatabe and D. Kitamura, "Time-frequency-masking-based determined bss with application to sparse iva," in *ICASSP 2019-2019 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)*. IEEE, 2019, pp. 715–719.

[29] J. Bolte, S. Sabach, and M. Teboulle, "Proximal alternating linearized minimization for nonconvex and nonsmooth problems," *Mathematical Programming*, vol. 146, no. 1, pp. 459–494, 2014.

[30] T. Kim, I. Lee, and T.-W. Lee, "Independent vector analysis: Definition and algorithms," in *Proceedings of the 40th Asilomar Conference on Signals, Systems and Computers (ASILOMAR 2006)*, 2006, pp. 1393–1396.

[31] P. Ciuciu, G. Varoquaux, P. Abry, S. Sadaghiani, and A. Kleinschmidt, "Scale-free and multifractal time dynamics of fmri signals during rest and task," *Frontiers in Physiology*, vol. 3, p. 186, 06 2012.

[32] T. Vu, F. Laport, H. Yang, V. D. Calhoun, and T. Adali, "Constrained independent vector analysis with reference for multi-subject fmri analysis," *IEEE Transactions on Biomedical Engineering*, vol. 71, no. 12, pp. 3531–3542, 2024.

[33] E. Chouzenoux, T. T.-K. Lau, C. Lefort, and J.-C. Pesquet, "Optimal multivariate Gaussian fitting with applications to PSF modeling in two-photon microscopy imaging," *Journal of Mathematical Imaging and Vision*, vol. 61, no. 7, pp. 1037–1050, Sept. 2019.

[34] H. H. Bauschke and P. L. Combettes, *Convex Analysis and Monotone Operator Theory in Hilbert Spaces*, 2nd ed. New York: Springer, 2017.

[35] P. L. Combettes and J.-C. Pesquet, "Proximal splitting methods in signal processing," in *Fixed-Point Algorithms for Inverse Problems in Science and Engineering*, Bauschke, H. Burachik, R. Combettes, P. Elser, V. Luke, D. Wolkowicz, and H. (Eds.), Eds. Springer, 2011, pp. 185–212.

[36] A. Benfenati, E. Chouzenoux, and J.-C. Pesquet, "Proximal approaches for matrix optimization problems: Application to robust precision matrix estimation." *Signal Processing*, no. 69, p. 107417, 2020.

[37] E. Chouzenoux, J.-C. Pesquet, and A. Repetti, "A block coordinate variable metric forward-backward algorithm," *Journal of Global Optimization*, vol. 66, no. 3, pp. 457–485, 2016.

[38] L. T. K. Hien, D. Phan, and N. Gillis, "An inertial block majorization minimization framework for nonsmooth nonconvex optimization," *Journal of Machine Learning Research*, vol. 24, pp. 1–41, 2023.

[39] E. Chouzenoux, M.-C. Corbineau, J.-C. Pesquet, and G. Scrivanti, "A variational approach for joint image recovery and feature extraction based on spatially varying generalised Gaussian models," *Journal of Mathematical Imaging and Vision*, pp. 1–22, 2024.

[40] S. Amari, A. Cichocki, and H. Yang, "A new learning algorithm for blind signal separation," in *Advances in Neural Information Processing Systems*, vol. 8. MIT Press, 1995.

[41] K. B. Petersen and M. S. Pedersen, "The matrix cookbook," Oct. 2008, version 20081110.

[42] J.-C. Pesquet, "Prox of a function of the singular values of a matrix," internal report, CVN, CentraleSupélec, 2022.

[43] H. Attouch, J. Bolte, P. Redont, and A. Soubeyran, "Proximal alternating minimization and projection methods for nonconvex problems: An approach based on the Kurdyka-Lojasiewicz inequality," *Mathematics of Operations Research*, vol. 35, no. 2, pp. 438–457, 2010.

[44] L. van den Dries and C. Miller, "Geometric categories and o-minimal structures," *Duke Mathematical Journal*, vol. 84, no. 2, Aug. 1996.