# Leveraging Joint-Diagonalization in Transform-Learning NMF

Sixin Zhang, , Emmanuel Soubies, , and Cédric Févotte The authors are with
IRIT, Université de Toulouse, CNRS, Toulouse, France. Email:
<EMAIL_ADDRESS>This work was supported by the European Research
Council (ERC FACTORY-CoG-6681839), the French National Research Agency (ANR
ANITI-3IA) and the National Research Foundation, Prime Minister’s Office,
Singapore under its Campus for Research Excellence and Technological
Enterprise (CREATE) programme.

###### Abstract

Non-negative matrix factorization with transform learning (TL-NMF) is a recent
idea that aims at learning data representations suited to NMF. In this work,
we relate TL-NMF to the classical matrix joint-diagonalization (JD) problem.
We show that, when the number of data realizations is sufficiently large, TL-
NMF can be replaced by a two-step approach—termed as JD+NMF—that estimates the
transform through JD, prior to NMF computation. In contrast, we found that
when the number of data realizations is limited, not only is JD+NMF no longer
equivalent to TL-NMF, but the inherent low-rank constraint of TL-NMF turns out
to be an essential ingredient to learn meaningful transforms for NMF.

###### Index Terms:

Transform learning, Nonnegative matrix factorization, Joint-diagonalization,
Statistical estimation, Nonconvex optimization, Quasi-Newton method.

##  I Introduction

Let ${\mathbf{Y}}\in{\mathbb{R}}^{M\times N}$ be a random matrix that
satisfies the following moment conditions: $\forall(m^{\prime},m,n)$,

$\left\\{\begin{array}[]{ll}{\mathbb{E}}([{\bar{\bm{\Phi}}}{\mathbf{Y}}]_{mn})=0,&\text{(mean)}\\\
{\mathbb{E}}([{\bar{\bm{\Phi}}}{\mathbf{Y}}]_{mn}^{2})=[\bar{\mathbf{W}}\bar{\mathbf{H}}]_{mn},&\text{(variance)}\\\
{\mathbb{E}}([{\bar{\bm{\Phi}}}{\mathbf{Y}}]_{mn}[{\bar{\bm{\Phi}}}{\mathbf{Y}}]_{m^{\prime}n})=0,&\text{(covariance)}\end{array}\right.$
(1)

where ${\bar{\bm{\Phi}}}\in{\mathbb{R}}^{M\times M}$ is an orthogonal matrix
and $\bar{\mathbf{W}}\in{\mathbb{R}}_{+}^{M\times\bar{K}}$,
$\bar{\mathbf{H}}\in{\mathbb{R}}_{+}^{\bar{K}\times N}$ are two nonnegative
matrices such that $[\bar{\mathbf{W}}\bar{\mathbf{H}}]_{mn}>0$. Here $\bar{K}$
denotes the non-negative rank of $\bar{\mathbf{W}}\bar{\mathbf{H}}$ [1]. These
fairly general conditions encompass, for instance, the popular Gaussian
composite model (GCM) [2, 3], which reads as

$[{\bar{\bm{\Phi}}}{\mathbf{Y}}]_{mn}\underset{\text{i.i.d.}}{\sim}\mathcal{N}(0,[\bar{\mathbf{W}}\bar{\mathbf{H}}]_{mn}),\;\forall(m,n),$
(2)

where $\mathcal{N}$ denotes the real-valued Gaussian distribution.

In this paper, we are concerned with the problem of estimating the triplet
$({\bar{\bm{\Phi}}},\bar{\mathbf{W}},\bar{\mathbf{H}})$ from $S$ i.i.d.
realizations of ${\mathbf{Y}}$, which we denote by
$\\{{\mathbf{Y}}^{(s)}\\}_{s=1}^{S}$ (this includes the most frequent scenario
$S=1$). In the literature, this problem is referred to as transform-learning
nonnegative matrix factorization (TL-NMF) [4, 5] and is also a special case of
independent low-rank tensor analysis (ILRTA) [6]. When the columns of
${\mathbf{Y}}$ are composed of overlapping segments of a temporal signal and
${\bar{\bm{\Phi}}}$ is a fixed frequency-transform matrix such as the DCT (we
only consider real-valued transforms in this paper), the product
${\bar{\bm{\Phi}}}{\mathbf{Y}}$ defines a short-time frequency transform. Its
power values $[{\bar{\bm{\Phi}}}{\mathbf{Y}}]_{mn}^{2}$ define a
_spectrogram_. Equation (1) dictates that the spectrogram is given in
expectation by the nonnegative matrix factorization (NMF)
$\bar{\mathbf{W}}\bar{\mathbf{H}}$. The zero-mean assumption is valid for
centered temporal signals such as audio signals. The zero-covariance
assumption is a working assumption that reflects uncorrelation of the
transform coefficients. The dictionary matrix $\bar{\mathbf{W}}$ captures
meaningful spectral patterns while the activation matrix $\bar{\mathbf{H}}$
describes how the spectral samples are decomposed onto the dictionary.

NMF has found a large range of applications in audio signal processing, such
as source separation or music transcription [3, 7, 8, 9]. In TL-NMF, the fixed
transform assumption is relaxed and ${\bar{\bm{\Phi}}}$ is estimated together
with $\bar{\mathbf{W}}$ and $\bar{\mathbf{H}}$ [4]. This was shown to lead in
some cases to more meaningful representations. For example, it was shown in
[6] that better source separation performance can be achieved with adaptive
learned transform. In [5], it was shown that learning ${\bar{\bm{\Phi}}}$
allows to capture harmonic atoms with fundamental frequencies that precisely
match the music notes in the data, when a DCT transform can only abide by a
fixed frequency grid. TL-NMF was inspired by works about learning sparsifying
transforms [10] and to a lesser extent independent component analysis (ICA)
[11, 12].

TL-NMF can be tackled through the resolution of the following optimization
problem

$({\bm{\Phi}}^{\ast},{\mathbf{W}}^{\ast},{\mathbf{H}}^{\ast})\in\operatorname*{arg\,min}_{\begin{subarray}{c}{\bm{\Phi}}\in\mathcal{O}(M)\\\
({\mathbf{W}},{\mathbf{H}})\in
F_{K}\end{subarray}}C_{S}({\bm{\Phi}},{\mathbf{W}},{\mathbf{H}}),$ (3)

where $\mathcal{O}(M)$ denotes the real orthogonal matrix group on
${\mathbb{R}}^{M}$ and $F_{K}\subset{\mathbb{R}}^{M\times
K}_{+}\times{\mathbb{R}}^{K\times N}_{+}$ is a constraint set for the
nonnegative matrices ${\mathbf{W}}$ and ${\mathbf{H}}$ which will be made
explicit later. Here $K$ represents the factorization rank whose choice
depends on the considered dataset and application [9, 13]. The objective
function $C_{S}$ in (3) is defined by

$C_{S}({\bm{\Phi}},{\mathbf{W}},{\mathbf{H}})=\sum_{m,n=1}^{M,N}\frac{{\mathbb{E}}_{S}([{\bm{\Phi}}{\mathbf{Y}}]_{mn}^{2})+{\epsilon_{0}}}{[{\mathbf{W}}{\mathbf{H}}]_{mn}+{\epsilon_{0}}}\\\
+\log\left([{\mathbf{W}}{\mathbf{H}}]_{mn}+{\epsilon_{0}}\right)$ (4)

where ${\mathbb{E}}_{S}([{\bm{\Phi}}{\mathbf{Y}}]_{mn}^{2})$ is a consistent
estimator of ${\mathbb{E}}([{\bm{\Phi}}{\mathbf{Y}}]_{mn}^{2})$ computed from
the $S$ i.i.d. realizations $\\{{\mathbf{Y}}^{(s)}\\}_{s=1}^{S}$, and
${\epsilon_{0}}\geq 0$.

The rationale behind the choice of Problem (3) is that, when

  * i)

$S=\infty$ (i.e.,
${\mathbb{E}}_{S}([{\bm{\Phi}}{\mathbf{Y}}]_{mn}^{2})={\mathbb{E}}([{\bm{\Phi}}{\mathbf{Y}}]_{mn}^{2})$),
and ${\epsilon_{0}}=0$,

  * ii)

$K\geq\bar{K}$ (i.e., the NMF rank $K$ in (3) is larger than $\bar{K}$, the
true rank in model (1)),

each global minimizer
$({\bm{\Phi}}^{\ast},{\mathbf{W}}^{\ast},{\mathbf{H}}^{\ast})$ of $C_{S}$ is
such that the rows of ${\bm{\Phi}}^{\ast}$ and ${\bar{\bm{\Phi}}}$ span the
same subspaces (we say that ${\bar{\bm{\Phi}}}$ is identifiable by solving
(3)) and
${\mathbf{W}}^{\ast}{\mathbf{H}}^{\ast}=\bar{\mathbf{W}}\bar{\mathbf{H}}$ [5].
Moreover, under the GCM, $C_{S}$ corresponds to the expected negative log-
likelihood of ${\mathbf{Y}}$ conditioned to
$({\bm{\Phi}},{\mathbf{W}},{\mathbf{H}})$.

##### Contributions and Outline

This identifiability result relies on a key property of $C_{S}$ that is, in
the regime $S=\infty$ and $K\geq\bar{K}$,

$\forall{\bm{\Phi}}\in\mathcal{O}(M),\;\min_{({\mathbf{W}},{\mathbf{H}})\in
F_{K}}C_{S}({\bm{\Phi}},{\mathbf{W}},{\mathbf{H}})=L_{S}({\bm{\Phi}}),$ (5)

where $L_{S}$ is related to the joint-diagonalization (JD) criterion derived
in [14] (and formally defined in Section II).

This suggests the following two-step alternative to Problem (3)

$\begin{array}[]{rl}{\bm{\Phi}}^{\bullet}\in&\displaystyle\mkern-3.0mu\operatorname*{arg\,min}_{{\bm{\Phi}}\in\mathcal{O}(M)}\,L_{S}({\bm{\Phi}}),\\\
({\mathbf{W}}^{\bullet},{\mathbf{H}}^{\bullet})\in&\displaystyle\mkern-15.0mu\operatorname*{arg\,min}_{({\mathbf{W}},{\mathbf{H}})\in
F_{K}}C_{S}({\bm{\Phi}}^{\bullet},{\mathbf{W}},{\mathbf{H}}),\end{array}$ (6)

which we refer to as JD+NMF. Although it is straightforward from (5) that
problems (3) and (6) are equivalent in the ideal case $S=\infty$ and
${\epsilon_{0}}=0$, their relation in the more practical situation $S<\infty$
and ${\epsilon_{0}}>0$ (usually required for numerical stability) needs to be
analyzed. In particular:

  * •

How does (6) deviate from (3) with respect to (w.r.t.) $S$?

  * •

Is one of these two formulations preferable from an optimization perspective?
(quality of the reached local minimizer, convergence speed)

The main contribution of this work is to provide theoretical and numerical
insights to these questions. In Section II, we set out the basic assumptions
and definitions required for our analysis and prove in Theorem 1 that both TL-
NMF (Problem (3)) and JD+NMF (Problem (6)) admit at least one solution. Then,
Section III is dedicated to the characterization of the closeness between the
solution sets of these two problems. We first identify in Theorem 2 quantities
that allow to characterize three regimes where the solution sets of TL-NMF and
JD+NMF are i) disjoint, ii) partially intersecting, and iii) equal,
respectively. Then, by analyzing the asymptotic behavior (w.r.t. $S$) of these
key quantities, we show that with high probability, the gaps between the two
solution sets converge to zero (in terms of the Itakua-Saito divergence [15,
2]) as $S$ grows (Proposition 2 and Theorem 3). Under the GCM, we further
prove in Theorem 4 that the rate of convergence is at least of the order of
$O(1/S)$. In Section IV, we adapt existing TL-NMF and JD minimization
algorithms to Problems (3) and (6) as they come with slightly different
objective functions and constraint sets. Finally, in Section V, we deploy
these algorithms to illustrate and complement numerically our theoretical
findings (tightness of the expected asymptotic behavior, behavior when $S$ is
small).

##### Notations

For a matrix ${\mathbf{X}}\in{\mathbb{R}}^{M\times N}$, we denote by
${\mathbf{x}}_{n}$, $\underline{{\mathbf{x}}}_{m}$ and $[{\mathbf{X}}]_{mn}$
(or $x_{mn}$) its $n$-th column, $m$-th row, and $(m,n)$-th element,
respectively. Moreover, we use the notation $|{\mathbf{X}}|^{\circ 2}$ for its
element-wise modulus square. For a vector $\mathbf{y}\in{\mathbb{R}}^{T}$,
$\mathbf{y}[t]$ denotes its $t$-th element and $\mathbf{y}^{\intercal}$
denotes its transpose. We write $\operatorname*{Diag}(\mathbf{y})$ for the
diagonal matrix formed out of the vector $\mathbf{y}$, and
$\operatorname*{Diag}(\mathbf{X})$ for the vector formed by the diagonal
elements of $\mathbf{X}$. Finally, $X_{n}\overset{p}{\rightarrow}X$ stands for
the convergence in probability of the sequence $\\{X_{n}\\}_{n}$ of random
variables toward the random variable $X$.

##  II Preliminaries

###  II-A Definitions and Useful Results

First, from now on, and throughout the paper, we set $\epsilon_{0}>0$. Then,
let us provide a reformulation of the moment conditions in (1) that will be
relevant in our analysis.

######  Proposition 1.

The moment conditions (1) can be equivalently expressed as

$\displaystyle{\mathbb{E}}(\mathbf{y}_{n})$ $\displaystyle=\mathbf{0},$ (7)
$\displaystyle{\bar{\bm{\Phi}}}\bm{\Sigma}_{n}{\bar{\bm{\Phi}}}^{\intercal}$
$\displaystyle=\mathrm{Diag}(\bar{\mathbf{W}}\bar{\mathbf{h}}_{n}),$ (8)

where $\bm{\Sigma}_{n}={\mathbb{E}}(\mathbf{y}_{n}\mathbf{y}_{n}^{\intercal})$
denotes the covariance matrix of the $n$-th column of ${\mathbf{Y}}$.

###### Proof.

The mean condition in (1) implies that for any $n$,
${\mathbb{E}}({\bar{\bm{\Phi}}}\mathbf{y}_{n})=\mathbf{0}$, thus we have
${\mathbb{E}}(\mathbf{y}_{n})=\mathbf{0}$ due to the invertibility of
${\bar{\bm{\Phi}}}$. The variance condition implies that the diagonal terms of
${\bar{\bm{\Phi}}}\bm{\Sigma}_{n}{\bar{\bm{\Phi}}}^{\intercal}$ equal to the
vector $\bar{\mathbf{W}}\bar{\mathbf{h}}_{n}$. This is because, by definition,
$[{\bar{\bm{\Phi}}}\bm{\Sigma}_{n}{\bar{\bm{\Phi}}}^{\intercal}]_{mm}=[\bar{\mathbf{W}}\bar{\mathbf{H}}]_{mn}=(\bar{\mathbf{W}}\bar{\mathbf{h}}_{n})[m]$.
Similarly, the covariance condition implies that the off-diagonal terms of
${\bar{\bm{\Phi}}}\bm{\Sigma}_{n}{\bar{\bm{\Phi}}}^{\intercal}$ are zero. ∎

Finally, to complete the formulation of Problems (3) and (6) we define below
the constraint set $F_{K}$ as well as the empirical expectation operator
${\mathbb{E}}_{S}(\cdot)$ used in the definition of $C_{S}$ in (3).

##### Constraint Set

we define the constraint set $F_{K}$ as

$F_{K}=\\{({\mathbf{W}},{\mathbf{H}})\in{\mathbb{R}}^{M\times
K}_{+}\times{\mathbb{R}}^{K\times N}_{+}\;|\;\forall
k,\|\mathbf{w}_{k}\|_{1}=1\\},$ (9)

on which $C_{S}$ is well-defined (i.e., no singularity) when
${\epsilon_{0}}>0$. The normalization constraint on each column of
${\mathbf{W}}$ allows to break the scaling ambiguity between ${\mathbf{W}}$
and ${\mathbf{H}}$. It is worth mentioning that there are other types
ambiguities in NMF [16, 9]. These include permutations of the columns of
${\mathbf{W}}$ and rows of ${\mathbf{H}}$, or local nonnegativity-preserving
rotations. Yet, considering the scaling ambiguity turns out to be an important
ingredient to prove the existence of a solution to TL-NMF and JD+NMF (see the
proof of Theorem 1 in Appendix B). This motivates the definition of $F_{K}$ in
(9).

##### Empirical Expectation

We consider the empirical expectation of $[{\bm{\Phi}}{\mathbf{Y}}]_{mn}^{2}$,
that is

${\mathbb{E}}_{S}([{\bm{\Phi}}{\mathbf{Y}}]_{mn}^{2})=\frac{1}{S}\sum_{s=1}^{S}[{\bm{\Phi}}{\mathbf{Y}}^{(s)}]_{mn}^{2}.$
(10)

Moreover, we use the notation
${\mathbb{E}}_{S}(|{\bm{\Phi}}{\mathbf{Y}}|^{\circ 2})$ to refer to the point-
wise empirical expectation, i.e.,
$[{\mathbb{E}}_{S}(|{\bm{\Phi}}{\mathbf{Y}}|^{\circ
2})]_{mn}={\mathbb{E}}_{S}([{\bm{\Phi}}{\mathbf{Y}}]_{mn}^{2})$. Finally, we
introduce the following empirical covariance matrix:

$\bm{\Sigma}_{n,S}=\frac{1}{S}\sum_{s=1}^{S}{\mathbf{y}}_{n}^{(s)}({\mathbf{y}}_{n}^{(s)})^{\intercal}.$
(11)

###  II-B From TL-NMF to JD+NMF

For any integer $S>0$, one can easily verify the following decomposition of
the TL-NMF objective in (4),

$C_{S}({\bm{\Phi}},{\mathbf{W}},{\mathbf{H}})=L_{S}({\bm{\Phi}})+I_{S}({\bm{\Phi}},{\mathbf{W}},{\mathbf{H}}),$
(12)

with $L_{S}$ and $I_{S}$ defined by

$\displaystyle
L_{S}({\bm{\Phi}})=MN+\sum_{m,n=1}^{M,N}\log({\mathbb{E}}_{S}([{\bm{\Phi}}{\mathbf{Y}}]_{mn}^{2})+{\epsilon_{0}}),$
(13) $\displaystyle
I_{S}({\bm{\Phi}},{\mathbf{W}},{\mathbf{H}})=D_{\epsilon_{0}}({\mathbb{E}}_{S}(|{\bm{\Phi}}{\mathbf{Y}}|^{\circ
2})|{\mathbf{W}}{\mathbf{H}}),$ (14)

where $D_{\epsilon_{0}}(\cdot|\cdot)$ denotes a regularized form of the
Itakura-Saito (IS) divergence. For two non-negative matrices $\mathbf{A}$ and
$\mathbf{B}$, it is defined by

$D_{\epsilon_{0}}(\mathbf{A}|\mathbf{B})=\sum_{m,n}\left(\frac{[\mathbf{A}]_{mn}+{\epsilon_{0}}}{[\mathbf{B}]_{mn}+{\epsilon_{0}}}-\log\left(\frac{[\mathbf{A}]_{mn}+{\epsilon_{0}}}{[\mathbf{B}]_{mn}+{\epsilon_{0}}}\right)-1\right).$
(15)

The term $I_{S}$ in (14) is related to Itakura-Saito NMF (IS-NMF) [2]. Indeed,
when ${\bm{\Phi}}$ is fixed, the minimization of $I_{S}$ with respect to
$({\mathbf{W}},{\mathbf{H}})\in F_{K}$ produces a NMF of
${\mathbb{E}}_{S}(|{\bm{\Phi}}{\mathbf{Y}}|^{\circ 2})$ where the fit is
measured by the IS divergence.

The term $L_{S}$ is related to joint-diagonalization (JD), as specified by
Lemma 1.

######  Lemma 1.

For all ${\bm{\Phi}}\in\mathcal{O}(M)$,

$L_{S}({\bm{\Phi}})=MN+\sum_{n=1}^{N}\log\det\operatorname*{Diag}({\bm{\Phi}}(\bm{\Sigma}_{n,S}+{\epsilon_{0}}\mathbf{I}){\bm{\Phi}}^{\intercal}),$
(16)

which corresponds, up to a constant term, to the JD criterion derived in [14].

###### Proof.

See Appendix A. ∎

The minimization of $L_{S}$ over $\mathcal{O}(M)$ leads to an orthogonal
matrix ${\bm{\Phi}}$ such that
${\bm{\Phi}}\bm{\Sigma}_{n,S}{\bm{\Phi}}^{\intercal}$ is as diagonal as
possible. The additional diagonal matrix ${\epsilon_{0}}\mathbf{I}$ in (16)
can be interpreted (up to a constant normalization) as a linear shrinkage
estimator for covariance matrix estimation [17].

Hence, from the decomposition (12), TL-NMF can be interpreted as a trade-off
between the JD of the covariance matrices $\\{\bm{\Sigma}_{n,S}\\}_{n}$ and
the NMF of the variance matrix
${\mathbb{E}}_{S}(|{\bm{\Phi}}{\mathbf{Y}}|^{\circ 2})$. In contrast, the
alternative scheme (6) attempts to achieve the same goal in a two-step fashion
by first solving a JD problem to obtain ${\bm{\Phi}}^{\bullet}$ and then
solving an IS-NMF problem given ${\bm{\Phi}}^{\bullet}$ to get
$({\mathbf{W}}^{\bullet},{\mathbf{H}}^{\bullet})$.

######  Remark 1.

From Lemma 1, one can see that the “regularization” parameter
${\epsilon_{0}}>0$ ensures that $L_{S}$ is continuous over $\mathcal{O}(M)$.
Indeed, when $S<M$, $\bm{\Sigma}_{n,S}$ have vanishing eigenvalues meaning
that some diagonal entries of
${\bm{\Phi}}\bm{\Sigma}_{n,S}{\bm{\Phi}}^{\intercal}$ may be zero. Hence,
taking ${\epsilon_{0}}>0$ ensures that all the log-terms of (16) are non-
degenerate.

######  Remark 2.

The JD criterion derived in [14] and essentially given by (16) does not assume
${\bm{\Phi}}$ to be explicitly orthogonal but merely non-singular. Many
orthogonal JD algorithms were designed in the early age of ICA when whitening
was still a standard data pre-processing step. A seminal example is the Jacobi
algorithm by Cardoso and Souloumiac [18] for JD with a least-squares
criterion. However, to ensure optimal one-step performance [19], the whitening
step was eventually dropped in the ICA community. Non-orthogonal JD became the
mainstream, and many algorithms ensued for various criteria (such as based on
least-squares, maximum likelihood or information measures), see, e.g., [20,
21, 22, 23, 24]. More recently, [25] leverages Riemannian optimization to
unify many existing methods and introduce new ones under various constraints
for ${\bm{\Phi}}$.

In this work, we restrict ${\bm{\Phi}}$ to be orthogonal because TL-NMF aims
to generalize short-time frequency transforms for which orthogonality is
somehow natural or desired. Moreover when $S=\infty$, we showed that the row
subspaces of $\bar{{\bm{\Phi}}}$ are identifiable [5]. Such an analysis is
more difficult in the general case (see [26, 27] for related discussions).
Still, removing the orthogonality assumption of ${\bm{\Phi}}$ could be
beneficial in practice and forms a relevant research direction (see further
discussion in Section VI).

##  III Relationship between JD+NMF and TL-NMF

###  III-A Existence of a Solution

Before analyzing the relation between TL-NMF and JD+NMF, it is worth checking
that they both admit at least one solution. To the best of our knowledge, this
question has never been addressed in the TL-NMF literature, nor in NMF
literature (i.e., existence of global minimizer(s) for $I_{S}$, see Lemma 3).

######  Theorem 1.

The solution sets $\Omega^{*}$ of TL-NMF (Problem (3)) and $\Omega^{\bullet}$
of JD+NMF (Problem (6)) are nonempty and compact.

###### Proof.

See Appendix B. ∎

######  Remark 3.

The assumption ${\epsilon_{0}}>0$ is necessary for the functions $L_{S}$ and
$I_{S}$ to be well-defined and to show the existence of solutions to TL-NMF
and JD+NMF. We chose to include ${\epsilon_{0}}$ in the definitions of $L_{S}$
and $I_{S}$ but we could alternatively have set ${\epsilon_{0}}=0$ in $L_{S}$
and $I_{S}$, and add an inequality of the form
${\mathbf{W}}{\mathbf{H}}\geq{\epsilon_{0}}$ in the constraint set $F_{K}$.
This is related to the approach followed by [28] to study the convergence of
NMF with a large range of divergences (including the generalized Kullback-
Leibler divergence [29]) under constraints of the form
${\mathbf{W}}\geq{\epsilon_{0}}$, ${\mathbf{H}}\geq{\epsilon_{0}}$.

###  III-B When does JD+NMF meet TL-NMF?

In Theorem 2, we characterize the closeness between the solutions of JD+NMF
and the solutions of TL-NMF. A graphical illustration of this result is
depicted in Figure 1.

######  Theorem 2.

Let $\Omega^{*}$ and $\Omega^{\bullet}$ be defined as in Theorem 1. Define

$\displaystyle\underline{\lambda}^{*}=\min_{({\bm{\Phi}},{\mathbf{W}},{\mathbf{H}})\in\Omega^{*}}I_{S}({\bm{\Phi}},{\mathbf{W}},{\mathbf{H}}),$
(17)
$\displaystyle\bar{\lambda}^{*}=\max_{({\bm{\Phi}},{\mathbf{W}},{\mathbf{H}})\in\Omega^{*}}I_{S}({\bm{\Phi}},{\mathbf{W}},{\mathbf{H}}),$
(18)

and similarly $\underline{\lambda}^{\bullet}$ and $\bar{\lambda}^{\bullet}$ by
replacing $\Omega^{*}$ by $\Omega^{\bullet}$. Then,

$0\leq\underline{\lambda}^{*}\leq\bar{\lambda}^{*}\leq\underline{\lambda}^{\bullet}\leq\bar{\lambda}^{\bullet}<+\infty.$
(19)

Moreover,

$\displaystyle\Omega^{\bullet}\cap\Omega^{*}\neq\emptyset\;\Longleftrightarrow\;\bar{\lambda}^{*}=\underline{\lambda}^{\bullet},$
(20)
$\displaystyle\Omega^{\bullet}=\Omega^{*}\;\Longleftrightarrow\;\underline{\lambda}^{*}=\bar{\lambda}^{\bullet}.$
(21)

###### Proof.

See Appendix C. ∎

Figure 1:  Graphical illustration of Theorem 2. The relationship between the
solutions sets $\Omega^{*}$ and $\Omega^{\bullet}$ of TL-NMF and JD+NMF is
described in terms of specific level lines of $I_{S}$. Note that we
represented here $\Omega^{*}$ and $\Omega^{\bullet}$ as countable sets for the
sake of illustration, but they may not be.

There are at least two situations where JD+NMF is equivalent to TL-NMF (i.e.,
(21) is satisfied):

  * •

$S=\infty$ and $K\geq\bar{K}$,

  * •

$S<\infty$ and $K\geq\min\\{M,N\\}$.

Indeed, in these two cases we have that,
$\forall{\bm{\Phi}}\in\mathcal{O}(M)$,

$\exists({\mathbf{W}},{\mathbf{H}})\in F_{K},\text{ such that
}\;{\mathbb{E}}_{S}(|{\bm{\Phi}}{\mathbf{Y}}|^{\circ
2})={\mathbf{W}}{\mathbf{H}},$ (22)

i.e., there exists an exact factorization of
${\mathbb{E}}_{S}(|{\bm{\Phi}}{\mathbf{Y}}|^{\circ 2})$. When
$K\geq\min\\{M,N\\}=N$ it suffices to populate ${\mathbf{W}}$ with
${\mathbb{E}}_{S}(|{\bm{\Phi}}{\mathbf{Y}}|^{\circ 2})$, ${\mathbf{H}}$ with
${\mathbf{I}}_{N}$, and to complete with zeros. A similar exact factorization
can be obtained when $K\geq\min\\{M,N\\}=M$ by switching the roles of
${\mathbf{W}}$ and ${\mathbf{H}}$. For the case $S=\infty$, we refer the
reader to the proof of Proposition 2 in Appendix D. Consequently, in such
cases, we get from (22) that

$\forall{\bm{\Phi}}\in\mathcal{O}(M),\;\min_{({\mathbf{W}},{\mathbf{H}})\in
F_{K}}I_{S}({\bm{\Phi}},{\mathbf{W}},{\mathbf{H}})=0.$ (23)

It then follows from (3) (resp., (6)) and (12) that

$\forall({\bm{\Phi}},{\mathbf{W}},{\mathbf{H}})\in\Omega^{*}\,(\mathrm{resp.,}\,\Omega^{\bullet}),\;I_{S}({\bm{\Phi}},{\mathbf{W}},{\mathbf{H}})=0$
(24)

and thus $\underline{\lambda}^{*}=\bar{\lambda}^{\bullet}=0$.

The situation is more involved when $S<\infty$ and $K<\min\\{M,N\\}$. In this
case, we get from Theorem 2 that the closeness between the solutions of JD+NMF
and TL-NMF is controlled by the inner
$\underline{\lambda}^{\bullet}-\bar{\lambda}^{*}$ and the outer
$\bar{\lambda}^{\bullet}-\underline{\lambda}^{*}$ gaps. Indeed, by definition,
we get that for any
$({{\bm{\Phi}}}^{*},{{\mathbf{W}}}^{*},{{\mathbf{H}}}^{*})\in\Omega^{*}$ and
$({\bm{\Phi}}^{\bullet},{\mathbf{W}}^{\bullet},{\mathbf{H}}^{\bullet})\in\Omega^{\bullet}$,

$0\leq\underline{\lambda}^{\bullet}-\bar{\lambda}^{*}\leq
I_{S}^{\bullet}-I_{S}^{*}\leq\bar{\lambda}^{\bullet}-\underline{\lambda}^{*},$
(25)

where $\displaystyle
I_{S}^{\bullet}=I_{S}({\bm{\Phi}}^{\bullet},{\mathbf{W}}^{\bullet},{\mathbf{H}}^{\bullet})$
and $\displaystyle
I_{S}^{*}={I_{S}}({{\bm{\Phi}}}^{*},{{\mathbf{W}}}^{*},{{\mathbf{H}}}^{*})$.
Moreover, one can see that the outer gap also controls the proximity between
minimizers in terms of $C_{S}$ values because

$0\leq C_{S}^{\bullet}-C_{S}^{*}\leq
I_{S}^{\bullet}-I_{S}^{*}\leq\bar{\lambda}^{\bullet}-\underline{\lambda}^{*},$
(26)

where the first inequality is obtained by combining (12) with the fact that
$\displaystyle L_{S}^{\bullet}\leq L_{S}^{*}$. Note that here, $\displaystyle
C_{S}^{\bullet}$ and $\displaystyle C_{S}^{*}$ (resp., $\displaystyle
L_{S}^{\bullet}$, and $\displaystyle L_{S}^{*}$) are the counterparts of
$\displaystyle I_{S}^{\bullet}$ and $\displaystyle I_{S}^{*}$ for the
objectives $C_{S}$ (resp., $L_{S}$). In Section III-C, we shall analyse the
asymptotic behaviour of the inner and outer gaps.

Finally, let us emphasize that the inequalities (19) also reveal that the
transforms ${\bm{\Phi}}^{\bullet}$ obtained by JD+NMF are at best as amenable
to NMF (in terms of IS-divergence) as any ${\bm{\Phi}}^{*}$ learned by TL-NMF.

######  Remark 4.

The above discussion about the quality of the learned ${\bm{\Phi}}$ in terms
of $I_{S}$ values may raise the following question: why not directly
minimizing $I_{S}$ with respect to $({\bm{\Phi}},{\mathbf{W}},{\mathbf{H}})$?
The latter approach was in particular the one considered in the original TL-
NMF paper [4] (with an additional sparsity-enforcing term for ${\mathbf{H}}$),
an arbitrary choice inherited from IS-NMF [30]. As a matter of fact, in light
of our results, several arguments play in favor of minimizing $C_{S}$ rather
than $I_{S}$. First, under the GCM, $C_{S}$ comes with a probabilistic
interpretation (log-likelihood function). Second, when $S=\infty$, there exist
$({\mathbf{W}},{\mathbf{H}})$ such that
$I_{S}({\bm{\Phi}},{\mathbf{W}},{\mathbf{H}})=0$ for all ${\bm{\Phi}}$ (see
(23)). This means that $I_{S}$ does not constitutes a good measure to
discriminate the ${\bm{\Phi}}$. In contrast, the minimization of $C_{S}$
allows to identify the true transform ${\bar{\bm{\Phi}}}$ (the term
$L_{S}({\bm{\Phi}})$ acts somehow as a penalization term) [5]. Lastly, this
ability to learn meaningful transform (e.g., close to ${\bar{\bm{\Phi}}}$) by
minimizing $C_{S}$ appears to be also true when $S$ is small (see Section
V-D).

###  III-C Asymptotic Analysis

In this section we analyze the closeness between JD+NMF and TL-NMF when
$S<\infty$ and $K\in[\bar{K},\min\\{M,N\\})$.111The more intricate case
$K<\bar{K}$ is not considered in the paper. We derive a sufficient condition
under which the outer gap $\bar{\lambda}^{\bullet}-\underline{\lambda}^{*}$
(and consequently the inner gap
$\underline{\lambda}^{\bullet}-\bar{\lambda}^{*}$) converges (in probability)
to zero as $S$ grows. Given that

$\underline{\lambda}^{\bullet}-\bar{\lambda}^{*}\leq\bar{\lambda}^{\bullet}-\underline{\lambda}^{*}\leq\bar{\lambda}^{\bullet},$
(27)

we first provide in Proposition 2 an upper bound of $\bar{\lambda}^{\bullet}$
that is independent of the solution sets $\Omega^{*}$ and $\Omega^{\bullet}$.

######  Proposition 2.

Under condition (1) and $K\geq\bar{K}$, we have for all $S>0$,

$\bar{\lambda}^{\bullet}\leq\max_{{\bm{\Phi}}\in\mathcal{O}(M)}\,Q_{S}({\bm{\Phi}})$
(28)

where

$Q_{S}({\bm{\Phi}})=D_{\epsilon_{0}}({\mathbb{E}}_{S}(|{\bm{\Phi}}{\mathbf{Y}}|^{\circ
2})\,|\,{\mathbb{E}}(|{\bm{\Phi}}{\mathbf{Y}}|^{\circ 2})).$ (29)

###### Proof.

See Appendix D. ∎

We now derive in Theorem 3 a sufficient condition for the upper bound
$\max_{{\bm{\Phi}}}Q_{S}({\bm{\Phi}})$ to converge to zero in probability.

######  Theorem 3.

Assume that the empirical estimator
${\mathbb{E}}_{S}(|{\bm{\Phi}}{\mathbf{Y}}|^{\circ 2})$ converges uniformly
toward ${\mathbb{E}}(|{\bm{\Phi}}{\mathbf{Y}}|^{\circ 2})$ in probability,
i.e.,

$\max_{(m,n)}\max_{{\bm{\Phi}}\in\mathcal{O}(M)}\bigg{|}{\mathbb{E}}_{S}([{\bm{\Phi}}{\mathbf{Y}}]_{mn}^{2})-{\mathbb{E}}([{\bm{\Phi}}{\mathbf{Y}}]^{2}_{mn})\bigg{|}\overset{p}{\to}0,$
(30)

Then

$\max_{{\bm{\Phi}}\in\mathcal{O}(M)}Q_{S}({\bm{\Phi}})\overset{p}{\to}0\;\text{
as }\;S\to\infty.$ (31)

###### Proof.

See Appendix E. ∎

Under the GCM, we can go one step further than Theorem 3 by deriving the
convergence rate of $\max_{{\bm{\Phi}}}Q_{S}({\bm{\Phi}})$.

######  Theorem 4.

Under GCM defined by (2), condition (30) is always satisfied. Moreover, for
$t>0$ and $S$ large enough such that $h_{S}=3\frac{\sqrt{M}+t}{\sqrt{S}}<1$,
we have

$\max_{{\bm{\Phi}}\in\mathcal{O}(M)}Q_{S}({\bm{\Phi}})<MN\frac{h_{S}^{2}}{1-h_{S}},$
(32)

with probability at least $(1-2e^{-t^{2}/2})^{N}$.

###### Proof.

See Appendix E. ∎

It is worth mentioning that the proof of Theorem 4 makes use of existing
results on covariance matrix estimation [31] which can be generalized to the
case where the entries of ${\mathbf{Y}}$ follow some sub-Gaussian
distributions [32].

######  Remark 5.

From Theorem 4, we get that with high probability, as $S$ grows,
$\max_{\bm{\Phi}}Q_{S}({\bm{\Phi}})$ converges at least at a rate of $O(1/S)$.
Indeed, for any $t>0$, there exists $S^{*}$ such that for all $S\geq S^{*}$,
$h_{S}<1$. This means that for $S$ sufficiently large, (32) holds with
probability almost one. Then, the claim comes by observing that, for $S$
large, $h_{S}=O(1/\sqrt{S})$.

######  Remark 6.

We can provide insights on the tightness of the uniform bound (32) by
analyzing a point-wise bound such as the asymptotic decay rate of
$Q_{S}({\bar{\bm{\Phi}}})$ with $S$. Let us first observe that

$\mkern-15.0muQ_{S}({\bar{\bm{\Phi}}})=\mkern-10.0mu\sum_{m,n=1}^{M,N}\mkern-10.0muf\left(X_{mn}^{\bar{\bm{\Phi}}}\right)\text{
for
}X_{mn}^{\bar{\bm{\Phi}}}=\frac{{\mathbb{E}}_{S}([{\bar{\bm{\Phi}}}{\mathbf{Y}}]^{2}_{mn})+{\epsilon_{0}}}{{\mathbb{E}}([{\bar{\bm{\Phi}}}{\mathbf{Y}}]^{2}_{mn})+{\epsilon_{0}}}$
(33)

where $f(x)=x-\log(x)-1$. Moreover, from (1), we get that
${\mathbb{E}}([{\bar{\bm{\Phi}}}{\mathbf{Y}}]^{2}_{mn})>0$ and, for $S$
sufficiently large, we can assume that
${\mathbb{E}}_{S}([{\bar{\bm{\Phi}}}{\mathbf{Y}}]^{2}_{mn})>0$ (see discussion
in Remark 1). Hence, we can set (for the purpose of this remark)
${\epsilon_{0}}=0$. We then get from the definition of ${\mathbb{E}}_{S}$ in
(10) that $SX_{mn}^{\bar{\bm{\Phi}}}$ follows a Chi-squared distribution
$\mathcal{X}_{S}$ of degree $S$. We deduce that

$\displaystyle{\mathbb{E}}(f\left(X_{mn}^{\bm{\Phi}}\right))$
$\displaystyle=\frac{1}{S}{\mathbb{E}}(SX_{mn}^{\bm{\Phi}})+\log(S)-{\mathbb{E}}\left(\log\left(SX_{mn}^{\bm{\Phi}}\right)\right)-1$
$\displaystyle=\log(S)-\psi(S/2)-\log(2)$ (34)

using the facts that ${\mathbb{E}}(SX_{mn}^{\bm{\Phi}})=S$ and
${\mathbb{E}}(\log(SX_{mn}^{\bm{\Phi}}))=\psi(S/2)+\log(2)$ where $\psi$
denotes the di-gamma function. Then, from the weak law of large number, we
obtain that

$\frac{1}{NM}Q_{S}({\bar{\bm{\Phi}}})\overset{p}{\underset{NM\to\infty}{\to}}(\log(S)-\psi(S/2)-\log(2)).$
(35)

Noticing that, when $S$ is large, $\psi(S/2)\approx\log(S/2)-1/S$, we conclude
(for $M$, $N$ large and ${\epsilon_{0}}=0$) that the decay rate of
$Q_{S}({\bar{\bm{\Phi}}})$ is of the order of $O(1/S)$ which matches the
uniform bound (32).

##  IV Optimization Methods

Algorithm 1 TL-NMF Solver

In: $(J,J_{\mathrm{NMF}},J_{\mathrm{TL}})\in\mathbb{N}^{3}$,
$({\bm{\Phi}}_{0},{\mathbf{W}}_{0},{\mathbf{H}}_{0})\in\mathcal{O}(M)\times
F_{K}$

1:for $j=0:J-1$ do

2:
$({\mathbf{W}}_{j+1},{\mathbf{H}}_{j+1})\leftarrow\textsf{MU}\left(C_{S}({\bm{\Phi}}_{j},\cdot,\cdot);J_{\mathrm{NMF}},{\mathbf{W}}_{j},{\mathbf{H}}_{j}\right)$

3:
${\bm{\Phi}}_{j+1}\leftarrow\textsf{QN}(C_{S}(\cdot,{\mathbf{W}}_{j+1},{\mathbf{H}}_{j+1});J_{\mathrm{TL}},{\bm{\Phi}}_{j})$

4:end for

Out: ${\bm{\Phi}}_{J},{\mathbf{W}}_{J},{\mathbf{H}}_{J}$

Algorithm 2 JD+NMF Solver

In: $(J,J_{\mathrm{NMF}},J_{\mathrm{TL}})\in\mathbb{N}^{3}$,
$({\bm{\Phi}}_{0},{\mathbf{W}}_{0},{\mathbf{H}}_{0})\in\mathcal{O}(M)\times
F_{K}$

1:${\bm{\Phi}}_{J}\leftarrow\textsf{QN}(L_{S};J\times
J_{\mathrm{TL}},{\bm{\Phi}}_{0})$

2:$({\mathbf{W}}_{J},{\mathbf{H}}_{J})\leftarrow\textsf{MU}(I_{S}({\bm{\Phi}}_{J},\cdot,\cdot);J\times
J_{\mathrm{NMF}},{\mathbf{W}}_{0},{\mathbf{H}}_{0})$

Out: ${\bm{\Phi}}_{J},{\mathbf{W}}_{J},{\mathbf{H}}_{J}$ Note that we
introduced the redundant parameter $J$ in Algorithm 2 so that the input
parameters are exactly the same as for the TL-NMF solver in Algorithm 1. This
allows us to ease the computational comparison of the two solvers (see Section
V-E). The sub-procedures of the form $\textsf{AA}(g;J,\mathbf{X})$ should be
read as: executing the AA method on the function $g$ from the initial point
$\mathbf{X}$ for $J$ iterations.

To numerically address TL-NMF and JD+NMF problems, we deploy the alternating
optimization methods outlined in Algorithms 1 and 2, respectively. They both
rely on the standard multiplicative updates (MU) for the NMF factors
${\mathbf{W}},{\mathbf{H}}$, and a quasi-Newton (QN) method for the update of
the transform ${\bm{\Phi}}$. Observing from (12) that, given ${\bm{\Phi}}$,
the minimization of $C_{S}({\bm{\Phi}},\cdot,\cdot)$ boils down to the
minimization of $I_{S}({\bm{\Phi}},\cdot,\cdot)$, the MU steps at line 2 of
Algorithms 1 and 2 are the same and are implemented according to [30], as
summarized in Algorithm 3. The latter is a block-descent algorithm in which
the factors are updated in turn. The update of each factor is carried out by
one step of majorization-equalization, a variant of majorization-minimization
[33] that produces an acceleration while ensuring nonincreasingness of the
objective function [30]. This results in the celebrated MU algorithm that is
customary in NMF [9].

Concerning the update of ${\bm{\Phi}}$, we adapt the quasi-Newton (QN)
algorithms proposed in [34, 35]. These methods were initially developed to
solve TL-NMF and JD problems with slightly different objective functions and
constraints than those considered in the present work. Numerically, they have
been shown to improve the convergence rate of gradient-descent or Jacobi-based
methods [18, 36, 20].

Algorithm 3 Multiplicative updates to minimize
$I_{S}({\bm{\Phi}},\cdot,\cdot)$ [30]

In: $J_{\mathrm{NMF}}\in\mathbb{N}$, $({\mathbf{W}}_{0},{\mathbf{H}}_{0})\in
F_{K}$, ${\mathbf{V}}={\mathbb{E}}_{S}(|{\bm{\Phi}}{\mathbf{Y}}|^{\circ 2})$

1:for $j=0:J_{\mathrm{NMF}}-1$ do

2:
$\hat{{\mathbf{V}}}\leftarrow{\mathbf{W}}_{j}{\mathbf{H}}_{j}+{\epsilon_{0}}\mathbf{1}_{M\times
N}$

3:
$\displaystyle{\mathbf{H}}_{j+1}\leftarrow{\mathbf{H}}_{j}\circ\left[\frac{{\mathbf{W}}^{\intercal}_{j}((\hat{{\mathbf{V}}})^{\circ-2}\circ{\mathbf{V}})}{{\mathbf{W}}^{\intercal}_{j}(\hat{{\mathbf{V}}})^{\circ-1}}\right]$

4:
$\hat{{\mathbf{V}}}\leftarrow{\mathbf{W}}_{j}{\mathbf{H}}_{j+1}+{\epsilon_{0}}\mathbf{1}_{M\times
N}$

5:
$\displaystyle{\mathbf{W}}_{j+1}\leftarrow{\mathbf{W}}_{j}\circ\left[\frac{((\hat{{\mathbf{V}}})^{\circ-2}\circ{\mathbf{V}}){\mathbf{H}}^{\intercal}_{j+1}}{(\hat{{\mathbf{V}}})^{\circ-1}{\mathbf{H}}^{\intercal}_{j+1}}\right]$

6: Normalize ${\mathbf{W}}_{j+1}$, ${\mathbf{H}}_{j+1}$ to remove scale
ambiguity.

7:end for

Out: ${\mathbf{W}}_{J_{\mathrm{NMF}}},{\mathbf{H}}_{J_{\mathrm{NMF}}}$

###  IV-A General Principle of QN Methods over $\mathcal{O}(M)$

For both problems, the main difficulty in deriving a QN method for updating
${\bm{\Phi}}$ comes from the handling of the orthogonality constraint (i.e.,
${\bm{\Phi}}\in\mathcal{O}(M)$). In this work, we follow the standard approach
that consists in defining a local parameterization
$\rho_{{\bm{\Phi}}_{j}}:\Upsilon\subset{\mathbb{R}}^{M\times
M}\rightarrow\mathcal{O}(M)$ of the neighborhood of an iterate
${\bm{\Phi}}_{j}$ (see for instance [37]). The idea is then to compute a QN
direction based on a quadratic approximation of the local objective function
$g\circ\rho_{{\bm{\Phi}}_{j}}$ (here $g$ stands for either
$C_{S}(\cdot,{\mathbf{W}},{\mathbf{H}})$ or $L_{S}$). From the second-order
Taylor expansion around $\mathbf{0}$, we get

$g(\rho_{{\bm{\Phi}}_{j}}({\mathbf{E}}))=g({\bm{\Phi}}_{j})+{\langle}{\mathbf{G}},{\mathbf{E}}{\rangle}+\frac{1}{2}{\langle}{\mathbf{E}}|{\mathcal{H}}|{\mathbf{E}}{\rangle}+O\left(\|{\mathbf{E}}\|^{3}\right),$
(36)

using the fact that $\rho_{{\bm{\Phi}}_{j}}(\mathbf{0})={\bm{\Phi}}_{j}$, and
denoting respectively by ${\mathbf{G}}\in{\mathbb{R}}^{M\times M}$ and
${\mathcal{H}}\in{\mathbb{R}}^{M\times M\times M\times M}$ the gradient matrix
and the Hessian tensor of $g\circ\rho_{{\bm{\Phi}}_{j}}$ at $\mathbf{0}$. In
(36), the inner products should be read as
${\langle}{\mathbf{G}},{\mathbf{E}}{\rangle}=\sum_{a,b}[{\mathbf{G}}]_{ab}[{\mathbf{E}}]_{ab}$
and
${\langle}{\mathbf{E}}|{\mathcal{H}}|{\mathbf{E}}{\rangle}=\sum_{a,b,c,d}[{\mathcal{H}}]_{abcd}[{\mathbf{E}}]_{ab}[{\mathbf{E}}]_{cd}$.

Because the Hessian is usually costly to compute, we will consider an
approximation of it, denoted $\tilde{\mathcal{H}}$, and define a QN direction
through the resolution of

${\mathbf{E}}_{j}=\operatorname*{arg\,min}_{{\mathbf{E}}\in\Upsilon}{\langle}{\mathbf{G}},{\mathbf{E}}{\rangle}+\frac{1}{2}{\langle}{\mathbf{E}}|\tilde{{\mathcal{H}}}|{\mathbf{E}}{\rangle}.$
(37)

Given this QN direction, the current estimate ${\bm{\Phi}}_{j}$ is updated
according to

${\bm{\Phi}}_{j+1}=\rho_{{\bm{\Phi}}_{j}}(\eta{\mathbf{E}}_{j}),$ (38)

where $\eta>0$ is obtained by line search (e.g., backtracking [34] or Wolfe
[38]). This generic QN approach is summarized in Algorithm 4. In the next two
sections, we provide details on its instantiation to our problems. More
precisely, the main task is to compute the gradient ${\mathbf{G}}$, the
Hessian approximation $\tilde{{\mathcal{H}}}$, as well as the QN direction in
(37).

Algorithm 4 Generic quasi-Newton method to minimize a function $g$ over
$\mathcal{O}(M)$

In: $J_{\mathrm{QN}}\in\mathbb{N}$, ${\bm{\Phi}}_{0}\in\mathcal{O}(M)$,
$\rho_{{\bm{\Phi}}}(\cdot)$

1:for $j=0:J_{\mathrm{QN}}-1$ do

2: Compute the gradient ${\mathbf{G}}$ and the Hessian approximation
$\tilde{\mathcal{H}}$ of the local function $g\circ\rho_{{\bm{\Phi}}_{j}}$

3:
${\mathbf{E}}_{j}\leftarrow\operatorname*{arg\,min}_{{\mathbf{E}}\in\Upsilon}{\langle}{\mathbf{G}},{\mathbf{E}}{\rangle}+\frac{1}{2}{\langle}{\mathbf{E}}|\tilde{{\mathcal{H}}}|{\mathbf{E}}{\rangle}$

4: Compute $\eta>0$ via line search

5: ${\bm{\Phi}}_{j+1}\leftarrow\rho_{{\bm{\Phi}}_{j}}(\eta{\mathbf{E}}_{j})$

6:end for

Out: ${\bm{\Phi}}_{J_{\mathrm{TL}}}$

######  Remark 7.

In the following, we adopt two different strategies to parametrize the
neighborhood of an iterate ${\bm{\Phi}}_{j}$. The reason is that each one has
been motivated independently by previous works on TL-NMF [35] and JD [34].
Yet, let us emphasize that for both problems, a QN method could be derived
with each of the two local parameterizations described below.

###  IV-B QN Method for TL-NMF

As in [35], we consider the local parameterization

$\rho_{\bm{\Phi}}({\mathbf{E}})=\exp({\mathbf{E}}){\bm{\Phi}},\;\forall{\mathbf{E}}\in\Upsilon,$
(39)

where $\Upsilon:=\left\\{{\mathbf{E}}\in{\mathbb{R}}^{M\times
M}:{\mathbf{E}}^{\intercal}=-{\mathbf{E}}\right\\}$ is the set of anti-
symmetric matrices (leading to $\exp({\mathbf{E}})\in\mathcal{O}(M)$). Then,
denoting ${\mathbf{X}}^{(s)}={\bm{\Phi}}{\mathbf{Y}}^{(s)}$, we follow the
same steps as in [35] to obtain the gradient of
$C_{S}(\cdot,{\mathbf{W}},{\mathbf{H}})\circ\rho_{{\bm{\Phi}}}$ at
$\mathbf{0}$,

$[{\mathbf{G}}]_{ab}=\frac{2}{S}\sum_{n,s=1}^{N,S}\frac{[{\mathbf{X}}^{(s)}]_{an}[{\mathbf{X}}^{(s)}]_{bn}}{[{\mathbf{W}}{\mathbf{H}}]_{an}+{\epsilon_{0}}},$
(40)

as well as the following diagonal Hessian approximation

$[\tilde{{\mathcal{H}}}]_{abcd}=\delta_{ac}\delta_{bd}[\bm{\Gamma}]_{ab}$ (41)

where

$[\bm{\Gamma}]_{ab}=\frac{2}{S}\sum_{n,s=1}^{N,S}\frac{[{\mathbf{X}}^{(s)}]_{bn}^{2}}{[{\mathbf{W}}{\mathbf{H}}]_{an}+{\epsilon_{0}}}.$
(42)

In Proposition 3, we provide the closed-form expression of the solution of
(37) for such a diagonal $\tilde{{\mathcal{H}}}$. Again, this result is an
extension of [35] to our setting (i.e., minimization of
$C_{S}(\cdot,{\mathbf{W}},{\mathbf{H}})$ instead of
$I_{S}(\cdot,{\mathbf{W}},{\mathbf{H}})$ in [35]).

######  Proposition 3.

For $\tilde{{\mathcal{H}}}$ in (41), a solution of Problem (37) is given by
${\mathbf{E}}\in\Upsilon$ such that

$[{\mathbf{E}}]_{ab}=\left\\{\begin{array}[]{ll}\displaystyle-\frac{[{\mathbf{G}}^{(\mathrm{anti})}]_{ab}}{[\bm{\Gamma}^{(\mathrm{sym})}]_{ab}}&\text{
if }[\bm{\Gamma}^{(\mathrm{sym})}]_{ab}\neq 0\\\ 0&\text{ if
}[\bm{\Gamma}^{(\mathrm{sym})}]_{ab}=0\end{array}\right.$ (43)

where
${\mathbf{G}}^{(\mathrm{anti})}=({\mathbf{G}}-{\mathbf{G}}^{\intercal})/2$ and
$\bm{\Gamma}^{(\mathrm{sym})}=(\bm{\Gamma}+\bm{\Gamma}^{\intercal})/2$.
Moreover, it is a descent direction, i.e.,
${\mathbf{G}}^{(\mathrm{anti})}\neq\mathbf{0}$ implies that
${\langle}{\mathbf{G}},{\mathbf{E}}{\rangle}<0$.

###### Proof.

See Appendix F-B. ∎

###  IV-C QN Method for JD

For the minimization of $L_{S}$, we rely on the local parameterization
discussed in [37], that is

$\rho_{\bm{\Phi}}({\mathbf{E}})=\pi({\bm{\Phi}}+{\mathbf{E}}{\bm{\Phi}}),\;\forall{\mathbf{E}}\in\Upsilon,$
(44)

where $\pi$ stands for the projection operator on $\mathcal{O}(M)$ and
$\Upsilon$ is again the set of anti-symmetric matrices. The gradient of
$L_{S}\circ\rho_{{\bm{\Phi}}}$ at $\mathbf{0}$ is then given by

$[{\mathbf{G}}]_{ab}=\frac{1}{N}\sum_{n=1}^{N}\left(\frac{[{\bm{\Phi}}(\bm{\Sigma}_{n,S}+{\epsilon_{0}}{\mathbf{I}}){\bm{\Phi}}^{\intercal}]_{ab}}{[{\bm{\Phi}}(\bm{\Sigma}_{n,S}+{\epsilon_{0}}{\mathbf{I}}){\bm{\Phi}}^{\intercal}]_{aa}}-\delta_{ab}\right).$
(45)

Similarly to [34] (where $L_{S}$ is optimized over the set of invertible
matrices), we consider the Hessian approximation

$[\tilde{{\mathcal{H}}}]_{abcd}=\delta_{ac}\delta_{bd}[{\bm{\Gamma}}]_{ab}+\delta_{ad}\delta_{bc}-2\delta_{abcd}$
(46)

where

$[{\bm{\Gamma}}]_{ab}=\frac{1}{N}\sum_{n=1}^{N}\frac{[{\bm{\Phi}}(\bm{\Sigma}_{n,S}+{\epsilon_{0}}{\mathbf{I}}){\bm{\Phi}}^{\intercal}]_{bb}}{[{\bm{\Phi}}(\bm{\Sigma}_{n,S}+{\epsilon_{0}}{\mathbf{I}}){\bm{\Phi}}^{\intercal}]_{aa}}.$
(47)

It is noteworthy to mention that each element of ${\mathbf{G}}$ and
$\bm{\Gamma}$ is well-defined as $\forall a\leq
M,{\bm{\Phi}}\in\mathcal{O}(M)$,
$[{\bm{\Phi}}(\bm{\Sigma}_{n,S}+{\epsilon_{0}}{\mathbf{I}}){\bm{\Phi}}^{\intercal}]_{aa}\geq{\epsilon_{0}}>0$.

Because the Hessian approximation in (46) is not diagonal, we cannot use
Proposition 3 to get the QN direction. Yet, the specific structure of
$\tilde{{\mathcal{H}}}$ in (46) allows to express the solution of (37) in
closed form, as stated in Proposition 4.

######  Proposition 4.

For $\tilde{{\mathcal{H}}}$ in (46), a solution of Problem (37) is given by
${\mathbf{E}}\in\Upsilon$ such that

$[{\mathbf{E}}]_{ab}=\left\\{\begin{array}[]{ll}\displaystyle-\frac{[{\mathbf{G}}^{(\mathrm{anti})}]_{ab}}{[{\bm{\Gamma}}^{(\mathrm{sym})}]_{ab}-1}&\text{
if }[{\bm{\Gamma}}^{{(\mathrm{sym})}}]_{ab}\neq 1\\\ 0&\text{ if
}[{\bm{\Gamma}}^{{(\mathrm{sym})}}]_{ab}=1\end{array}\right.$ (48)

where
${\mathbf{G}}^{(\mathrm{anti})}=({\mathbf{G}}-{\mathbf{G}}^{\intercal})/2$ and
${\bm{\Gamma}}^{(\mathrm{sym})}=({\bm{\Gamma}}+{\bm{\Gamma}}^{\intercal})/2$.
Moreover, it is a descent direction, i.e.,
${\mathbf{G}}^{(\mathrm{anti})}\neq\mathbf{0}$ implies that
${\langle}{\mathbf{G}},{\mathbf{E}}{\rangle}<0$.

###### Proof.

See Appendix F-C. ∎

As opposed to [34], here ${\mathbf{E}}$ is constrained to be anti-symmetric.
This is to ensure that, for ${\bm{\Phi}}\in\mathcal{O}(M)$,
${\mathbf{E}}{\bm{\Phi}}$ belongs to the tangent space of $\mathcal{O}(M)$ at
${\bm{\Phi}}$. This makes $\tilde{{\mathcal{H}}}$ invertible, discarding the
need to rely on a pseudo-inverse.

###  IV-D Computational Complexity

We now briefly discuss the computational complexity of the two solvers in
Algorithm 1 and 2. First of all, one can see from Algorithm 3 that an
iteration of MU has a complexity of the order of $O(KMN)$ (the main cost being
the products of matrices with sizes $M\times K$ and $K\times N$).

Concerning the QN step, the main cost comes from the computation of the
gradient and the Hessian approximation which lead to a complexity per
iteration of the order of $O(SNM^{2})$ for TL-NMF (see (40) and (41)) and
$O(NM^{3})$ for JD (see (45) and (46)).

The overall complexities of the two solvers are thus of the order of

  * •

$O\left(J(J_{\mathrm{NMF}}KMN+J_{\mathrm{TL}}SNM^{2})\right)$ for TL-NMF,

  * •

$O\left(J(J_{\mathrm{NMF}}KMN+J_{\mathrm{TL}}NM^{3})\right)$ for JD+NMF.

In complement of Remark 7, we see that the QN strategy deployed for TL-NMF is
preferable when $S\leq M$, whereas the one used for JD when $S\geq M$. Yet,
for the latter, the storage of the $N$ covariance matrices
$\\{\bm{\Sigma}_{n,S}\\}_{n\leq N}$ is required, which results in a memory
overload of $O(NM^{2})$.

##  V Numerical Experiments

From the theoretical analysis conducted in Section III, we know that the
solutions of TL-NMF and JD+NMF are getting closer as $S$ grows. More
precisely, the closeness between the solution sets of these two problems is
controlled by the inner and outer $\lambda$’s gaps (namely
$\underline{\lambda}^{\bullet}-\bar{\lambda}^{*}$ and
$\bar{\lambda}^{\bullet}-\underline{\lambda}^{*}$) which converge to zero with
high probability as $S$ grows (Proposition 2 and Theorem 3). Moreover, under
the GCM, they converge at least at a rate of $O(1/S)$ (Theorem 4).

Yet, these theoretical results only provide a partial answer to the questions
we identified in the introduction. For instance, the situation where $S$ is
small (including the extreme but very relevant case $S=1$) remains unclear. If
TL-NMF and JD+NMF are not equivalent, does one of them provide “better”
solutions (from an application-oriented perspective) than the other? On the
other hand, when they are equivalent (i.e., when the gap is approaching zero),
is one of them preferable from the point of view of computational complexity?
And lastly, how tight are the derived theoretical bounds? The goal of this
section is to provide numerical insights on these questions so as to
illustrate and complement our theoretical findings. The results presented in
this section can be reproduced using the Python code available at
https://github.com/sixin-zh/tlnmf-tsp.

###  V-A Datasets

We consider two simulated datasets obtained respectively from a pure GCM and a
blend of two synthetic music notes. The latter departs from GCM while roughly
satisfying the conditions (1).

####  V-A1 GCM

In this dataset, the dimensions are set as $M=10$, $N=50$, and $\bar{K}=5$.
The ground-truth transform ${\bar{\bm{\Phi}}}$ is fixed to the type-II DCT.
The entries of the ground-truth NMF factors $\bar{{\mathbf{W}}}$ and
$\bar{{\mathbf{H}}}$ are independently drawn from a Gamma distribution

$\displaystyle[\bar{\mathbf{W}}]_{mk}\sim\text{Gamma}(a,\theta),$ (49)
$\displaystyle[\bar{\mathbf{H}}]_{kn}\sim\text{Gamma}(a,\theta),$ (50)

with shape parameter $a=1$ and scale parameter $\theta=2$. The realizations
$\\{{\mathbf{Y}}^{(s)}\\}_{s=1}^{S}$ are generated such that (2) holds. Note
that, in the following experiments, when $S$ is changed the $S$ realizations
are re-sampled independently using the same $\bar{\mathbf{W}}$ and
$\bar{\mathbf{H}}$.

####  V-A2 Synthetic Music Notes

We follow the construction described in [5, Section III-A] which we recall for
completeness. For all $s\in\\{1,\ldots,S\\}$, we define the signal
${\mathbf{y}}^{(s)}\in{\mathbb{R}}^{T}$ as

${\mathbf{y}}^{(s)}[t]=\sum_{i=1}^{I=2}\sum_{r=1}^{R=2}0.5^{r}\cos\left(r\left(2\pi\frac{f_{i}}{f_{0}}t+\theta_{i}^{(s)}\right)\right)\mathbf{g}_{i}[t],$
(51)

where $\theta_{i}^{(s)}\in[0,2\pi)$ is a uniform random phase, $f_{i}>0$ is a
fundamental frequency, $f_{0}>0$ denotes the sampling frequency, $r$ is an
integer number, and $\mathbf{g}_{i}\in{\mathbb{R}}^{T}$ is a positive envelop
that varies slowly over $t\in\\{1,\ldots,T\\}$. In other words,
${\mathbf{y}}^{(s)}$ is the sum of $I=2$ pure music notes with frequencies
$f_{1}$ and $f_{2}$ that each contain $R=2$ harmonics. As in [5], we set the
fundamental frequencies to $f_{1}=440$Hz (corresponding to the note A4) and
$f_{2}=466.16$Hz (A4#), and fix $T=15000$ and $f_{0}=5000$Hz. This leads to a
signal of duration 3s. The envelops $\mathbf{g}_{1}$ and $\mathbf{g}_{2}$ are
such that the two notes are played separately in the first two seconds and
then simultaneously in the third second. Finally, the random nature of the
phases $\theta_{i}^{(s)}$ allows us to generate $S$ independent signals
$\\{{\mathbf{y}}^{(s)}\\}_{s=1}^{S}$.

As analyzed in [5, Section III-B], the underlying random matrix ${\mathbf{Y}}$
roughly follows the conditions (1). A subset of the atoms of the ground-truth
matrix ${\bar{\bm{\Phi}}}$ captures the fundamental frequencies
$\\{f_{1},f_{2}\\}$ and their harmonics $\\{2f_{1},2f_{2}\\}$ (two orthogonal
harmonic atoms per frequency to account for phase shifts). The ground-truth
NMF factors $\bar{\mathbf{W}}$ and $\bar{\mathbf{H}}$ capture the polyphonic
spectra and the activations of the two notes (hence $\bar{K}=2$).

TABLE I: Parameters used in the experiments: $K$ is the common dimension of the NMF factors (${\mathbf{W}}$,${\mathbf{H}}$); $\epsilon_{0}$ is defined in (4); $J$, $J_{\mathrm{TL}}$, and $J_{\mathrm{NMF}}$ are the numbers of iterations in Algorithms 1 and 2; $P$ denotes the number of random initializations (see Appendix G).  Dataset | GCM | Music Notes
---|---|---
$K$ |  $5$ ($=\bar{K}$) |  $2$ (see [5])
${\epsilon_{0}}$ | $10^{-8}$ | $5\times 10^{-7}$
$J_{\mathrm{TL}}$ | $1$ | $1$
$J_{\mathrm{NMF}}$ | $10$ | $10$
$J$ | $1000$ | $100$
$P$ | 100 | $10$

###  V-B Parameter Settings

The algorithm and model parameter values for all experiments are set to those
listed in Table I. Let us comment these choices. For the GCM dataset, we set
$K=\bar{K}$, and
${\epsilon_{0}}\ll\min_{m,n}[\bar{\mathbf{W}}\bar{\mathbf{H}}]_{mn}\simeq
1.57$. It respects the conditions of our theoretical results in Section III.
For the music notes dataset, we set $K=2$ to learn a NMF that can separate the
two music notes. Regarding ${\epsilon_{0}}$, we fixed it manually, as a trade-
off between numerical stability and quality of the solution in terms of
recovering the fundamental frequencies of the music notes.

The number of iterations $J_{\mathrm{TL}}$ and $J_{\mathrm{NMF}}$ have been
selected so as to maximize the convergence speed of TL-NMF, as they have no
effect on the convergence of JD+NMF due to its two-step nature. In particular,
we observed that TL-NMF enjoys a faster convergence when more NMF than TL
updates are performed per outer iteration.

Finally, because they address non-convex problems, we ran the algorithms from
$P$ different initializations and we retained the best solution (in terms of
objective function) for each problem. As such, we assume that we reached
points that belong to (or at least are close to) the solution sets
$\Omega^{*}$ and $\Omega^{\bullet}$. More details on our multi-initialization
strategy are provided in Appendix G.

###  V-C Analysis of the Gap

1101001000$10^{-6}$$10^{-2}$$10^{2}$$S$GCM
Dataset$\hat{I}^{\bullet}_{S}$$\hat{I}_{S}^{\ast}$$\hat{I}^{\bullet}_{S}-\hat{I}_{S}^{\ast}$$O\left(\frac{1}{S}\right)$$O\left(\frac{1}{S^{2}}\right)$1101001000$10^{-6}$$10^{-2}$$10^{2}$$S$Synthetic
Music Notes
Dataset$\hat{I}^{\bullet}_{S}$$\hat{I}_{S}^{\ast}$$\hat{I}^{\bullet}_{S}-\hat{I}_{S}^{\ast}$$O\left(\frac{1}{S}\right)$$O\left(\frac{1}{S^{2}}\right)$
Figure 2:  Evolution of the empirical quantities $\hat{I}^{\bullet}_{S}$,
$\hat{I}_{S}^{\ast}$, and $\hat{I}^{\bullet}_{S}-\hat{I}_{S}^{\ast}$ as
functions of the number of realizations $S$.

Next we will study the evolution of $\hat{I}^{\bullet}_{S}$,
$\hat{I}_{S}^{\ast}$, and $\hat{I}^{\bullet}_{S}-\hat{I}_{S}^{\ast}$ as
functions of the number of data realizations $S$. These quantities correspond
to the evaluation of the function $I_{S}$ at the solution points
$(\hat{{\bm{\Phi}}}^{\ast},{\hat{\mathbf{W}}}^{\ast},{\hat{\mathbf{H}}}^{\ast})$
and
$(\hat{{\bm{\Phi}}}^{\bullet},{\hat{\mathbf{W}}}^{\bullet},{\hat{\mathbf{H}}}^{\bullet})$
obtained through the numerical resolution of TL-NMF and JD+NMF, respectively.
Given our multi-initialization strategy, we may assume that
$(\hat{{\bm{\Phi}}}^{\ast},{\hat{\mathbf{W}}}^{\ast},{\hat{\mathbf{H}}}^{\ast})$
and
$(\hat{{\bm{\Phi}}}^{\bullet},{\hat{\mathbf{W}}}^{\bullet},{\hat{\mathbf{H}}}^{\bullet})$
belong to (or at least are close to) the solution sets $\Omega^{*}$ and
$\Omega^{\bullet}$. Note that from Theorem 2, Proposition 2, and by
construction of $\hat{I}_{S}^{\ast}$ and $\hat{I}^{\bullet}_{S}$, we expect
the following inequalities

$\underline{\lambda}^{*}\leq\hat{I}_{S}^{\ast}\leq\bar{\lambda}^{*}\leq\underline{\lambda}^{\bullet}\leq\hat{I}^{\bullet}_{S}\leq\bar{\lambda}^{\bullet}\leq\max_{{\bm{\Phi}}\in\mathcal{O}(M)}\,Q_{S}({\bm{\Phi}}).$
(52)

Figure 2 displays the values of $\hat{I}^{\bullet}_{S}$, $\hat{I}_{S}^{\ast}$,
and $\hat{I}^{\bullet}_{S}-\hat{I}_{S}^{\ast}$ for the GCM and music notes
datasets. For the GCM dataset, we observe that $\hat{I}_{S}^{\ast}$ decays as
$O(1/S)$ which implies that all the quantities in (52) (except
$\underline{\lambda}^{*}$) enjoy at best the same rate of convergence. In
particular, this constitutes an additional evidence that the bound derived in
Remark 5 is tight (also complementing Remark 6). It further shows that the
bound in Proposition 2 is also tight (under GCM). Concerning the gap
$\hat{I}^{\bullet}_{S}-\hat{I}_{S}^{\ast}$, we observe in Figure 2 that it
converges at the faster rate of $O(1/S^{2})$. From (25) and (27), this implies
that the convergence rate of the inner gap
$\underline{\lambda}^{\bullet}-\bar{\lambda}^{*}$ is at least of the order of
$O(1/S^{2})$ while the outer gap
$\bar{\lambda}^{\bullet}-\underline{\lambda}^{*}$ converges at a rate that
lies between $O(1/S^{2})$ and $O(1/S)$. We can draw two possible hypothesis
from this. The first one is that both gaps would enjoy the same rate of
convergence of $O(1/S^{2})$ which would mean that our upper bound in (27) is
loose (as we have seen that bounds in Proposition 2 and Theorem 4 are tights).
A direct analysis of $\bar{\lambda}^{\bullet}-\underline{\lambda}^{*}$ rather
than $\bar{\lambda}^{\bullet}$ would thus be of interest. The second possible
interpretation is that the two gaps would indeed not converge at the same rate
and thus that a finer analysis would require to consider both gaps separately.

Concerning the music notes dataset, we observe that all the reported
quantities are decaying up to $S\leq 100$. However, as $S$ further grows, the
quantity $\hat{I}^{\bullet}_{S}$ seems to tend toward a constant.

Although a theoretical explanation of this phenomenon remains an open problem,
it is likely due to the fact that the data do not exactly follow the
conditions (1). Nevertheless, the solutions returned by TL-NMF and JD+NMF on
this dataset remain very close when $S$ is large, as illustrated next.

TL-NMF$S=1$$S=100$JD+NMF Figure 3: Plots of the eight most significant atoms
(rows) of the learnt ${\bm{\Phi}}$ as well as the two columns of the estimated
${\mathbf{W}}$ and two rows of the estimated ${\mathbf{H}}$ for the music
notes dataset. The most significant atoms are chosen as those that maximize
${\mathbb{E}}_{S}\|\underline{\bm{\phi}}_{k}\mathbf{Y}\|_{2}^{2}$, where
$\underline{\bm{\phi}}_{k}$ is the $k^{th}$ row of ${\bm{\Phi}}$.

###  V-D Analysis of the Solutions

In Figure 3, we display the solutions obtained by TL-NMF and JD+NMF on the
music notes dataset for both $S=1$ and $S=100$. As expected, when $S$ is large
both approaches are able to learn a transform ${\bm{\Phi}}$ that captures the
exact fundamental and harmonic frequencies of the two musical notes, allowing
to significantly surpass the source separation performance of standard NMF
that must abide to arbitrary frequency grid of the chosen DCT or Fourier
frequency transform (see [5]).

However, the picture is quite different for $S=1$. Although the solution
provided by TL-NMF appears to be as good as for $S=100$ (i.e., it enjoys the
same favorable properties), the solution of JD+NMF is visibly degraded. To
further quantify the differences between the atoms learnt by TL-NMF and
JD+NMF, we performed a nonlinear least-square regression of the learnt atoms
with the harmonic model $(a,f,\theta)\mapsto
a\cos(2\pi\frac{f}{f_{0}}\cdot+\theta)$, as in [5]. We report in Table II the
estimated frequency $f$ for an atom $\underline{\bm{\phi}}\in\mathbb{R}^{M}$,
together with the regression error
$\|\underline{\bm{\phi}}-a\cos(2\pi\frac{f}{f_{0}}\cdot+\theta)\|^{2}$. We
observe that the eight most significant atoms of TL-NMF ($S=1$) lead to a
better fit to the fundamental frequencies of the two music notes due to
smaller regression errors. These errors could also explain why the NMF factors
${\mathbf{W}}^{\bullet}$ and ${\mathbf{H}}^{\bullet}$ are less regular in
JD+NMF. These observations are in line with Figure 2 where we see that, when
$S$ is small, the numerical gap $\hat{I}^{\bullet}_{S}-\hat{I}_{S}^{\ast}$ is
quite large, meaning that the transform ${\bm{\Phi}}^{*}$ learnt by TL-NMF
performs better in making
${\mathbb{E}}_{S}([{\bm{\Phi}}{\mathbf{Y}}]_{mn}^{2})$ more amenable to a low-
rank NMF approximation than the transform ${\bm{\Phi}}^{\bullet}$ obtained by
JD+NMF.

TABLE II: Frequencies and regression errors of the TL-NMF and JD+NMF atoms
when $S=1$ (same order as in Fig. 3).

TL-NMF
---
Freq. | Error
440.10 | 0.03
466.35 | 0.03
466.11 | 0.04
439.74 | 0.04
932.32 | 0.03
932.39 | 0.04
879.94 | 0.04
879.99 | 0.03

JD+NMF
---
Freq. | Error
466.12 | 0.19
439.72 | 0.37
879.97 | 0.51
879.99 | 0.45
879.97 | 0.44
932.46 | 0.42
932.28 | 0.40
932.31 | 0.39

###  V-E Computational Complexity

When JD+NMF provides solutions that are very close to those of TL-NMF (i.e.,
when $S$ is large), it is of interest to compare the computational load of the
two methods. To that end, we fix $J_{\mathrm{TL}}=1$, $J_{\mathrm{NMF}}=10$
and compute, for several values $J\in\mathbb{N}$,

  * •

$({\bm{\Phi}}^{*}_{J},{\mathbf{W}}^{*}_{J},{\mathbf{H}}^{*}_{J})$: the
solution obtained after $J$ iterations of TL-NMF (each composed of
$J_{\mathrm{TL}}$ updates of ${\bm{\Phi}}$ and $J_{\mathrm{NMF}}$ updates of
(${\mathbf{W}},{\mathbf{H}}$)),

  * •

$({\bm{\Phi}}^{\bullet}_{J},{\mathbf{W}}^{\bullet}_{J},{\mathbf{H}}^{\bullet}_{J})$:
the solution obtained after $(J\times J_{\mathrm{TL}})$ iterations of JD and
$(J\times J_{\mathrm{NMF}})$ NMF updates.

By doing so, given $J\in\mathbb{N}$, the number of cumulative updates of
${\bm{\Phi}}$ and $({\mathbf{W}},{\mathbf{H}})$ for both methods is the same
and the comparison is fair.

We report in Figure 4 the evolution of
$C_{S}({\bm{\Phi}}^{*}_{J},{\mathbf{W}}^{*}_{J},{\mathbf{H}}^{*}_{J})$ and
$C_{S}({\bm{\Phi}}^{\bullet}_{J},{\mathbf{W}}^{\bullet}_{J},{\mathbf{H}}^{\bullet}_{J})$
as a function of $J$ for $S=100$. While on the music notes dataset both
methods converge at a similar speed, we can see that JD+NMF converges faster
than TL-NMF on the GCM dataset. This suggests that, in the situation where
JD+NMF meets TL-NMF, the former should be preferred from the perspective of
the computational complexity.

$0$$20$$40$$60$$80$$100$$1{,}840$$1{,}860$$1{,}880$$J$$C_{S}$GCMTL-
NMFJD+NMF$0$$20$$40$$-4$$-3.9$$-3.8$$\cdot 10^{5}$$J$Music Notes Figure 4:
Evolution of $C_{S}$ (for $S=100$) as a function of the parameter $J$ in
Algorithms 1 and 2.

##  VI Discussion and Concluding Remarks

From our theoretical and numerical analysis, we can describe the relationship
between TL-NMF and JD+NMF by distinguishing two situations. On the one hand,
when $S$ is large and $K\geq\bar{K}$, the two problems are equivalent
(Theorems 2 and 3) and—at least on the reported numerical experiments—JD+NMF
seems to enjoy a faster convergence (Figure 4). On the other hand, when $S$ is
small (in particular $S=1$), they are not equivalent anymore and the solutions
obtained by TL-NMF appear to be more meaningful from an application point of
view (Figure 3). One explanation of this phenomenon comes from the inherent
low-rank constraint in the ${\bm{\Phi}}$ update of TL-NMF which favors the
learning of orthogonal transforms that are better suited to a low-rank NMF.
This is clearly visible when comparing the values of $I_{S}$ at the obtained
solutions (see Figure 2 and Theorem 2).

In this work, we fixed ${\mathbb{E}}_{S}(|{\bm{\Phi}}{\mathbf{Y}}|^{\circ 2})$
to the empirical expectation given by (10), which naturally led us to analyze
the relationship between the two problems as a function of the number of
realizations $S$. However, we can see from Proposition 2 that the main
ingredient to make TL-NMF and JD+NMF equivalent is to build a good estimator
of the variance ${\mathbb{E}}(|{\bm{\Phi}}{\mathbf{Y}}|^{\circ 2})$, in the
sense that the Itakura-Saito divergence between the estimator and the true
expectation is small (i.e., $Q_{S}({\bm{\Phi}})$ uniformly small). The
empirical expectation is a natural choice when several signal realizations are
available, but other choices are possible. Indeed, even for $S=1$, the
variance ${\mathbb{E}}(|{\bm{\Phi}}{\mathbf{Y}}|^{\circ 2})$ can be estimated
by a local moving average, under suitable local stationary assumptions. An
example is when ${\mathbf{H}}$ is given a block structure (with identical
columns inside the blocks), similarly to the block Gaussian model of [14].
Actually, such a block structure can be seen alternatively as having access to
several realizations of $\mathbf{Y}$ (with $S>1$ the size of the block).

In this paper we only considered real-valued transform for mathematical
convenience. Complex-valued transform such as the short-time Fourier transform
(STFT) are more customary in audio signal processing (in particular for shift-
invariance). Extending our work to such complex-valued transform is an
appealing research direction. Note also that the complex-valued Fourier
transform is an orthogonal transform endowed with a specific form of Hermitian
symmetry. In particular, when applied to real-valued signals, only half of the
power spectrum needs to be considered. Imposing such a structure to
$\bm{\Phi}$, to mimic some of the properties of the Fourier transform, would
be a very interesting topic. This is however a challenging problem, both in
terms of mathematical analysis and design of optimisation methods. Another
exciting research direction would be to lift the orthogonal constraint imposed
on $\bm{\Phi}$ and only assume invertibility. However, for TL-NMF to be useful
to signal processing tasks, it is important that ${\bm{\Phi}}$ can be inverted
for temporal reconstruction and thus that its inverse be well-conditioned.
Existing JD algorithms such that [24, 25] might be useful for this purpose.

##  Appendix A Proof of Lemma 1

The objective $L_{S}$ in (13) can be rewritten as

$\displaystyle L_{S}({\bm{\Phi}})$
$\displaystyle=MN+\sum_{n=1}^{N}\log\left(\prod_{m=1}^{M}\left({\mathbb{E}}_{S}\left([{\bm{\Phi}}{\mathbf{Y}}]_{mn}^{2}\right)+{\epsilon_{0}}\right)\right)$
(53)

Then from the definition of the empirical expectation in (10) and the
orthogonality of ${\bm{\Phi}}$, we have that for every $(m,n)$,

$\displaystyle{\mathbb{E}}_{S}\left([{\bm{\Phi}}{\mathbf{Y}}]_{mn}^{2}\right)+{\epsilon_{0}}$
$\displaystyle=\left[\frac{1}{S}\sum_{s=1}^{S}({\bm{\Phi}}\mathbf{y}_{n}^{(s)})({\bm{\Phi}}\mathbf{y}_{n}^{(s)})^{\intercal}+{\epsilon_{0}}{\mathbf{I}}\right]_{mm}$
$\displaystyle\underset{\eqref{eq:EmpiCov}}{=}\left[{\bm{\Phi}}(\bm{\Sigma}_{n,S}+{\epsilon_{0}}{\mathbf{I}}){\bm{\Phi}}^{\intercal}\right]_{mm}.$
(54)

Therefore combining the formula (54) with (53), we have

$L_{S}({\bm{\Phi}})=MN+\sum_{n=1}^{N}\log\det\operatorname*{Diag}\left({\bm{\Phi}}(\bm{\Sigma}_{n,S}+{\epsilon_{0}}{\mathbf{I}}){\bm{\Phi}}^{\intercal}\right).$

Finally, the JD criterion derived in [14], that is

$\sum_{n=1}^{N}\log\det\operatorname*{Diag}\left({\bm{\Phi}}\left(\bm{\Sigma}_{n,S}+{\epsilon_{0}}{\mathbf{I}}\right){\bm{\Phi}}^{\intercal}\right)-\\\
\log\det\left({\bm{\Phi}}\left((\bm{\Sigma}_{n,S}+{\epsilon_{0}}{\mathbf{I}}){\bm{\Phi}}^{\intercal}\right)\right),$
(55)

equals $L_{S}$ up to a constant term because, for ${\bm{\Phi}}$ orthogonal,
$\det({\bm{\Phi}}(\bm{\Sigma}_{n,S}+{\epsilon_{0}}{\mathbf{I}}){\bm{\Phi}}^{\intercal})=\det(\bm{\Sigma}_{n,S}+{\epsilon_{0}}{\mathbf{I}})$.

##  Appendix B Proof of Theorem 1

###  B-A Existence of a Solution for JD+NMF (Problem (6))

Given the two-step nature of Problem (6), the existence of a solution
$({\bm{\Phi}}^{\bullet},{\mathbf{W}}^{\bullet},{\mathbf{H}}^{\bullet})\in\Omega^{\bullet}$
is a direct consequence of the following two lemmas.

######  Lemma 2.

The set of global minimizers of $L_{S}$ over $\mathcal{O}(M)$ is nonempty and
compact.

###### Proof.

As $\epsilon_{0}>0$,
$\log({\mathbb{E}}_{S}([{\bm{\Phi}}{\mathbf{Y}}]_{mn}^{2})+{\epsilon_{0}})$ is
well defined and $L_{S}$ is continuous on the compact manifold
$\mathcal{O}(M)$. Therefore
$\inf_{{\bm{\Phi}}\in\mathcal{O}(M)}\;L_{S}({\bm{\Phi}})$ is attained and
finite (Weierstrass Theorem). Furthermore, the solution set is closed (it is a
level set of the continuous function $L_{S}$) and bounded because it is
included in the bounded set $\mathcal{O}(M)$. This proves the compactness and
completes the proof. ∎

######  Lemma 3.

For all ${\bm{\Phi}}\in\mathcal{O}(M)$, the set of global minimizers of
$I_{S}({\bm{\Phi}},\cdot,\cdot)$ over $F_{K}$ is nonempty and compact.

###### Proof.

By definition, for all ${\bm{\Phi}}\in\mathcal{O}(M)$,
$I_{S}({\bm{\Phi}},\cdot,\cdot)$ is continuous over $F_{K}$. As $F_{K}$ is
closed but unbounded, a sufficient argument to complete the proof is to show
that, for all ${\bm{\Phi}}\in\mathcal{O}(M)$, $I_{S}({\bm{\Phi}},\cdot,\cdot)$
is coercive over $F_{K}$, i.e.,

$\forall{\bm{\Phi}}\in\mathcal{O}(M),\;\lim_{\begin{subarray}{c}\|({\mathbf{W}},{\mathbf{H}})\|_{1,1}\rightarrow\infty\\\
({\mathbf{W}},{\mathbf{H}})\in
F_{K}\end{subarray}}I_{S}({\bm{\Phi}},{\mathbf{W}},{\mathbf{H}})=+\infty.$
(56)

Since all norms being equivalent in finite dimension, we use the norm
$\|({\mathbf{W}},{\mathbf{H}})\|_{1,1}=\|{\mathbf{W}}\|_{1,1}+\|{\mathbf{H}}\|_{1,1}$,
where $\|\cdot\|_{1,1}$ denotes the entrywise $\ell_{1}$-norm for matrices.
Because of the normalization constraint on ${\mathbf{W}}$ in $F_{K}$, we have

$\forall({\mathbf{W}},{\mathbf{H}})\in
F_{K},\;\|({\mathbf{W}},{\mathbf{H}})\|_{1,1}=K+\|{\mathbf{H}}\|_{1,1}.$ (57)

Now, let ${\bm{\Phi}}\in\mathcal{O}(M)$, $\nu>0$ and set

$\rho=KMN\mathrm{e}^{\nu+\log(D)+1}+K>K,$ (58)

where
$D=\max_{m,n}{\mathbb{E}}_{S}([{\bm{\Phi}}{\mathbf{Y}}]_{m,n}^{2})+{\epsilon_{0}}>0$.
Then we have the following implications:

$({\mathbf{W}},{\mathbf{H}})\in F_{K}\,\text{ and
}\,\|({\mathbf{W}},{\mathbf{H}})\|_{1,1}\geq\rho\\\ \quad\
\Longrightarrow\,\|{\mathbf{H}}\|_{1,1}\geq\rho-K>0\\\
\Longrightarrow\,\exists(k_{0},n_{0})\text{ s.t.
}h_{k_{0}n_{0}}\geq\frac{\rho-K}{KN},$ (59)

where the last implication comes from the fact that
${\mathbf{H}}\in{\mathbb{R}}_{+}^{K\times N}$. Similarly, using the fact that
$\mathbf{w}_{k_{0}}\in{\mathbb{R}}_{+}^{M}$ and
$\|\mathbf{w}_{k_{0}}\|_{1}=1$, we get that $\exists m_{0}$ s.t.
$w_{m_{0}k_{0}}\geq 1/M$. Combining the latter with (59), we obtain

$[{\mathbf{W}}{\mathbf{H}}]_{m_{0}n_{0}}\geq\frac{\rho-K}{MKN}\underset{\eqref{eq:Proof_ExistenceMinIs-3}}{\geq}\mathrm{e}^{\nu+\log(D)+1}.$
(60)

Finally, denoting $f(x)=x-\log(x)-1$ we have

$\displaystyle I_{S}({\bm{\Phi}},{\mathbf{W}},{\mathbf{H}})$
$\displaystyle=\sum_{m,n=1}^{M,N}f\left(\frac{{\mathbb{E}}_{S}([{\bm{\Phi}}{\mathbf{Y}}]_{mn}^{2})+\epsilon_{0}}{[{\mathbf{W}}{\mathbf{H}}]_{mn}+\epsilon_{0}}\right)$
$\displaystyle\geq
f\left(\frac{{\mathbb{E}}_{S}([{\bm{\Phi}}{\mathbf{Y}}]_{m_{0}n_{0}}^{2})+{\epsilon_{0}}}{[{\mathbf{W}}{\mathbf{H}}]_{m_{0}n_{0}}+\epsilon_{0}}\right)\quad\text{(as
$f\geq 0$ on ${\mathbb{R}}_{+}$)}$
$\displaystyle\geq-\log\left(\frac{{\mathbb{E}}_{S}([{\bm{\Phi}}{\mathbf{Y}}]_{m_{0}n_{0}}^{2})+{\epsilon_{0}}}{[{\mathbf{W}}{\mathbf{H}}]_{m_{0}n_{0}}+\epsilon_{0}}\right)-1$
$\displaystyle\geq\log([{\mathbf{W}}{\mathbf{H}}]_{m_{0}n_{0}}+\epsilon_{0})-\log(D)-1\underset{\eqref{eq:Proof_ExistenceMinIs-5}}{\geq}\nu.$

Hence, we have shown that for all ${\bm{\Phi}}\in\mathcal{O}(M)$ and all
$\nu>0$, there exists $\rho>0$ (i.e., (58)) such that

$({\mathbf{W}},{\mathbf{H}})\in F_{K}\text{ and
}\|({\mathbf{W}},{\mathbf{H}})\|_{1,1}\geq\rho\\\ \Longrightarrow
I_{S}({\bm{\Phi}},{\mathbf{W}},{\mathbf{H}})\geq\nu,$ (61)

which completes the proof. ∎

###  B-B Existence of a Solution for the TL-NMF Problem (3)

From Lemma 3, we get that the function

$O_{S}({\bm{\Phi}})=\min_{({\mathbf{W}},{\mathbf{H}})\in
F_{K}}\,I_{S}({\bm{\Phi}},{\mathbf{W}},{\mathbf{H}}),$ (62)

is well defined on $\mathcal{O}(M)$ in the sense that the $\min$ exists for
any ${\bm{\Phi}}\in\mathcal{O}(M)$. As such, from the decomposition of $C_{S}$
in (12) and given that $\mathcal{O}(M)$ is compact and $L_{S}$ continuous, it
is sufficient to prove that $O_{S}$ is continuous over $\mathcal{O}(M)$ (done
by Lemma 4) and invoke the Weierstrass Theorem to complete the proof that
$\Omega^{*}$ is non-empty.

######  Lemma 4.

The function $O_{S}$ defined in (62) is continuous over $\mathcal{O}(M)$.

###### Proof.

We need to show that, for any ${\bm{\Phi}}\in\mathcal{O}(M)$ and
${\mathbf{D}}\in\mathcal{T}_{\bm{\Phi}}$ (the tangent space of
$\mathcal{O}(M)$ at ${\bm{\Phi}}$),

$\lim_{\delta\rightarrow
0}\,|O_{S}(\pi({\bm{\Phi}}+\delta{\mathbf{D}}))-O_{S}({\bm{\Phi}})|=0,$ (63)

where $\pi$ stands for the projection operator on $\mathcal{O}(M)$ [37].

Let ${\bm{\Phi}}\in\mathcal{O}(M)$, ${\mathbf{D}}\in\mathcal{T}_{\bm{\Phi}}$
and $({\mathbf{W}}_{\delta},{\mathbf{H}}_{\delta})\in F_{K}$ be a global
minimizer of $I_{S}(\pi({\bm{\Phi}}+\delta{\mathbf{D}}),\cdot,\cdot)$ (there
exists at least one from Lemma 3). Then,

$\displaystyle
O_{S}(\pi({\bm{\Phi}}+\delta{\mathbf{D}}))-O_{S}({\bm{\Phi}})\leq
g_{\delta}({\mathbf{W}}_{0},{\mathbf{H}}_{0})$ (64) $\displaystyle
O_{S}(\pi({\bm{\Phi}}+\delta{\mathbf{D}}))-O_{S}({\bm{\Phi}})\geq
g_{\delta}({\mathbf{W}}_{\delta},{\mathbf{H}}_{\delta}),$ (65)

with
$g_{\delta}({\mathbf{W}},{\mathbf{H}})=I_{S}(\pi({\bm{\Phi}}+\delta{\mathbf{D}}),{\mathbf{W}},{\mathbf{H}})-I_{S}({\bm{\Phi}},{\mathbf{W}},{\mathbf{H}})$.
Moreover, from the definition of $I_{S}$ in (14)–(15), we get that for all
$({\mathbf{W}},{\mathbf{H}})\in F_{K}$,

$|g_{\delta}({\mathbf{W}},{\mathbf{H}})|\leq\sum_{m,n=1}^{M,N}\left|\log\left(\frac{{\mathbb{E}}_{S}[\pi({\bm{\Phi}}+\delta{\mathbf{D}}){\mathbf{Y}}]_{mn}^{2}+{\epsilon_{0}}}{{\mathbb{E}}_{S}[{\bm{\Phi}}{\mathbf{Y}}]_{mn}^{2}+{\epsilon_{0}}}\right)\right|\\\
+\frac{1}{{\epsilon_{0}}}\left|{\mathbb{E}}_{S}([\pi({\bm{\Phi}}+\delta{\mathbf{D}}){\mathbf{Y}}]_{mn}^{2})-{\mathbb{E}}_{S}([{\bm{\Phi}}{\mathbf{Y}}]_{mn}^{2})\right|.$
(66)

Because the right-hand side is independent of $({\mathbf{W}},{\mathbf{H}})$
and converges to $0$ as $\delta$ tends to $0$, we obtain

$\lim_{\delta\rightarrow
0}|g_{\delta}({\mathbf{W}}_{0},{\mathbf{H}}_{0})|=0,\quad\lim_{\delta\rightarrow
0}|g_{\delta}({\mathbf{W}}_{\delta},{\mathbf{H}}_{\delta})|=0.$ (67)

Finally, combining (67) with (64)–(65) proves (63). ∎

###  B-C Compactness of $\Omega^{\bullet}$

Given that $\Omega^{\bullet}$ is finite dimensional, it is sufficient to show,
according to Heine–Borel theorem, that it is a non-empty (already proved),
closed, and bounded set.

####  B-C1 Closedness

Let
$\\{({\bm{\Phi}}_{j},{\mathbf{W}}_{j},{\mathbf{H}}_{j})\in\Omega^{\bullet}\\}_{j}$
be a sequence of $\Omega^{\bullet}$ that converges toward
$(\widehat{{\bm{\Phi}}},\widehat{{\mathbf{W}}},\widehat{{\mathbf{H}}})$. Then,
from the continuity of $I_{S}$ and the fact that (by definition of
$\Omega^{\bullet}$),

$\mkern-10.0muI_{S}({\bm{\Phi}}_{j},{\mathbf{W}}_{j},{\mathbf{H}}_{j})\leq
I_{S}({\bm{\Phi}}_{j},{\mathbf{W}},{\mathbf{H}}),\;\forall({\mathbf{W}},{\mathbf{H}})\in
F_{K},$ (68)

we obtain

$I_{S}(\widehat{{\bm{\Phi}}},\widehat{{\mathbf{W}}},\widehat{{\mathbf{H}}})\leq
I_{S}(\widehat{{\bm{\Phi}}},{\mathbf{W}},{\mathbf{H}}),\;\forall({\mathbf{W}},{\mathbf{H}})\in
F_{K}.$ (69)

Noticing that
$\widehat{{\bm{\Phi}}}\in\operatorname*{arg\,min}_{{\bm{\Phi}}\in\mathcal{O}(M)}L_{S}({\bm{\Phi}})$
(compactness of the set of minimizers of $L_{S}$, Lemma 2), we conclude from
(69) that
$(\widehat{{\bm{\Phi}}},\widehat{{\mathbf{W}}},\widehat{{\mathbf{H}}})\in\Omega^{\bullet}$,
which shows that $\Omega^{\bullet}$ is closed.

####  B-C2 Boundedness

Let us assume that $\Omega^{\bullet}$ is unbounded. Hence, there exists a
sequence
$\\{({\bm{\Phi}}_{j},{\mathbf{W}}_{j},{\mathbf{H}}_{j})\in\Omega^{\bullet}\\}_{j}$
such that

$\displaystyle\lim_{j\rightarrow\infty}\|({\bm{\Phi}}_{j},{\mathbf{W}}_{j},{\mathbf{H}}_{j})\|_{1,1}=+\infty.$
(70)

Then, using the fact that $\mathcal{O}(M)$ is compact, we get that the norm of
${\bm{\Phi}}_{j}$ is bounded and thus, with (57) and (70), that
$\lim_{j\rightarrow\infty}\|{\mathbf{H}}_{j}\|_{1,1}=+\infty$. It follows from
the coercivity of $I_{S}$ (see Lemma 3) that

$\lim_{j\rightarrow\infty}O_{S}({\bm{\Phi}}_{j})\underset{\eqref{eq:Proof_ExistenceMinCs-1}}{=}\lim_{j\rightarrow\infty}I_{S}({\bm{\Phi}}_{j},{\mathbf{W}}_{j},{\mathbf{H}}_{j})=+\infty,$
(71)

which contradicts the fact that $O_{S}$ is a continuous function over the
compact set $\mathcal{O}(M)$ (see Lemma 4).

###  B-D Compactness of $\Omega^{*}$

Again, because $\Omega^{*}$ is finite dimensional, it is sufficient to show
that $\Omega^{*}$ is a non-empty (already proved), closed, and bounded set.
Here, we directly get that $\Omega^{*}$ is closed as a level set of the
continuous function $C_{S}$. Then, let us assume that $\Omega^{*}$ is
unbounded. Hence, using the same arguments as in the proof for
$\Omega^{\bullet}$, this means that there exists a sequence
$\\{({\bm{\Phi}}_{j},{\mathbf{W}}_{j},{\mathbf{H}}_{j})\in\Omega^{*}\\}_{j}$
for which $\lim_{j\rightarrow\infty}\|{\mathbf{H}}_{j}\|_{1,1}=+\infty$.
Hence, because

$C_{S}({\bm{\Phi}}_{j},{\mathbf{W}}_{j},{\mathbf{H}}_{j})\geq\sum_{m,n=1}^{M,N}\log\left([{\mathbf{W}}_{j}{\mathbf{H}}_{j}]_{mn}+{\epsilon_{0}}\right),$
(72)

we have

$\lim_{j\rightarrow\infty}C_{S}({\bm{\Phi}}_{j},{\mathbf{W}}_{j},{\mathbf{H}}_{j})=+\infty,$
(73)

which contradicts the fact that $\forall j$,
$({\bm{\Phi}}_{j},{\mathbf{W}}_{j},{\mathbf{H}}_{j})\in\Omega^{*}$.

##  Appendix C Proof of Theorem 2

First of all, the continuity of $I_{S}$ together with the compactness of
$\Omega^{*}$ and $\Omega^{\bullet}$ ensures the existence of the $\max$ and
$\min$ in (17)–(18) (Weierstrass Theorem). The following three proofs rely on
the key Equation (12) which states that $C_{S}=L_{S}+I_{S}$.

##### Proof of (19)

By definition, we have $\underline{\lambda}^{*}\leq\bar{\lambda}^{*}$ and
$\underline{\lambda}^{\bullet}\leq\bar{\lambda}^{\bullet}$. Moreover, the
lower bound $0$ in (19) is due to the fact that $I_{S}\geq 0$. It thus remains
to show that $\bar{\lambda}^{*}\leq\underline{\lambda}^{\bullet}$. Let us
assume that $\bar{\lambda}^{*}>\underline{\lambda}^{\bullet}$. Hence there
exist $({\bm{\Phi}}^{*},{\mathbf{W}}^{*},{\mathbf{H}}^{*})\in\Omega^{*}$ and
$({\bm{\Phi}}^{\bullet},{\mathbf{W}}^{\bullet},{\mathbf{H}}^{\bullet})\in\Omega^{\bullet}$
such that

$\displaystyle
I_{S}({\bm{\Phi}}^{*},{\mathbf{W}}^{*},{\mathbf{H}}^{*})>I_{S}({\bm{\Phi}}^{\bullet},{\mathbf{W}}^{\bullet},{\mathbf{H}}^{\bullet})$
(74) $\displaystyle\underset{\eqref{eq:decompCs}}{\Longleftrightarrow}\;$
$\displaystyle
C_{S}({\bm{\Phi}}^{*},{\mathbf{W}}^{*},{\mathbf{H}}^{*})>L_{S}({\bm{\Phi}}^{*})+I_{S}({\bm{\Phi}}^{\bullet},{\mathbf{W}}^{\bullet},{\mathbf{H}}^{\bullet})$
(75) $\displaystyle\Longrightarrow\;$ $\displaystyle
C_{S}({\bm{\Phi}}^{*},{\mathbf{W}}^{*},{\mathbf{H}}^{*})>C_{S}({\bm{\Phi}}^{\bullet},{\mathbf{W}}^{\bullet},{\mathbf{H}}^{\bullet}),$
(76)

where the last implication comes from the fact that
$L_{S}({\bm{\Phi}}^{*})\geq L_{S}({\bm{\Phi}}^{\bullet})$. This contradicts
the fact that $({\bm{\Phi}}^{*},{\mathbf{W}}^{*},{\mathbf{H}}^{*})$ is a
global minimizer of $C_{S}$ and completes the proof of (19).

##### Proof of (20) and (21)

The inclusion $\Longrightarrow$ is straightforward from (19) together with the
definition of $\bar{\lambda}^{*}$ and $\underline{\lambda}^{\bullet}$. Now, if
$\bar{\lambda}^{*}=\underline{\lambda}^{\bullet}$, then there exist
$({\bm{\Phi}}^{*},{\mathbf{W}}^{*},{\mathbf{H}}^{*})\in\Omega^{*}$ and
$({\bm{\Phi}}^{\bullet},{\mathbf{W}}^{\bullet},{\mathbf{H}}^{\bullet})\in\Omega^{\bullet}$
such that

$\displaystyle
I_{S}({\bm{\Phi}}^{*},{\mathbf{W}}^{*},{\mathbf{H}}^{*})=I_{S}({\bm{\Phi}}^{\bullet},{\mathbf{W}}^{\bullet},{\mathbf{H}}^{\bullet})$
(77) $\displaystyle\underset{\eqref{eq:decompCs}}{\Longleftrightarrow}\;$
$\displaystyle
C_{S}({\bm{\Phi}}^{*},{\mathbf{W}}^{*},{\mathbf{H}}^{*})=L_{S}({\bm{\Phi}}^{*})+I_{S}({\bm{\Phi}}^{\bullet},{\mathbf{W}}^{\bullet},{\mathbf{H}}^{\bullet})$
(78) $\displaystyle\Longrightarrow\;$ $\displaystyle
C_{S}({\bm{\Phi}}^{*},{\mathbf{W}}^{*},{\mathbf{H}}^{*})\geq
C_{S}({\bm{\Phi}}^{\bullet},{\mathbf{W}}^{\bullet},{\mathbf{H}}^{\bullet}).$
(79)

Given that $({\bm{\Phi}}^{*},{\mathbf{W}}^{*},{\mathbf{H}}^{*})$ is a global
minimizer of $C_{S}$, we have
$C_{S}({\bm{\Phi}}^{*},{\mathbf{W}}^{*},{\mathbf{H}}^{*})=C_{S}({\bm{\Phi}}^{\bullet},{\mathbf{W}}^{\bullet},{\mathbf{H}}^{\bullet})$.
Moreover, the latter equality with (77) lead to
$L_{S}({\bm{\Phi}}^{*})=L_{S}({\bm{\Phi}}^{\bullet})$. This completes the
proof.

##### Proof of (21)

can be obtained following the same steps.

##  Appendix D Proof of Proposition 2

First of all, one can easily see that

$\bar{\lambda}^{\bullet}\leq\max_{{\bm{\Phi}}\in\mathcal{O}(M)}\,O_{S}({\bm{\Phi}}),$
(80)

where $O_{S}$ is defined as in (62). To complete the proof, we shall show
that, $\forall{\bm{\Phi}}\in\mathcal{O}(M)$, $O_{S}({\bm{\Phi}})\leq
Q_{S}({\bm{\Phi}})$ where $Q_{S}$ is defined in (29).

To that end, let us show222This result was already shown in [5] without the
constraint set $F_{K}$. For completeness, we recall all the steps here. that,
for any ${\bm{\Phi}}\in\mathcal{O}(M)$, there exist
$({\mathbf{W}}_{\bm{\Phi}}^{\ast},{\mathbf{H}}_{\bm{\Phi}}^{\ast})\in F_{K}$
such that
${\mathbb{E}}([{\bm{\Phi}}{\mathbf{Y}}]_{mn}^{2})={\mathbf{W}}_{\bm{\Phi}}^{\ast}{\mathbf{H}}_{\bm{\Phi}}^{\ast}$.
Let $\mathbf{D}={\bm{\Phi}}{\bar{\bm{\Phi}}}^{\intercal}$,
$\bar{\mathbf{V}}=\bar{\mathbf{W}}\bar{\mathbf{H}}$ with
$\bar{v}_{mn}=[\bar{\mathbf{V}}]_{mn}$ and
$d_{mm^{\prime}}=[\mathbf{D}]_{mm^{\prime}}$. Then, from (7) and (8) we have

$\displaystyle{\mathbb{E}}([{\bm{\Phi}}{\mathbf{Y}}]_{mn}^{2})$
$\displaystyle={\mathbb{E}}\left(\left(\sum_{m^{\prime}=1}^{M}d_{mm^{\prime}}[{\bar{\bm{\Phi}}}{\mathbf{Y}}]_{m^{\prime}n}\right)^{2}\right)$
(81)
$\displaystyle=\sum_{m^{\prime}=1}^{M}d_{mm^{\prime}}^{2}\bar{v}_{m^{\prime}n}$
(82)
$\displaystyle=\sum_{k=1}^{\bar{K}}\sum_{m^{\prime}=1}^{M}d_{mm^{\prime}}^{2}\bar{w}_{m^{\prime}k}\bar{h}_{kn}$
(83)

For (81), we used the orthogonality of $\bar{\bm{\Phi}}$ to get
${\bm{\Phi}}=\mathbf{D}{\bar{\bm{\Phi}}}$. Defining
$({\mathbf{W}}_{\bm{\Phi}}^{\ast},{\mathbf{H}}_{\bm{\Phi}}^{\ast})\in{\mathbb{R}}^{M\times
K}_{+}\times{\mathbb{R}}^{K\times N}_{+}$ such that

$\displaystyle[{\mathbf{W}}_{\bm{\Phi}}^{\ast}]_{mk}=\left\\{\begin{array}[]{ll}\sum_{m^{\prime}}\frac{\bar{w}_{m^{\prime}k}}{\|\bar{\mathbf{w}}_{k}\|_{1}}d_{mm^{\prime}}^{2}&\text{
if }k\leq\bar{K},\\\ \frac{1}{M}&\text{ if }\bar{K}<k\leq
K,\end{array}\right.$ (86)
$\displaystyle[{\mathbf{H}}_{\bm{\Phi}}^{\ast}]_{kn}=\left\\{\begin{array}[]{ll}\bar{h}_{kn}\|\bar{\mathbf{w}}_{k}\|_{1}&\text{
if }k\leq\bar{K},\\\ 0&\text{ if }\bar{K}<k\leq K,\end{array}\right.$ (89)

we get from (83) that

${\mathbf{W}}_{\bm{\Phi}}^{\ast}{\mathbf{H}}_{\bm{\Phi}}^{\ast}={\mathbb{E}}(|{\bm{\Phi}}\mathbf{Y}|^{\circ
2}).$ (90)

Moreover, as $\sum_{m}d_{mm^{\prime}}^{2}=1$, we obtain that, for all
$k\leq\bar{K}$,

$\sum_{m=1}^{M}[{\mathbf{W}}_{\bm{\Phi}}^{\ast}]_{mk}=\sum_{m^{\prime}=1}^{M}\frac{\bar{w}_{m^{\prime}k}}{\|\bar{\mathbf{w}}_{k}\|_{1}}\sum_{m=1}^{M}d_{mm^{\prime}}^{2}=1.$
(91)

Given that, by definition in (86), we also have
$\sum_{m=1}^{M}[{\mathbf{W}}_{\bm{\Phi}}^{\ast}]_{mk}=1$ for $\bar{K}<k\leq
K$. Therefore we have shown that
$({\mathbf{W}}_{\bm{\Phi}}^{\ast},{\mathbf{H}}_{\bm{\Phi}}^{\ast})\in F_{K}$.

It follows from (90) that, $\forall{\bm{\Phi}}\in\mathcal{O}(M)$,

$\displaystyle O_{S}({\bm{\Phi}})$ $\displaystyle\leq
I_{S}({\bm{\Phi}},{\mathbf{W}}_{\bm{\Phi}}^{\ast},{\mathbf{H}}_{\bm{\Phi}}^{\ast})$
(92)
$\displaystyle=D_{\epsilon_{0}}\left({\mathbb{E}}_{S}(|{\bm{\Phi}}{\mathbf{Y}}|^{\circ
2})|{\mathbf{W}}_{\bm{\Phi}}^{\ast}{\mathbf{H}}_{\bm{\Phi}}^{\ast}\right)=Q_{S}({\bm{\Phi}}).$
(93)

##  Appendix E Proof of Theorems 3 and 4

###  E-A Preliminaries

In Lemma 5 we derive an equivalent formulation of condition (30) that is more
convenient to prove Theorems 3 and 4.

######  Lemma 5.

Condition (30) is satisfied if and only if ${U_{S}}\overset{p}{\to}0\;\text{
as }\;S\to\infty$ where

$\displaystyle{U_{S}}$
$\displaystyle=\max_{m,n}\max_{{\bm{\Phi}}\in\mathcal{O}(M)}\left|\frac{{\mathbb{E}}_{S}([{\bm{\Phi}}{\mathbf{Y}}]_{mn}^{2})+{\epsilon_{0}}}{{\mathbb{E}}([{\bm{\Phi}}{\mathbf{Y}}]^{2}_{mn})+{\epsilon_{0}}}-1\right|$
(94)
$\displaystyle\leq\max_{n}\|\bm{\Sigma}_{n}^{-1/2}\bm{\Sigma}_{n,S}\bm{\Sigma}_{n}^{-1/2}-\mathbf{I}\|.$
(95)

###### Proof.

Given the fact that $[\bar{\mathbf{W}}\bar{\mathbf{H}}]_{mn}>0$
$\forall(m,n)$, we introduce two constants $\epsilon_{1}>0$ and
$\epsilon_{2}>0$ such that for all $(m,n)$

$\epsilon_{1}\leq[\bar{\mathbf{W}}\bar{\mathbf{H}}]_{mn}\leq\epsilon_{2}.$
(96)

Then, combining (83) with $\sum_{m^{\prime}}d_{mm^{\prime}}^{2}=1$, one can
see that for all ${\bm{\Phi}}\in\mathcal{O}(M)$ and $(m,n)$,

$\epsilon_{1}\leq{\mathbb{E}}([{\bm{\Phi}}{\mathbf{Y}}]_{mn}^{2})\leq\epsilon_{2}.$
(97)

It follows that for any $\epsilon>0$

$\displaystyle\left|{\mathbb{E}}_{S}([{\bm{\Phi}}{\mathbf{Y}}]_{mn}^{2})-{\mathbb{E}}([{\bm{\Phi}}{\mathbf{Y}}]^{2}_{mn})\right|<\epsilon,$
(98) $\displaystyle\Longrightarrow\;$
$\displaystyle\left|\frac{{\mathbb{E}}_{S}([{\bm{\Phi}}{\mathbf{Y}}]^{2}_{mn})+{\epsilon_{0}}}{{\mathbb{E}}([{\bm{\Phi}}{\mathbf{Y}}]^{2}_{mn})+{\epsilon_{0}}}-1\right|<\frac{\epsilon}{\epsilon_{0}},$
(99)

and conversely

$\displaystyle\left|\frac{{\mathbb{E}}_{S}([{\bm{\Phi}}{\mathbf{Y}}]^{2}_{mn})+{\epsilon_{0}}}{{\mathbb{E}}([{\bm{\Phi}}{\mathbf{Y}}]^{2}_{mn})+{\epsilon_{0}}}-1\right|<\epsilon,$
(100) $\displaystyle\Longrightarrow\;$
$\displaystyle\left|{\mathbb{E}}_{S}([{\bm{\Phi}}{\mathbf{Y}}]_{mn}^{2})-{\mathbb{E}}([{\bm{\Phi}}{\mathbf{Y}}]^{2}_{mn})\right|<\epsilon(\epsilon_{2}+{\epsilon_{0}}).$
(101)

By definition of the convergence in probability, condition (30) is equivalent
to: $\forall\epsilon>0$ and $\delta\in(0,1)$, there exists $S^{\star}$ such
that $\forall S\geq S^{\star}$

$\mathrm{Pr}\left(\max_{m,n}\max_{{\bm{\Phi}}\in\mathcal{O}(M)}\left|{\mathbb{E}}_{S}([{\bm{\Phi}}{\mathbf{Y}}]_{mn}^{2})-{\mathbb{E}}([{\bm{\Phi}}{\mathbf{Y}}]^{2}_{mn})\right|<\epsilon\right)\\\
>1-\delta,$ (102)

which, using (98)–(99), implies that

$\mathrm{Pr}\left({U_{S}}<\frac{\epsilon}{\epsilon_{0}}\right)>1-\delta.$
(103)

This proves that (30) $\Longrightarrow\left({U_{S}}\overset{p}{\to}0\;\text{
as }\;S\to\infty\right)$. The converse can be proven similarly with
(100)–(101).

To complete the proof, it remains to show the inequality (95). First of all,
let us notice that, from (8) and (96), $\bm{\Sigma}_{n}$ is positive definite
and thus $\bm{\Sigma}_{n}^{-1/2}$ is well defined. Then, we have for all
$(m,n)$

$\displaystyle\max_{{\bm{\Phi}}\in\mathcal{O}(M)}\bigg{|}\frac{{\mathbb{E}}_{S}([{\bm{\Phi}}{\mathbf{Y}}]_{mn}^{2})+{\epsilon_{0}}}{{\mathbb{E}}([{\bm{\Phi}}{\mathbf{Y}}]^{2}_{mn})+{\epsilon_{0}}}-1\bigg{|}$
(104) $\displaystyle\underset{\eqref{eq:combin2}}{\leq}$
$\displaystyle\max_{\bm{\phi}\in\mathbb{R}^{M}:\|\bm{\phi}\|_{2}=1}\left|\frac{\bm{\phi}^{\intercal}(\bm{\Sigma}_{n,S}-\bm{\Sigma}_{n})\bm{\phi}}{\bm{\phi}^{\intercal}\bm{\Sigma}_{n}\bm{\phi}}\right|$
(105) $\displaystyle=$
$\displaystyle\max_{\tilde{\bm{\phi}}\in\mathbb{R}^{M}:\|\bm{\Sigma}_{n}^{-1/2}\tilde{\bm{\phi}}\|_{2}=1}\left|\frac{{\tilde{\bm{\phi}}}^{\intercal}(\bm{\Sigma}_{n}^{-1/2}\bm{\Sigma}_{n,S}\bm{\Sigma}_{n}^{-1/2}-\mathbf{I})\tilde{\bm{\phi}}}{\|\tilde{\bm{\phi}}\|_{2}^{2}}\right|$
(106) $\displaystyle=$
$\displaystyle\max_{\tilde{\bm{\phi}}\in\mathbb{R}^{M}:\|\tilde{\bm{\phi}}\|_{2}=1}\left|{\tilde{\bm{\phi}}}^{\intercal}(\bm{\Sigma}_{n}^{-1/2}\bm{\Sigma}_{n,S}\bm{\Sigma}_{n}^{-1/2}-\mathbf{I})\tilde{\bm{\phi}}\right|$
(107) $\displaystyle=$
$\displaystyle\;\|\bm{\Sigma}_{n}^{-1/2}\bm{\Sigma}_{n,S}\bm{\Sigma}_{n}^{-1/2}-\mathbf{I}\|,$
(108)

which completes the proof. ∎

In Lemma 6, we show that
$\max_{{\bm{\Phi}}\in\mathcal{O}(M)}Q_{S}({\bm{\Phi}})$ (defined in (29)) can
be controlled in terms of ${U_{S}}$ and that it converges to $0$ as ${U_{S}}$
tends to $0$.

######  Lemma 6.

Let ${U_{S}}$ be defined as in Lemma 5. Then,

${U_{S}}<1\Longrightarrow\max_{{\bm{\Phi}}\in\mathcal{O}(M)}Q_{S}({\bm{\Phi}})\leq
MN\frac{{U_{S}}^{2}}{1-{U_{S}}}.$ (109)

###### Proof.

Let us first recall that

$Q_{S}({\bm{\Phi}})=\sum_{m,n=1}^{M,N}f\left(\frac{{\mathbb{E}}_{S}([{\bm{\Phi}}{\mathbf{Y}}]^{2}_{mn})+{\epsilon_{0}}}{{\mathbb{E}}([{\bm{\Phi}}{\mathbf{Y}}]^{2}_{mn})+{\epsilon_{0}}}\right),$
(110)

for $f(x)=x-\log(x)-1$. Then, by definition of ${U_{S}}$ in (94) we have that,
for all ${\bm{\Phi}}\in\mathcal{O}(M)$ and $(m,n)$,

$\frac{{\mathbb{E}}_{S}([{\bm{\Phi}}{\mathbf{Y}}]^{2}_{mn})+{\epsilon_{0}}}{{\mathbb{E}}([{\bm{\Phi}}{\mathbf{Y}}]^{2}_{mn})+{\epsilon_{0}}}\in\left[1-{U_{S}},1+{U_{S}}\right].$
(111)

Hence, because $\forall x>0,f(x)\leq g(x):=x+1/x-2$, it follows from (110) and
(111) that, if ${U_{S}}<1$, $\forall{\bm{\Phi}}\in\mathcal{O}(M)$

$Q_{S}({\bm{\Phi}})\leq\sum_{m,n=1}^{M,N}g\left(\frac{{\mathbb{E}}_{S}([{\bm{\Phi}}{\mathbf{Y}}]^{2}_{mn})+{\epsilon_{0}}}{{\mathbb{E}}([{\bm{\Phi}}{\mathbf{Y}}]^{2}_{mn})+{\epsilon_{0}}}\right)$
(112)

and thus

$\displaystyle\max_{{\bm{\Phi}}\in\mathcal{O}(M)}Q_{S}({\bm{\Phi}})$
$\displaystyle\leq MN\max_{x\in[1-{U_{S}},1+{U_{S}}]}g(x)$ (113)
$\displaystyle=MN\max\\{g(1-{U_{S}}),g(1+{U_{S}})\\}$ (114)
$\displaystyle=MNg(1-{U_{S}})=MN\frac{{U_{S}}^{2}}{1-{U_{S}}}.$ (115)

Equality (114) comes from the fact that $g$ is a convex function and (115) is
due to $g(1-{U_{S}})\geq g(1+{U_{S}})$ for ${U_{S}}\in[0,1)$. ∎

###  E-B Proof of Theorem 3

The proof of Theorem 3 is a direct consequence of Lemmas 5 and 6. Indeed,
under condition (30), we get from Lemma 5 that, for all
$\epsilon\in(0,\epsilon_{0})$ and $\delta\in(0,1)$, there exists $S^{\ast}$
such that

$S\geq
S^{\ast}\quad\Longrightarrow\quad\mathrm{Pr}\left({U_{S}}<\frac{\epsilon}{\epsilon_{0}}\right)>1-\delta.$
(116)

(cf. (103)) which, with Lemma 6, completes the proof.

###  E-C Proof of Theorem 4

From Lemma 5, showing that condition (30) is satisfied under GCM amounts to
prove that under GCM we have ${U_{S}}\overset{p}{\to}0\;\text{ as
}\;S\to\infty$. Moreover, by making explicit a bound on the convergence rate
of ${U_{S}}$, we will get a bound on the convergence rate of
$\max_{{\bm{\Phi}}\in\mathcal{O}(M)}Q_{S}({\bm{\Phi}})$ thanks to Lemma 6.

According to Lemma 5 (Equation (95)), to show the convergence of ${U_{S}}$ in
probability, it is sufficient to show that each eigenvalue of
$\bm{\Sigma}_{n}^{-1/2}\bm{\Sigma}_{n,S}\bm{\Sigma}_{n}^{-1/2}$ converges to
one for each $n$.

To that end, let us first remark that the normalized matrix
$\bm{\Sigma}_{n}^{-1/2}\bm{\Sigma}_{n,S}\bm{\Sigma}_{n}^{-1/2}$ can be
regarded as an empirical covariance matrix obtained from the $S$ “whitened”
vectors
$\\{\mathbf{x}_{n}^{(s)}=\bm{\Sigma}_{n}^{-1/2}\mathbf{y}_{n}^{(s)}\\}_{s=1}^{S}$,
i.e.,

$\displaystyle\bm{\Sigma}_{n}^{-1/2}\bm{\Sigma}_{n,S}\bm{\Sigma}_{n}^{-1/2}=\frac{1}{S}\sum_{s=1}^{S}\mathbf{x}_{n}^{(s)}\mathbf{x}_{n}^{(s)\intercal}=\mathbf{A}_{S}^{\intercal}\mathbf{A}_{S},$
(117)

where
$\mathbf{A}_{S}=\frac{1}{\sqrt{S}}[\mathbf{x}_{n}^{1}\cdots\mathbf{x}_{n}^{S}]^{\intercal}\in{\mathbb{R}}^{S\times
M}$. Because the $\mathbf{x}_{n}^{(s)}$ are i.i.d. realizations of a whitened
Gaussian distribution whose covariance matrix is the identity, we obtain from
[31, Corollary 5.35] that, for any $t>0$

$\mbox{Pr}\left(1-\delta_{S}\leq\sigma_{\min}(\mathbf{A}_{S})\leq\sigma_{\max}(\mathbf{A}_{S})\leq
1+\delta_{S}\right)\geq 1-2e^{-t^{2}/2},$ (118)

where $\delta_{S}=\frac{\sqrt{M}+t}{\sqrt{S}}$ and
$\sigma_{\max}(\mathbf{A}_{S})$ (resp. $\sigma_{\min}(\mathbf{A}_{S})$)
denotes the largest (resp. smallest) singular value of $\mathbf{A}_{S}$. Using
[31, Lemma 5.36], we get from (118) that

$\mbox{Pr}\left(\|\mathbf{A}_{S}^{\intercal}\mathbf{A}_{S}-\mathbf{I}\|\leq
h_{S}\right)\geq 1-2e^{-t^{2}/2},$ (119)

with $h_{S}:=3\max(\delta_{S},\delta_{S}^{2})$. Then, from (117),

$\|\bm{\Sigma}_{n}^{-1/2}\bm{\Sigma}_{n,S}\bm{\Sigma}_{n}^{-1/2}-\mathbf{I}\|=\|\mathbf{A}_{S}^{\intercal}\mathbf{A}_{S}-\mathbf{I}\|.$
(120)

Combining the latter with (119), we obtain

$\mbox{Pr}\left(\|\bm{\Sigma}_{n}^{-1/2}\bm{\Sigma}_{n,S}\bm{\Sigma}_{n}^{-1/2}-\mathbf{I}\|\leq
h_{S}\right)\geq 1-2e^{-t^{2}/2}.$ (121)

It follows, using the fact that the vectors
$\\{\mathbf{x}_{n}^{(s)}\\}_{s=1}^{S}$ are independent of the vectors
$\\{\mathbf{x}_{n^{\prime}}^{(s)}\\}_{s=1}^{S}$ for $n\neq n^{\prime}$, that

$\displaystyle\mbox{Pr}\left({U_{S}}\leq h_{S}\right)$
$\displaystyle\geq\prod_{n=1}^{N}\mbox{Pr}\left(\|\bm{\Sigma}_{n}^{-1/2}\bm{\Sigma}_{n,S}\bm{\Sigma}_{n}^{-1/2}-\mathbf{I}\|\leq
h_{S}\right)$ $\displaystyle\geq(1-2e^{-t^{2}/2})^{N},$ (122)

which proves that, under GCM, ${U_{S}}\overset{p}{\to}0\;\text{ as
}\;S\to\infty$.

Finally, for $t>0$ and $S$ large enough such that $h_{S}=3\delta_{S}<1$, we
deduce from (122) and Lemma 6 that

Pr
$\displaystyle\left(\max_{{\bm{\Phi}}\in\mathcal{O}(M)}Q_{S}({\bm{\Phi}})\leq
MN\frac{h_{S}^{2}}{1-h_{S}}\right)$ (123)
$\displaystyle\geq\mbox{Pr}\left({U_{S}}\leq
h_{S}\right)\geq(1-2e^{-t^{2}/2})^{N}.$ (124)

##  Appendix F Proofs of Propositions 3 and 4

###  F-A Optimality Conditions of Problem (37)

In order to prove Propositions 3 and 4, we first explicit the optimality
conditions of the quadratic problem (37).

######  Lemma 7.

Let $\Upsilon$ be the set of anti-symmetric matrices and
$\tilde{{\mathcal{H}}}$ be positive semi-definite. Then,
${\mathbf{E}}\in\Upsilon$ is solution of problem (37) if and only if
${\mathbf{E}}=({\mathbf{M}}-{\mathbf{M}}^{\intercal})/2$ with
${\mathbf{M}}\in{\mathbb{R}}^{M\times M}$ solution of

$\frac{1}{4}\sum_{cd}[\mathcal{Z}]_{abcd}\left[\frac{{\mathbf{M}}-{\mathbf{M}}^{\intercal}}{2}\right]_{cd}=-[{\mathbf{G}}^{{(\mathrm{anti})}}]_{ab}$
(125)

where
$[\mathcal{Z}]_{abcd}=[\tilde{{\mathcal{H}}}]_{abcd}+[\tilde{{\mathcal{H}}}]_{cdab}+[\tilde{{\mathcal{H}}}]_{badc}+[\tilde{{\mathcal{H}}}]_{dcba}$.

###### Proof.

Using the parametrization of the anti-symmetric matrix
${\mathbf{E}}=({\mathbf{M}}-{\mathbf{M}}^{\intercal})/2$, solving the
constrained problem (37) is equivalent to solving the unconstrained problem

$\widehat{{\mathbf{M}}}=\operatorname*{arg\,min}_{{\mathbf{M}}\in{\mathbb{R}}^{M\times
M}}\;J({\mathbf{M}})$ (126)

with

$J({\mathbf{M}})=\left\langle{\mathbf{G}},\frac{{\mathbf{M}}-{\mathbf{M}}^{\intercal}}{2}\right\rangle+\frac{1}{2}\left\langle\frac{{\mathbf{M}}-{\mathbf{M}}^{\intercal}}{2}|\tilde{{\mathcal{H}}}|\frac{{\mathbf{M}}-{\mathbf{M}}^{\intercal}}{2}\right\rangle.$

where we recall that $\left\langle\cdot,\cdot\right\rangle$ and
$\left\langle\cdot|\cdot|\cdot\right\rangle$ are defined after equation (36).
Given that $J$ is convex (as $\tilde{{\mathcal{H}}}$ is positive semi-
definite), we get from the first-order optimality conditions that $\nabla
J(\widehat{{\mathbf{M}}})=\mathbf{0}$.

Because
${\langle}{\mathbf{G}},{\mathbf{M}}^{\intercal}{\langle}{\mathbf{G}}^{\intercal},{\mathbf{M}}$,
the first term of $J$ is equal to
${\langle}{\mathbf{G}}^{{(\mathrm{anti})}},{\mathbf{M}}$ and its gradient is
thus given by ${\mathbf{G}}^{{(\mathrm{anti})}}$. To compute the gradient of
the second term, let us first expand

$\displaystyle{\langle}{\mathbf{M}}+\bm{\Delta}|\tilde{{\mathcal{H}}}|{\mathbf{M}}+\bm{\Delta}$
$\displaystyle={\langle}{\mathbf{M}}|\tilde{{\mathcal{H}}}|{\mathbf{M}}o(\|\bm{\Delta}\|)$

Then, by definition of the inner product, we get

$\displaystyle{\langle}\bm{\Delta}|\tilde{{\mathcal{H}}}|{\mathbf{M}}\sum_{ab}[\bm{\Delta}]_{ab}\sum_{cd}[\tilde{{\mathcal{H}}}]_{abcd}[{\mathbf{M}}]_{cd},$
$\displaystyle{\langle}{\mathbf{M}}|\tilde{{\mathcal{H}}}|\bm{\Delta}\sum_{ab}[\bm{\Delta}]_{ab}\sum_{cd}[\tilde{{\mathcal{H}}}]_{cdab}[{\mathbf{M}}]_{cd}.$

Similar expressions can be obtained when transposing one or both variables of
the bilinear form in (LABEL:eq:Proof_QN_dir). We finally deduce from these
expansions that

$\displaystyle\left[\nabla{\langle}\cdot|\tilde{{\mathcal{H}}}|\cdot({\mathbf{M}})\right]_{ab}=\sum_{cd}([\tilde{{\mathcal{H}}}]_{abcd}+[\tilde{{\mathcal{H}}}]_{cdab})[{\mathbf{M}}]_{cd}$
$\displaystyle\left[\nabla{\langle}(\cdot)^{\intercal}|\tilde{{\mathcal{H}}}|\cdot({\mathbf{M}})\right]_{ab}=\sum_{cd}([\tilde{{\mathcal{H}}}]_{cdab}+[\tilde{{\mathcal{H}}}]_{badc})[{\mathbf{M}}^{\intercal}]_{cd}$
$\displaystyle\left[\nabla{\langle}\cdot|\tilde{{\mathcal{H}}}|(\cdot)^{\intercal}({\mathbf{M}})\right]_{ab}=\sum_{cd}([\tilde{{\mathcal{H}}}]_{abcd}+[\tilde{{\mathcal{H}}}]_{dcba})[{\mathbf{M}}^{\intercal}]_{cd}$
$\displaystyle\left[\nabla{\langle}(\cdot)^{\intercal}|\tilde{{\mathcal{H}}}|(\cdot)^{\intercal}({\mathbf{M}})\right]_{ab}=\sum_{cd}([\tilde{{\mathcal{H}}}]_{badc}+[\tilde{{\mathcal{H}}}]_{dcba})[{\mathbf{M}}]_{cd}$

Combining all these expressions, we obtain that $\nabla
J(\widehat{{\mathbf{M}}})=\mathbf{0}$ is equivalent to (125), which completes
the proof. ∎

###  F-B Proof of Proposition 3

By injecting the expression of
$[\tilde{\mathbf{{\mathcal{H}}}}]_{abcd}=\delta_{ac}\delta_{bd}[\bm{\Gamma}]_{ab}$
into (125) we get, for all $(a,b)$,

$[\bm{\Gamma}^{{(\mathrm{sym})}}]_{ab}\left[\frac{{\mathbf{M}}-{\mathbf{M}}^{\intercal}}{2}\right]_{ab}=-[{\mathbf{G}}^{{(\mathrm{anti})}}]_{ab}.$
(127)

If, for a couple $(a,b)$, we have $[{\bm{\Gamma}}^{{(\mathrm{sym})}}]_{ab}=0$,
we get from the definition of ${\bm{\Gamma}}$ in (42) that
$\forall(n,s),[{\mathbf{X}}^{(s)}]_{bn}=0$. It follows from (40) that
$[{\mathbf{G}}]_{ab}=0$ and thus that
$[{\mathbf{G}}^{{(\mathrm{anti})}}]_{ab}=0$. Hence, for all $(a,b)$ such that
$[{\bm{\Gamma}}^{{(\mathrm{sym})}}]_{ab}=0$, fixing
$\left[\frac{{\mathbf{M}}-{\mathbf{M}}^{\intercal}}{2}\right]_{ab}$ to any
real is solution of (127). In this work we choose the value $0$ (second line
in (43)). For all $(a,b)$ such that
$[{\bm{\Gamma}}^{{(\mathrm{sym})}}]_{ab}\neq 0$, we get from (127) the
expression given in the first line in (43).

Finally, let us show that the optimal solution is a descent-direction, i.e.,
${\langle}{\mathbf{G}},{\mathbf{E}}<0$. From the definition of
${\bm{\Gamma}}$, we get that $[{\bm{\Gamma}}^{{(\mathrm{sym})}}]_{ab}\geq 0$.
Moreover, for all $(a,b)$ such that
$[{\bm{\Gamma}}^{{(\mathrm{sym})}}]_{ab}=0$ we fixed $[{\mathbf{E}}]_{ab}=0$.
Therefore,

${\langle}{\mathbf{G}},{\mathbf{E}}{\langle}{\mathbf{G}}^{{(\mathrm{anti})}},{\mathbf{E}}\sum_{\begin{subarray}{c}(a,b)\\\
[{\mathbf{E}}]_{ab}\neq
0\end{subarray}}\frac{[{\mathbf{G}}^{{(\mathrm{anti})}}]_{ab}^{2}}{{\bm{\Gamma}}^{{(\mathrm{sym})}}_{ab}}<0.$
(128)

###  F-C Proof of Proposition 4

Again, by injecting the expression of

$[\tilde{{\mathcal{H}}}]_{abcd}=\delta_{ac}\delta_{bd}[{\bm{\Gamma}}]_{ab}+\delta_{ad}\delta_{bc}-2\delta_{abcd}.$
(129)

into (125) we obtain, for all $(a,b)$,

$([\bm{\Gamma}^{{(\mathrm{sym})}}]_{ab}-1)\left[\frac{{\mathbf{M}}-{\mathbf{M}}^{\intercal}}{2}\right]_{ab}=-[{\mathbf{G}}^{{(\mathrm{anti})}}]_{ab}.$
(130)

If, for a couple $(a,b)$, we have $[{\bm{\Gamma}}^{{(\mathrm{sym})}}]_{ab}=1$,
we get from the definition of ${\bm{\Gamma}}$ in (47) that $\forall n$,
$\frac{[{\bm{\Phi}}(\Sigma_{n,S}+{\epsilon_{0}}{\mathbf{I}}){\bm{\Phi}}^{\intercal}]_{bb}}{[{\bm{\Phi}}(\Sigma_{n,S}+{\epsilon_{0}}{\mathbf{I}}){\bm{\Phi}}^{\intercal}]_{aa}}=1$.
It follows from (45) that $[{\mathbf{G}}]_{ab}=[{\mathbf{G}}]_{ba}$ and thus
that $[{\mathbf{G}}^{{(\mathrm{anti})}}]_{ab}=0$. Hence, for all $(a,b)$ such
that $[{\bm{\Gamma}}^{{(\mathrm{sym})}}]_{ab}=1$, fixing
$\left[\frac{{\mathbf{M}}-{\mathbf{M}}^{\intercal}}{2}\right]_{ab}$ to any
real number is solution of (130). In this work we choose the value $0$ (second
line in (48)). For all $(a,b)$ such that
$[{\bm{\Gamma}}^{{(\mathrm{sym})}}]_{ab}\neq 1$, we get from (130) the
expression given in the first line in (48).

Finally, let us show that the optimal solution is a descent-direction, i.e.,
${\langle}{\mathbf{G}},{\mathbf{E}}<0$. From the definition of
${\bm{\Gamma}}$, we get that $[{\bm{\Gamma}}^{{(\mathrm{sym})}}]_{ab}\geq 1$.
Moreover, for all $(a,b)$ such that
$[{\bm{\Gamma}}^{{(\mathrm{sym})}}]_{ab}=1$ we fixed $[{\mathbf{E}}]_{ab}=0$.
Therefore,

${\langle}{\mathbf{G}},{\mathbf{E}}{\langle}{\mathbf{G}}^{{(\mathrm{anti})}},{\mathbf{E}}\sum_{\begin{subarray}{c}(a,b)\\\
[{\mathbf{E}}]_{ab}\neq
0\end{subarray}}\frac{[{\mathbf{G}}^{{(\mathrm{anti})}}]_{ab}^{2}}{{\bm{\Gamma}}^{{(\mathrm{sym})}}_{ab}-1}<0.$
(131)

##  Appendix G Multi-Initialization Strategy

We describe a multi-initialization strategy that allows to obtain numerical
solutions that are close to the solution sets $\Omega^{*}$ and
$\Omega^{\bullet}$. To that end, we first consider $P$ random initializations

$\\{({\bm{\Phi}}_{0}^{(p)},{\mathbf{W}}_{0}^{(p)},{\mathbf{H}}_{0}^{(p)})\\}_{p\leq
P}$ (132)

from which we run the TL-NMF and JD+NMF solvers in Algorithms 1 and 2. Given
an initialization, each method is run for $J$ iterations. In order to
robustify the process, the solutions obtained with JD+NMF (resp. TL+NMF)
$\\{({\bm{\Phi}}_{J}^{(p)},{\mathbf{W}}_{J}^{(p)},{\mathbf{H}}_{J}^{(p)})\\}_{p\leq
P}$ are used as extra $P$ initializations

$\\{({\bm{\Phi}}_{0}^{(p)},{\mathbf{W}}_{0}^{(p)},{\mathbf{H}}_{0}^{(p)})\\}_{P+1\leq
p\leq 2P}$ (133)

for TL-NMF (resp. JD+NMF). Finally, for TL-NMF, we preserve only the solution
$(\hat{{\bm{\Phi}}}^{\ast},{\hat{\mathbf{W}}}^{\ast},{\hat{\mathbf{H}}}^{\ast})$
that achieves the smallest $C_{S}$. For JD+NMF, we preserve only the solution
$\hat{{\bm{\Phi}}}^{\bullet}$ that achieves the smallest $L_{S}$. We then
minimize $I_{S}(\hat{{\bm{\Phi}}}^{\bullet},\cdot,\cdot)$, started from
$\\{({\mathbf{W}}_{0}^{(p)},{\mathbf{H}}_{0}^{(p)})\\}_{p\leq 2P}$, and
preserve the best solution
$({\hat{\mathbf{W}}}^{\bullet},{\hat{\mathbf{H}}}^{\bullet})$ that achieves
the smallest loss $I_{S}$.

## References

  * [1] A. Berman and R. J. Plemmons, _Nonnegative matrices in the mathematical sciences_. SIAM, 1994.
  * [2] C. Févotte, N. Bertin, and J. L. Durrieu, “Nonnegative matrix factorization with the Itakura-Saito divergence: With application to music analysis,” _Neural Computation_ , vol. 21, no. 3, pp. 793–830, mar 2009\.
  * [3] P. Smaragdis, C. Févotte, G. J. Mysore, N. Mohammadiha, and M. Hoffman, “Static and dynamic source separation using nonnegative factorizations: A unified view,” _IEEE Signal Processing Magazine_ , vol. 31, no. 3, pp. 66–75, may 2014.
  * [4] D. Fagot, H. Wendt, and C. Févotte, “Nonnegative matrix factorization with transform learning,” in _Proc. ICASSP_ , Calgary, Canada, 2018, pp. 2431–2435.
  * [5] S. Zhang, E. Soubies, and C. Févotte, “On the identifiability of transform learning for non-negative matrix factorization,” _IEEE Signal Processing Letters_ , vol. 27, pp. 1555–1559, 2020.
  * [6] K. Yoshii, K. Kitamura, Y. Bando, E. Nakamura, and T. Kawahara, “Independent low-rank tensor analysis for audio source separation,” in _European Signal Processing Conference_ , vol. 2018-Septe, Roma, Italy, 2018, pp. 1657–1661.
  * [7] E. Vincent, T. Virtanen, and S. Gannot, _Audio source separation and speech enhancement_. John Wiley & Sons, 2018.
  * [8] S. Makino, _Audio source separation_. Springer, 2018.
  * [9] N. Gillis, _Nonnegative Matrix Factorization_. SIAM, 2020.
  * [10] S. Ravishankar and Y. Bresler, “Learning Sparsifying Transforms,” _IEEE Transactions on Signal Processing_ , vol. 61, no. 5, pp. 1072–1086, mar 2013.
  * [11] A. Hyvärinen, J. Karhunen, and E. Oja, _Independent Component Analysis_. John Wiley & Sons, 2001.
  * [12] P. Comon and C. Jutten, _Handbook of blind source separation_. Academic press, 2010.
  * [13] X. Fu, K. Huang, N. D. Sidiropoulos, and W.-K. Ma, “Nonnegative Matrix Factorization for Signal and Data Analytics: Identifiability, Algorithms, and Applications,” _IEEE Signal Processing Magazine_ , vol. 36, no. 2, pp. 59–80, 2019.
  * [14] D. T. Pham and J.-F. Cardoso, “Blind separation of instantaneous mixtures of nonstationary sources,” _IEEE Transactions on Signal Processing_ , vol. 49, no. 9, pp. 1837–1848, sep 2001.
  * [15] F. Itakura and S. Saito, “Analysis synthesis telephony based on the maximum likelihood method,” in _International Congress on Acoustics_ , Los Alamitos, CA, USA, 1968, pp. 17–20.
  * [16] X. Fu, K. Huang, and N. D. Sidiropoulos, “On identifiability of nonnegative matrix factorization,” _IEEE Signal Processing Letters_ , vol. 25, no. 3, pp. 328–332, 2018.
  * [17] O. Ledoit and M. Wolf, “A well-conditioned estimator for large-dimensional covariance matrices,” _Journal of Multivariate Analysis_ , vol. 88, no. 2, pp. 365–411, feb 2004.
  * [18] J.-F. Cardoso and A. Souloumiac, “Jacobi angles for simultaneous diagonalization,” _SIAM Journal on Matrix Analysis and Applications_ , vol. 17, no. 1, pp. 161–164, jan 1996.
  * [19] J.-F. Cardoso, “On the performance of orthogonal source separation algorithms,” in _Proc. EUSIPCO_ , Edinburgh, Scotland, sep. 1994, pp. 776–779.
  * [20] D. T. Pham, “Joint approximate diagonalization of positive definite Hermitian matrices,” _SIAM Journal on Matrix Analysis and Applications_ , vol. 22, no. 4, pp. 1136–1152, may 2001.
  * [21] A. Yeredor, “Non-orthogonal joint diagonalization in the least-squares sense with application in blind source separation,” _IEEE Transactions on Signal Processing_ , vol. 50, no. 7, pp. 1545–1553, jul 2002.
  * [22] A. Ziehe, P. Laskov, G. Nolte, and K.-R. Müller, “A Fast Algorithm for Joint Diagonalization with Non-Orthogonal Transformations and Its Application to Blind Source Separation,” _Journal of Machine Learning Research_ , vol. 5, no. 1, pp. 777–800, jul 2004.
  * [23] P.-A. Absil and K. A. Gallivan, “Joint Diagonalization on the Oblique Manifold for Independent Component Analysis,” in _Proc. ICASSP_ , Toulouse, France, 2006.
  * [24] A. Souloumiac, “Nonorthogonal Joint Diagonalization by Combining Givens and Hyperbolic Rotations,” _IEEE Transactions on Signal Processing_ , vol. 57, no. 6, pp. 2222–2231, mar 2009.
  * [25] F. Bouchard, B. Afsari, J. Malick, and M. Congedo, “Approximate Joint Diagonalization with Riemannian Optimization on the General Linear Group,” _SIAM Journal on Matrix Analysis and Applications_ , vol. 41, no. 1, pp. 152–170, jan 2020.
  * [26] B. Afsari, “Sensitivity Analysis for the Problem of Matrix Joint Diagonalization,” _SIAM Journal on Matrix Analysis and Applications_ , vol. 30, no. 3, pp. 1148–1171, sep 2008.
  * [27] M. Kleinsteuber and H. Shen, “Uniqueness Analysis of Non-Unitary Matrix Joint Diagonalization,” _IEEE Transactions on Signal Processing_ , vol. 61, no. 7, pp. 1786–1796, jan 2013.
  * [28] N. Takahashi, J. Katayama, M. Seki, and J. Takeuchi, “A unified global convergence analysis of multiplicative update rules for nonnegative matrix factorization,” _Computational Optimization and Applications_ , vol. 71, no. 1, pp. 221–250, mar 2018.
  * [29] L. T. K. Hien and N. Gillis, “Algorithms for Nonnegative Matrix Factorization with the Kullback–Leibler Divergence,” _Journal of Scientific Computing_ , vol. 87, no. 3, p. 93, may 2021.
  * [30] C. Févotte and J. Idier, “Algorithms for nonnegative matrix factorization with the $\beta$-divergence,” _Neural Computation_ , vol. 23, no. 9, pp. 2421–2456, sep 2011.
  * [31] R. Vershynin, “Introduction to the non-asymptotic analysis of random matrices,” in _Compressed Sensing: Theory and Applications_. Cambridge University Press, 2012, pp. 210–268.
  * [32] ——, _High-Dimensional Probability: An Introduction with Applications in Data Science_. Cambridge University Press, 2018.
  * [33] Y. Sun, P. Babu, and D. P. Palomar, “Majorization-minimization algorithms in signal processing, communications, and machine learning,” _IEEE Transactions on Signal Processing_ , vol. 65, no. 3, pp. 794–816, 2017.
  * [34] P. Ablin, J.-F. Cardoso, and A. Gramfort, “Beyond Pham’s algorithm for joint diagonalization,” in _Proc. ESANN_ , Bruges, Belgium, 2019, pp. 607–612.
  * [35] P. Ablin, D. Fagot, H. Wendt, A. Gramfort, and C. Févotte, “A quasi-Newton algorithm on the orthogonal manifold for NMF with transform learning,” in _Proc. ICASSP_ , Brighton, United Kingdom, 2019, pp. 700–704.
  * [36] H. Wendt, D. Fagot, and C. Févotte, “Jacobi algorithm for nonnegative matrix factorization with transform learning,” in _European Signal Processing Conference_ , vol. 2018-Septe, Roma, Italy, 2018, pp. 1062–1066.
  * [37] J. H. Manton, “Optimization algorithms exploiting unitary constraints,” _IEEE Transactions on Signal Processing_ , vol. 50, no. 3, pp. 635–650, aug 2002.
  * [38] P. Wolfe, “Convergence Conditions for Ascent Methods,” _SIAM Review_ , vol. 11, no. 2, pp. 226–235, apr 1969.