# Efficient Algorithms for Regularized Nonnegative Scale-invariant Low-rank Approximation Models

Jérémy E. Cohen, Valentin Leplat

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

## To cite this version:

Jérémy E. Cohen, Valentin Leplat. Efficient Algorithms for Regularized Nonnegative Scale-invariant Low-rank Approximation Models. SIAM Journal on Mathematics of Data Science, In press. hal-04922736

HAL Id: hal-04922736

https://hal.science/hal-04922736v1

Submitted on 31 Jan 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

**Jeremy E. Cohen***
CREATIS UMR 5220, U1294
INSA Lyon, UCBL, CNRS, Inserm
Lyon, France
jeremy.cohen@cnrs.fr

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

**Valentin Leplat**
Faculty of Computer Science and Engineering
Innopolis University
Innopolis, Russia
V.Leplat@innopolis.ru

ABSTRACT

Regularized nonnegative low-rank approximations, such as sparse Nonnegative Matrix Factorization or sparse Nonnegative Tucker Decomposition, form an important branch of dimensionality reduction models known for their enhanced interpretability. From a practical perspective, however, selecting appropriate regularizers and regularization coefficients, as well as designing efficient algorithms, remains challenging due to the multifactor nature of these models and the limited theoretical guidance available. This paper addresses these challenges by studying a more general model, the Homogeneous Regularized Scale-Invariant model. We prove that the scale-invariance inherent to low-rank approximation models induces an implicit regularization effect that balances solutions. This insight provides a deeper understanding of the role of regularization functions in low-rank approximation models, informs the selection of regularization hyperparameters, and enables the design of balancing strategies to accelerate the empirical convergence of optimization algorithms.

Additionally, we propose a generic Majorization-Minimization (MM) algorithm capable of handling $\ell_p$-regularized nonnegative low-rank approximations with non-Euclidean loss functions, with convergence guarantees. Our contributions are demonstrated on sparse Nonnegative Matrix Factorization, ridge-regularized Nonnegative Canonical Polyadic Decomposition, and sparse Nonnegative Tucker Decomposition.

**Keywords** Constrained low-rank approximations, non-convex optimization, sparse approximation

# 1 Introduction

## 1.1 Context

Unsupervised machine learning heavily relies on Low-Rank Approximation (LRA) models, which are fundamental tools for dimensionality reduction. These models exploit the inherent multilinearity in datasets represented as matrices or tensors. One of the most commonly used LRA methods is Principal Component Analysis (PCA), which identifies orthogonal principal components (i.e., rank-one matrices) within a dataset. However, the orthogonality constraint can be detrimental when a physical interpretation of the principal components is required Harshman [1970], Comon and Jutten [2010].

Over the past two decades, researchers have explored various constrained LRA models with enhanced interpretability. Notably, Nonnegative Matrix Factorization (NMF) replaces PCA's orthogonality constraint with an elementwise nonnegativity constraint, allowing for a part-based decomposition of the data Lee and Seung [1999], Gillis [2020]. For tensorial data, models such as Nonnegative Canonical Polyadic Decomposition (NCPD) have received significant attention. The Nonnegative Tucker Decomposition (NTD), a multilinear dictionary-learning model, has also gained popularity Mørup et al. [2008]. Unconstrained LRA models face interpretation challenges due to inherent rotation

<footnote>*Jeremy Cohen is supported by ANR JCJC LoRAiA ANR-20-CE23-0010.</footnote>
ambiguities. Consequently, the interpretability of LRA models depends largely on the constraints or regularizations imposed on their parameters.

While the theoretical and algorithmic aspects of regularized Low-Rank Approximation have been extensively explored, several practical challenges remain. As regularized LRA are inherently multi-factor problems, these challenges include selecting appropriate regularization functions and hyperparameters for each parameter matrix, as well as developing efficient algorithms for non-Euclidean loss functions that offer convergence guarantees. In this work, we address these challenges by introducing a broader framework, which we term the Homogeneous Regularized Scale-Invariant (HRSI) problem. Through this framework, we examine the implicit regularization effects arising from the scale invariance of explicitly regularized LRA.

The following sections of this introduction formalize the general framework of our study, outline its current limitations, summarize our contributions, and provide an overview of the manuscript's structure.

## 1.2 Notations

Vectors, matrices and tensors are denoted $x$, $X$ and $\mathcal{X}$ respectively. The set $\mathbb{R}_+^{\times_{i=1}^d m_i}$ contains $d$-order tensors with dimensions $m_1 \times \dots \times m_d$ and nonnegative entries. Product $\mathcal{X} \times_i U$ is the n-mode product along the $i$-th mode between tensor $\mathcal{X}$ and matrix $U$. The Frobenius norm of matrix or tensor $X$ is denoted as $\|X\|_F$. The entry indexed by $(i, j)$ and the column indexed by $q$ in matrix $X$ are respectively written as $X[i, j]$ and $X[:, q]$. $\{X_i\}_{i \le n}$ is a set of matrices or tensors $X_i$ indexed by $i$ from 1 to $n$. The set of vectors $x$ such that $g(x) < +\infty$ is denoted by $\text{dom}(g)$, and $\eta_C$ is the characteristic function of convex set $C$. Matrix $I_r$ denotes the identity matrix of size $r \times r$, $x \otimes y$ the tensor (outer) product, $A \odot B$ and $A \boxtimes B$ respectively the Hadamard and Kronecker product between two tensors $A$ and $B$ and $X^\alpha$ the elementwise power of tensor $X$ with exponent $\alpha$. More details on multi-linear algebra can be found in Kolda and Bader [2009].

## 1.3 The Homogeneous Regularized Scale-Invariant model

Consider the following optimization problem

$$ \underset{\forall i \le n, X_i \in \mathbb{R}^{m_i \times r}}{\operatorname{argmin}} \quad f(\{X_i\}_{i \le n}) + \sum_{i=1}^{n} \mu_i \sum_{q=1}^{r} g_i(X_i[i, q]), \quad (1) $$

where $f$ is a continuous map from the Cartesian product $\times_{i=1}^n \mathbb{R}^{m_i \times r} \cap \text{dom}(g_i)$ to $\mathbb{R}_+$, $\{\mu_i\}_{i \le n}$ is a set of positive regularization hyperparameters, and $\{g_i\}_{i \le n}$ is a set of lower semi-continuous regularization maps from $\mathbb{R}^{m_i}$ to $\mathbb{R}_+$. We assume that the total cost is coercive<fnref n="2" />, ensuring the existence of a minimizer. Furthermore, we assume the following assumptions hold:

**A1** The function $f$ is invariant to column-wise scaling of the parameter matrices. For any sequence of positive diagonal matrices $\{\Lambda_i\}_{i \le n}$,

$$ f(\{X_i\}_{i \le n}) = f(\{X_i\}_{i \le n}) \quad \text{if} \quad \prod_{i=1}^{n} \Lambda_i = I_r, \quad \Lambda_i > 0. \quad (2) $$

In other words, scaling any two columns $X_{i_1}[\cdot, q]$ and $X_{i_2}[\cdot, q]$ inversely proportionally has no effect on the value of $f$.

**A2** Each regularization function $g_i$ is positive homogeneous of degree $p_i$, that is for all $i \le n$ there exists $p_i > 0$ such that for any $\lambda > 0$ and any $x \in \mathbb{R}^{m_i}$,

$$ g_i(\lambda x) = \lambda^{p_i} g_i(x). \quad (3) $$

This property holds in particular for any $\ell_p$ norm in the ambient vector space.

**A3** Each function $g_i$ is positive-definite, meaning that for any $i \le n$, $g_i(x) = 0$ if and only if $x$ is null.<fnref n="3" />

We refer to the approximation problem (1) with assumptions **A1**, **A2** and **A3** the Homogeneous Regularized Scale-Invariant problem (HRSI).

<footnote n="2">²Coercivity is satisfied by the regularization functions considered in this work.</footnote>

<footnote n="3">³This assumption limits the framework's applicability, excluding positive homogeneous regularizations like Total Variation.</footnote>
HRSI is quite general and encompasses many instances of regularized LRA problems. To illustrate this, consider the computation of solutions to a (double) sparse NMF (sNMF) model, where an $l_1$-norm regularization promotes sparsity in the matrices $X_1$ and $X_2$ and the Kullback-Leibler divergence serves as the loss function. This problem can be formalized as the following optimization problem:

$$ \underset{X_1 \ge 0, X_2 \ge 0}{\operatorname{argmin}} \mathrm{KL}(M, X_1 X_2^T) + \mu_1 \|X_1\|_1 + \mu_2 \|X_2\|_1, \quad (4) $$

where $M \in \mathbb{R}^{m_1 \times m_2}$ is the input data matrix. This formulation for sNMF fits within the HRSI framework because the cost function is scale-invariant and the regularizations are homogeneous, positive-definite, and columnwise separable. Specifically, sNMF as defined in Equation(4) is an HRSI problem with $f(\{X_i\}_{i \le n}) = \mathrm{KL}(M, X_1 X_2^T)$, $g_1(x) = \|x\|_1 + \eta_{\mathbb{R}^{m_1}}^*(x)$ and $g_2(x) = \|x\|_1 + \eta_{\mathbb{R}^{m_2}}^*(x)$ where $\|x\|_1$ denotes the $l_1$ norm in $\mathbb{R}^{m_1}$ applied to each column of $X_i$.

## 1.4 Current limitations

Let us consider conventional variational problems, focusing on the well-known LASSO problem Tibshirani [2011]. A single regularization hyperparameter is used for this problem, although more complex models have been proposed Bogdan et al. [2015]. Representer theorems state that solutions to such problems are combinations of extreme points within the regularization unit ball Boyer et al. [2019], Unser [2021]. A maximum value for the regularization hyperparameter exists beyond which all solutions become degenerate Tibshirani [2013], and it is simple to compute. Additionally, homotopy methods establish a connection between the regularity of solutions and the values of the regularization hyperparameters Efron et al. [2004]. In simpler terms, the relationship between the penalty function choice, the regularization hyperparameter, and the solution set is well understood in standard cases. For scale-invariant multifactor models however, there are no known representer theorems or maximum regularization hyperparameter characterizing the regularity of the solutions to the best of our knowledge. The presence of multiple regularization hyperparameters also complicates model selection and evaluation.

Furthermore, the design of algorithms for HRSI with non-Euclidean loss functions — particularly in the context of nonnegative parameter matrices and beta-divergence loss Cichocki et al. [2006], Févotte and Idier [2011] — has been underexplored. The design of efficient optimization algorithms for NMF with non-Euclidean loss is a relatively mature area. Many algorithms, including the popular Majorization-Minimization (MM) algorithm, have been proposed to solve the unregularized problem Lee and Seung [1999], Févotte and Idier [2011], but these methods do not trivially extend for regularized factorization. Many regularized LRA models rely on normalization constraints to simplify the tuning of multiple hyperparameters and mitigate degeneracy (see Section 2.1). This approach makes the problem more challenging to solve due to nonconvexity issues. To the best of our knowledge, MM algorithms have no known convergence guarantees in this context.

## 1.5 Contributions and Impact

Our main contribution is to show that solutions to the HRSI problem inherently benefit from implicit regularization due to their scale-invariance. We show the following:

* Solutions to HRSI are balanced with respect to the regularization functions, as detailed in Section 2.2.

* Explicit regularization on matrices $X_i$ induces implicit regularization on the rank-one components, providing insight into the properties of the solutions. For instance, applying ridge regularization to all parameter matrices in a ridge CPD implicitly enforces a group LASSO regularization on the rank-one components, revealing the tensor CP rank.

* The hyperparameters for each matrix $X_i$ can be chosen simultaneously, as discussed in Section 2.3.

These results extend the known properties of specific regularized LRA models. For further details, see Section 2.2.2 for two examples, and Section 2.4 for an overview of existing results.

On the front of algorithm design, focusing on the case where the loss function $f$ is a $\beta$-divergence, we make the following contributions.

* A general alternating MM Meta-algorithm is introduced in Section 3 which provably computes a stationary point of the regularized LRA problem and accommodates various regularization functions (in particular any $\ell_p^p$ regularization). In particular, our proposed MM update rules for various combinations of $\beta$-divergences and $\ell_p^p$ regularizations are an original contribution, as they guarantee convergence. To the best of our knowledge, existing works typically propose heuristic adaptations of unregularized MM updates in the presence of regularizations, but without convergence guarantees for the iterates Le Roux and Hershey [2015].
* This Meta-algorithm features a balancing strategy that scales the columns of the parameter matrices to minimize the regularization terms. We show formally on a toy problem and experimentally on several practical examples that this procedure improves the convergence speed for the iterates, which otherwise may be sublinear.

We illustrate these implicit regularization effects and the proposed Meta-algorithm through the lens of three models: sNMF, ridge-regularized NCPD (rNCPD), and sparse NTD (sNTD). These models are tested on simple synthetic and more complex real-world datasets in Section 4.

These contributions have the potential to improve the usability of regularized LRA models across a wide range of applications, such as spectral unmixing in remote sensing Bioucas-Dias et al. [2012] or medical microscopy Caredda et al. [2023], music information retrieval Marmoret et al. [2020], source separation in chemometric Bro [1998], neu- roimaging Roald et al. [2022] and antenna array processing Sidiropoulos et al. [2000]. See Gillis [2020], Kolda and Bader [2009] for surveys about matrix and tensor decomposition applications, respectively. In these applications, regularization is often crucial for ensuring or enhancing the interpretability of the results. While CPD enjoys uniqueness properties under mild conditions Kruskal [1977], both NMF and NTDare generally not unique without further regular- ization Cichocki et al. [2009], Gillis [2020]. Homogeneous regularizations also reshape the optimization landscape by making the cost function coercive, thereby avoiding issues such as swamps that can arise in tensor decompositions Lim and Comon [2009], Mohlenkamp [2019]. This work guides how to choose such regularization functions, select individual regularization hyperparameters, and design efficient optimization algorithms with convergence guarantees for a broad range of loss functions, accommodating various noise models. Finally, we show that in many scenarios, explicit regularization of all matrices in regularized LRA leads to group sparsity on the components. This helps in selecting the rank, a significant practical issue with these models. This work also includes Python code<fnref n="4" /> providing off-the-shelf algorithms for sNMF, rCPD and sNTD.

## 2 Scale invariance in HRSI leads to implicit balancing

We begin by discussing the ill-posedness of HRSI, emphasizing the importance of regularizing all the modes. After addressing this issue, we move on to our main contribution: demonstrating that the scale-invariance in the cost function leads to optimal scaling with respect to the regularization, effectively balancing the solutions. An underlying optimization problem emerges that is essentially equivalent<fnref n="5" /> to HRSI, exhibiting an implicit regularization effect. This implicit HRSI formulation provides valuable insights for selecting regularization functions {$g_i$}$_{i \le n}$ and regularization hyperparameters.

### 2.1 Ill-posedness of partially regularized homogeneous scale-invariant problems

In Equation (1) which defines the HRSI problem, we assumed that the regularization hyperparameters are nonnegative. It is known in the literature Marmin et al. [2023] that if the parameter $\mu_j$ for a specific index $j$ is zero, the scale- indeterminacy of the loss leads to degenerate solutions for HRSI. A stronger result is shown in this work: a partially regularized HRSI problem yields the same optimal cost as an unregularized problem, but with unattainable solutions. In other words, the regularization functions do not change the minimum cost unless all parameter matrices are regularized.

**Proposition 1.** Let $\mu_j = 0$ for at least one index $j \le n$. Then the infimum of the function

$$f(\{X_i\}_{i \le n}) + \sum_{i=1}^{n} \mu_i \sum_{q=1}^{r} g_i(X_i[:q]) \quad (5)$$

is equal to $\inf_{X_i \in \mathbb{R}^{m_i \times r}} f(\{X_i\}_{i \le n})$. Moreover, the infimum is not attained unless $f$ attains its minimum at $X_i = 0$ for all $i \le n$.

A proof is provided in Appendix A. In the following, we assume that all matrices $X_i$ are regularized, with associated regularization hyperparameter $\mu_i > 0$.

### 2.2 Optimal balancing of parameter matrices in HRSI

We now turn to the central discussion of this section: the effect of scaling invariance on regularization. Since the loss function $f$ is scale-invariant, it follows that one can minimize the whole cost by adjusting the scaling of the columns

<footnote n="4">⁴https://github.com/vleplat/NTD-Algorithms</footnote>

<footnote n="5">⁵We use “essential equivalence” to describe two optimization problems where any solution to one is also a solution to the other, up to a scaling that complies with assumption A1.</footnote>
with respect to the regularization functions only. A closed-form solution for the optimal scaling derives from the homogeneity and the positive-definiteness of the regularization functions. This solution offers valuable insights into the implicit regularizations and invariances present within the HRSI framework.

### 2.2.1 Optimal scaling of regularization terms

Let {$X_i$}$_{i \le n}$ be a set where all $X_i$ have nonzero columns, and $q \le r$ be a column index. To simplify the notation, we let $a_i = \mu_i g_i(X_i[;, q])$ where $a_i > 0$ due to assumption A3. Consider the optimization problem obtained by considering the HRSI problem (1) with respect to only the scales {$\lambda_i$}$_{i \le n}$ of the q-th columns of all matrices $X_i$,

$$ \min_{\forall i \le n, \lambda_i \ge 0} \sum_{i=1}^{n} \lambda_i^{p_i} a_i \quad \text{such that} \quad \prod_{i=1}^{n} \lambda_i = 1. \qquad (6) $$

The scale-invariant loss function $f$ and the regularizations over other columns $g_i(X_i[;, p])$ for $p \ne q$, which are constant with respect to a scaling of the q-th columns represented by the parameters {$\lambda_i$}$_{i \le n}$, have been removed. Our objective is to determine the optimal scales.

One can observe that Problem (6) admits a unique optimal solution. The optimal scales can provably be computed from the weighted geometrical mean $\beta$ of {$p_i a_i, \frac{1}{p_i}$} as follows:

$$ \lambda_i^* = \left( \frac{\beta}{p_i a_i} \right)^{\frac{1}{p_i}}, \quad \text{with } \beta = \left( \prod_{i \le n} (p_i a_i)^{\frac{1}{p_i}} \right)^{\frac{1}{\sum_{i \le n} \frac{1}{p_i}}}. \qquad (7) $$

Here, the set {$\lambda_i^*$}$_{i \le n}$ contains the solutions to the scaling problem (6).

Since at optimality the balancing procedure should have no effect, given an optimal solution $X_i^*$ to HRSI, it holds that $\lambda_i^* = 1$ for all $i \le n$, where $a_i = \mu_i g_i(X_i[;, q^*])$. In that case, we can infer that all additive terms in (6) must be equal up to a multiplicative factor. Specifically, for any $i \le j \le n$, we have $p_i a_i = p_j a_j = \beta$. We can now state the following proposition characterizing the solutions of HRSI.

**Proposition 2.** If $\mu_i > 0$ for all $i$, any solution {$X_i^*$}$_{i \le n}$ to Problem (1) satisfies $p_i \mu_i g_i(X_i^*[, q]) = \beta_q$ for all $i \le n$ and $q \le r$, where

$$ \beta_q = \left( \prod_{i \le n} (p_i \mu_i g_i(X_i^*[, q]))^{\frac{1}{p_i}} \right)^{\frac{1}{\sum_{i \le n} \frac{1}{p_i}}}. \qquad (8) $$

A proof for Proposition 2 is provided in Appendix B. If all regularization hyperparameters $\mu_i$ are nonzero but a column $X_{i_0}[;, q]$ is null for a given index $i_0 \le n$, by scaling to zero columns $X_i[;, q]$ for $i \ne i_0$, the penalty functions are trivially minimized with a constant cost. Therefore, the optimal balancing in the presence of a null column $X_{i_0}[;, q]$ is to set the corresponding columns $X_{i \ne i_0}[;, q]$ to zero as well. Hence, Proposition 2 allows for solutions with zero columns.

### 2.2.2 Balancing algorithm and essentially equivalent scale-free implicit HRSI

In the previous section, we demonstrated that it is possible to compute the optimal scaling of parameter matrices to minimize the regularization terms while keeping the cost function $f$ constant. In simpler terms, at any specific point within the parameter space, we can reduce the cost function in (1) by applying a scaling procedure, which we refer to as the balancing procedure.

Formally, given a set of parameter matrices {$X_i$}$_{i \le n}$ with nonzero columns, and assuming that all regularization hyperparameters $\mu_i$ are positive, the balancing procedure is defined as follows:

$$ \tilde{X}_i[;, q] = \left( \frac{\beta_q}{p_i \mu_i g_i(X_i[;, q])} \right)^{1/p_i} X_i[;, q], \qquad (9) $$

where $\beta_q$ is given by Proposition 2. The balancing procedure results in a new set of parameters {$\tilde{X}_i$}$_{i \le n}$ which solves the optimal scaling problem

$$ \min_{\forall i \le n, \Lambda_i \in \mathbb{R}_+^{r \times r}} f(\{X_i \Lambda_i\}_{i \le n}) + \sum_{i=1}^{n} \mu_i \sum_{q=1}^{r} g_i(\Lambda_i[q, q] X_i[;, q]) \quad \text{s.t. } \Lambda_i \text{ are diagonal and } \prod_{i \le n} \Lambda_i = I_r. \qquad (10) $$

The balancing procedure is summarized in Algorithm 1.
**Algorithm 1 Balancing algorithm**

**Require:** A collection of columns {$X_i[;, q]}_{i \le n}$ for a fixed index $q \le r$, regularization hyperparameters $\mu_i > 0$, regularization functions $g_i$ and their homogeneity degree $p_i$.

1: **if** $X_i[;, q] = 0$ for any $i \le n$ **then**

2: Set $X_i[;, q] = 0$ for all $i \le n$.

3: **else**

4: Compute the geometric mean $\beta_q = \left( \prod_{i \le n} (p_i \mu_i g_i(X_i[;, q]))^{\frac{1}{p_i}} \right)^{\frac{1}{\sum_{i \le n} \frac{1}{p_i}}}$

5: **for** $i=1:n$ **do**

6: Set $\tilde{X}_i[;, q] = \left( \frac{\beta_q}{p_i \mu_i g_i(X_i[;, q])} \right)^{1/p_i} \tilde{X}_i[;, q]$

7: **end for**

8: **end if**

9: **return** {$X_i[;, q]}_{i \le n}$

As noted earlier, a global solution {$X_i^*}_{i \le n}$ of HRSI must adhere to the balancing condition. Consequently, the solution space of HRSI can be restricted to balanced solutions only. The core insight of our work is that this restriction transforms the columnwise sum of regularization terms into a columnwise product of regularization.

Since all regularization terms evaluated at a balanced solution are equal up to factors $p_i$, the regularization terms $\sum_{q=1}^r \sum_{i=1}^n \mu_i g_i(X_i[;, q])$ after optimal balancing becomes $\sum_{i=1}^n \frac{1}{p_i} \sum_{q=1}^r \beta_q$. This expression can be rewritten as an implicit regularization term

$$ \tilde{\mu} \sum_{q=1}^{r} \left( \prod_{i=1}^{n} g_i(X_i[;, q])^{\frac{1}{p_i}} \right)^{\frac{1}{\sum_{i=1}^{n} \frac{1}{p_i}}} \quad (11) $$

where the penalty weight $\tilde{\mu}$ is defined as $\tilde{\mu} = \left( \prod_{i=1}^n (p_i \mu_i)^{\frac{1}{p_i}} \right)^{\frac{1}{\sum_{i=1}^n \frac{1}{p_i}}} \left( \sum_{i=1}^n \frac{1}{p_i} \right)$. When all $p_i$-s are equal, the implicit regularization coefficient conveniently simplifies to $\tilde{\mu} = n \left( \prod_{i=1}^n \mu_i \right)^{\frac{1}{n}}$, and the regularization term simplify to $\tilde{\mu} \sum_{q=1}^r \left( \prod_{i=1}^n g_i(X_i[;, q]) \right)^{\frac{1}{n}}$.

The implicit regularization (11) is of significant interest. The optimization obtained by replacing the explicit regularization with the implicit one,

$$ \underset{\forall i \le n, X_i \in \mathbb{R}^{m_i \times r}}{\operatorname{argmin}} f(\{X_i\}_{i \le n}) + \tilde{\mu} \sum_{q=1}^{r} \left( \prod_{i=1}^{n} g_i(X_i[;, q])^{\frac{1}{p_i}} \right)^{\frac{1}{\sum_{i=1}^{n} \frac{1}{p_i}}} \quad (12) $$

is fully scale-invariant in the sense of assumption A1. We can show that the explicit and implicit regularized problems are in fact essentially equivalent.

**Proposition 3.** *Given strictly positive regularization hyperparameters $\lambda_i$, the HRSI problem (1) is essentially equivalent to the implicit HRSI problem (12).*

The proof of this proposition relies on the fact that, denoting $\phi(\{X_i\}_{i \le n}) = f(\{X_i\}_{i \le n}) + \sum_{i=1}^n \mu_i \sum_{q=1}^r g_i(X_i[;, q])$,

$$ \min_{\forall i \le n, X_i \in \mathbb{R}^{m_i \times r}} \phi(\{X_i\}_{i \le n}) = \min_{\forall i \le n, X_i \in \mathbb{R}^{m_i \times r}, \prod_{i \le n} \Lambda_i = 1} \phi(\{X_i \Lambda_i\}_{i \le n}). \quad (13) $$

This shows that there is a trivial mapping between the minimizers of the explicit and implicit HRSI problems, where the solution to one can be transformed into the solution to the other through a scaling of the columns of the parameter matrices.

The implicit HRSI formulation provides an alternative approach to the original HRSI problem and could be leveraged in practice to design more efficient optimization algorithms. While we have not pursued this direction in this work, implicit HRSI has already served as the foundation for the development of an efficient algorithm for NMF, as proposed in recent work Marmin et al. [2023].

**Example 1: Ridge penalization applied to all parameter matrices induces low-rankness.** It has been established in the matrix low-rank approximation literature that applying squared $\ell_2$ norm penalizations to the Burer-Montero representation of the low-rank approximation matrix induces implicit low-rank regularization Srebro et al. [2004].
Consider the matrix LRA problem with ridge penalizations, using notations from Equation (4):

$$ \min_{X_1 \in \mathbb{R}^{m_1 \times r}, X_2 \in \mathbb{R}^{m_2 \times r}} \|M - X_1 X_2^T\|_F^2 + \mu (\|X_1\|_F^2 + \|X_2\|_F^2). \quad (14) $$

Proposition 3 shows that this problem is essentially equivalent to

$$ \min_{X_1 \in \mathbb{R}^{m_1 \times r}, X_2 \in \mathbb{R}^{m_2 \times r}} \|M - X_1 X_2^T\|_F^2 + 2\mu \sum_{q=1}^{r} \|X_1[:q]\|_2 \|X_2[:q]\|_2. \quad (15) $$

Noticing that the product of the $\ell_2$ norms of two vectors is the $\ell_2$ norm of their tensor product, we define $L_q = X_1$[:, $q] \otimes X_2$[:, $q]$ where $a \otimes b = ab^T$ is the outer product. The problem can thus be rewritten as

$$ \min_{L_q \in \mathbb{R}^{m_1 \times m_2}, \mathrm{rank}(L_q) \le 1} \|M - \sum_{q=1}^{r} L_q\|_F^2 + 2\mu \sum_{q=1}^{r} \|L_q\|_2. \quad (16) $$

This is a group LASSO regularized low-rank approximation problem. Group LASSO regularization is known to promote sparsity in groups of variables Yuan and Lin [2006], in this case, the rank-one terms $L_q$. The expected behavior of this regularization is semi-automatic rank selection, driven by the value of $\mu$.

These derivations are not specific to matrix problems. The squared $\ell_2$ norm penalties also lead to a group norm implicit regularization on the rank-one terms of tensor low-rank approximations. This is illustrated in the case of rNCPD in Section 4, where ridge regularization induces low nonnegative CP-rank approximations. However, since the group lasso regularization applies atop the rank-one constraints, we were unable to extend existing guarantees on the sparsity level of the solution. Additionally, to make the nuclear norm apparent in the matrix LRA problem and recover the result of Srebro and co-authors Srebro et al. [2004], we would need to rely on the rotation ambiguity of the $\ell_2$ norm, which falls outside of the framework proposed here. In other words, whether there are scenarios where the (implicit) group lasso penalty does not induce low-rank solutions remains an open problem.

**Example 2: Sparse $\ell_1$ penalizations may yield sparse and group sparse low-rank approximations, but individual sparsity levels cannot be tuned** In certain applications of regularized LRA, promoting sparsity across all parameter matrices $\{X_i\}_{i \le n}$ is desirable. For example, in the context of NMF with $n = 2$, a common formulation to achieve this is:

$$ \min_{X_1 \in \mathbb{R}_+^{m_1 \times r}, X_2 \in \mathbb{R}_+^{m_2 \times r}} \|M - X_1 X_2^T\|_F^2 + \mu_1 \|X_1\|_1 + \mu_2 \|X_2\|_1. \quad (17) $$

This model has been studied with a perspective similar to ours Papalexakis and Bro [2013], hence the observations we make below are known for sNMF.

One might intuitively expect that $\mu_1$ and $\mu_2$ control the sparsity levels of $X_1$ and $X_2$, respectively. However, this expectation must be tempered. Solutions to (17) are equal, up to scaling, to minimizers of the implicit problem

$$ \min_{L_q \in \mathbb{R}_+^{m_1 \times m_2}, \mathrm{rank}(L_q) \le 1} \|M - \sum_{q=1}^{r} L_q\|_F^2 + 2\sqrt{\mu_1 \mu_2} \sum_{q=1}^{r} \sqrt{\|L_q\|_1}. \quad (18) $$

We may observe that only the product of parameters $\mu_1$ and $\mu_2$ is significant at the implicit level. Consequently, one cannot independently tune $\mu_1$ or $\mu_2$ to adjust the sparsity levels of $X_1$ and $X_2$ individually. While the intuition that stronger regularization leads to sparser solutions holds, it applies only to the product $\mu_1\mu_2$, affecting the sparsity level of rank-one components $L_q$, as confirmed experimentally in Section 4. Further details on hyperparameter design are provided in Section 2.3.

Another observation is that the implicit regularization term in (18) is a group sparsity regularization with a mixed $\ell_{0.5,1}^0$ norm. While this regularization is unconventional, it intuitively promotes both sparsity in the parameter matrices and group sparsity in the rank-one components, i.e. perform automatic rank detection. However, these two effects are not controlled independently. This interdependence creates ambiguity in the resulting sparsity patterns: it is unclear whether individual sparsity within parameter matrices or group sparsity across rank-one components will dominate.

**What happens with $\ell_1 + \ell_2$ regularizations** Another commonly encountered pair of regularization functions in the context of regularized LRA is the $\ell_1$ and the squared $\ell_2$ norms, penalizing respectively $X_1$ and $X_2$. In this case Proposition 3 shows that the implicit regularization is a group norm on the rows or columns of rank-one matrices $L_q$, therefore inducing sparsity on the implicit problem as intended in the explicit formulation.
## 2.3 Choosing the hyperparameters

For a given set of hyperparameters {$\mu_i$}$_{i\le n}$, proportional scaling of all hyperparameters does not alter the minimum of the cost function, although it impacts the specific column scaling of the optimal parameter matrices. Thus it appears that hyperparameters {$\mu_i$}$_{i\le n}$ may be chosen equal: $\mu_i =: \mu$. This choice simplifies the problem formulation and ensures that, at optimality, the equality $p_ig_i(X_i[,q]) = p_jg_j(X_j[,q])$ holds for all $i < j \le n$ and for all $q \le r$. Consequently, the solutions are balanced independently of the hyperparameter value $\mu$. When the functions $g_i$ are identical (e.g., all $\ell_1$ norm), this choice results in uniform penalization and scaling of all parameter matrices. Furthermore, as $\mu$ goes to zero, all regularization terms vanish. This is practically desirable because the limit case $\mu = 0$ is consistent with the unregularized problem. Conversely, fixing a single hyperparameter $\mu_j$ while allowing the others to approach zero renders the problem ill-posed, as demonstrated in Section 2.1. When the dimensions $m_1$ and $m_2$ are vastly different, it may be of interest to choose $\mu_i = \frac{1}{m_i}\mu$ to scale the parameter matrices and avoid large discrepancies in their entries numerical values.

## 2.4 Related literature

The implicit regularization effects of explicit regularization in scale-invariant models have been previously recognized in the literature. An early contribution by Benichoux et al. Benichoux et al. [2013] highlighted the pitfalls of employing $\ell_1$ and $\ell_2$ regularizations in blind deconvolution. In particular, they demonstrated that solutions to such regularized model can be degenerate and non-informative. Their analysis is however tailored to convolutive mixtures, and similar degeneracies were not observed in the context of low-rank approximations.

Another closely related work studies the effect of scale-invariance for sNMF Papalexakis and Bro [2013]. Their study derived the same implicit problem as Equation (18) and identified the pitfall of independently tuning hyperparameters. They also notice empirically the existence of "scaling swamps" which we analyse in Section 3.3. However, their work is confined to sNMF, lacks a discussion of implicit regularization effects, and does not generalize the findings to broader classes of regularized low-rank approximations.

Further contributions within the sNMF framework include the work of Marmin et. al. that established an equivalence between sNMF with normalization constraints and a scale-invariant problem Marmin et al. [2023]. This work then proceeds to minimize the scale-invariant loss. Tan and Févotte, by generalizing the work of Mørup Mørup and Hansen [2009], have studied automatic rank selection using $\ell_1$ or $\ell_2$ penalizations, where the scales are controlled via the variance of the parameter matrices Tan and Févotte [2013]. Srebro et. al. demonstrated that ridge penalization in two-factor low-rank matrix factorization is equivalent to nuclear-norm minimization, a foundational result that initiated broader investigations into implicit regularization effects Srebro et al. [2004]. Despite their significance, these works are confined to specific models and do not capture the systemic influence of scale-invariance on regularized LRA or address the potentially misleading intuitions associated with explicit regularization.

Beyond sNMF, several works have examined implicit regularization in regularized low-rank factorization problems. Uschmajew observed that ridge regularization balances solutions in canonical polyadic decomposition Uschmajew [2012], while Roald et. al. showed that regularizations hyperparameters can be collectively tuned in regularized low-rank approximations with homogeneous regularizations Roald et al. [2022]. Parallel insights have emerged in the context of deep learning, where scale-invariance also plays a significant role. Ridge penalization on the linear layers of a feedforward ReLU network is equivalent to a group sparse LASSO problem Neyshabur et al. [2015], as recalled in a recent note by Tibshirani Tibshirani [2021]. Practical implications of explicit and implicit regularizations in deep neural networks have also been explored Ergen and Pilanci [2021], Stock et al. [2019], Gribonval et al. [2023], see this recent overview for further references Parhi and Nowak [2023].

Finally, it is worth noting that Proposition 2 is reminiscent of the quadratic variational formulation of norms Bach et al. [2011] albeit adapted to the multifactor setting.

# 3 Algorithms

In this section, we focus on developing practical optimization strategies tailored to regularized low-rank approximations with elementwise $\beta$-divergence loss functions. While the HRSI model is more general, a broad derivation of optimization strategies for it would be impractical.

More precisely, we consider nonnegative decompositions of input tensors or matrices, using the elementwise $\beta$-divergence Févotte et al. [2009] as the loss function $f$. The regularization terms are chosen from a family of functions that can be upper-bounded by separable quadratics, which includes $\ell_p^\beta$ norms. Let us recall that the $\beta$-divergence
between two matrices A and B is defined as

$$D_{\beta}(A, B) = \sum_{i,j} d_{\beta}(A[i, j], B[i, j]),$$

where, for scalars x and y,

$$d_{\beta}(x, y) = \begin{cases} \frac{x}{y} - \log \frac{x}{y} - 1 & \text{for } \beta = 0, \\ x \log \frac{x}{y} - x + y & \text{for } \beta = 1, \\ \frac{1}{\beta(\beta-1)} (x^{\beta} + (\beta-1)y^{\beta} - \beta xy^{\beta-1}) & \text{for } \beta \neq 0, 1, \end{cases} \quad (19)$$

see Basu et al. [1998], Eguchi and Kano [2001]. When $\beta = 2$, the $\beta$-divergence corresponds to the least-squares measurement. For $\beta = 1$ and $\beta = 0$, it corresponds to the Kullback-Leibler (KL) divergence and the Itakura-Saito (IS) divergence, respectively.<fnref n="6" /> The choice of $\beta$ is usually driven by two main motivations:

1. **Degree of homogeneity:** The $\beta$-divergence $d_{\beta}(x, y)$ is homogeneous of degree $\beta$, meaning that $d_{\beta}(\lambda x, \lambda y) = \lambda^{\beta}(x, y)$. This implies that factorizations obtained with $\beta > 0$ (such as the Euclidean distance or the KL divergence) will give more weight to the largest data values and less precision is expected in estimating of the low-power components. The IS divergence ($\beta = 0$) is scale-invariant, i.e., $d_{IS}(\lambda x|\lambda y) = d_{IS}(x|y)$. This is the only divergence in the $\beta$-divergences family with this property, meaning that low-power entries are treated as equally important as large entries in the divergence computation. This property is useful in applications such as audio source separation, where low-power frequency bands contribute perceptually as much as high-power frequency bands. For further details, see Lefèvre [2012] and Gillis [2020].

2. **Noise statistics:** The choice of $\beta$ is also influenced by the assumed noise statistics in the data. Specifically, there is an equivalence between the choice of the parameter $\beta$ and the maximum likelihood estimators of factors in the decomposition, given specific statistical distributions for the noise. In particular, the values $\beta = 2$, $\beta = 1$, and $\beta = 0$ correspond to the assumptions of Gaussian, Poisson, and exponential noise distributions, respectively. For further details, the reader is referred to Gillis [2020] and references therein.

Our approach is based on three main components: Cyclic Block-coordinate Descent (BCD), Majorization-Mimization (MM), and Optimal balancing. The latter consists in performing a balancing step for all blocks to minimize the regularization terms, as described in Section 2.2.2. In the following, we discuss these three components in more detail. We then introduce our general algorithm, referred to as the *Meta-Algorithm*, and conclude by presenting its convergence guarantees.

## 3.1 BCD: which splitting ?

We begin by splitting the variables into disjoint blocks<fnref n="7" /> and sequentially minimize the objective function over each block. The updates for each block will be similar, and the objective function in (1) is separable columnwise. Therefore to simplify the notations, we consider the update of the block of variables $X_i$[;, q] (one column of $X_i$, dropping the i index), denoted $x \in \mathbb{R}_+^m$, while the other blocks of variables are fixed. Recall that $x[k]$ represents the entry in vector x at index k. We are interested in solving the following subproblem:

$$\min_{x \in \mathbb{R}_+^m} f(x) + \mu g(x). \quad (20)$$

## 3.2 Majorization-Minimization (MM)

To address Problem (20), we employ the widely-used MM framework. In this approach, the objective function $\phi$ is iteratively globally majorized by a locally tight approximation $\tilde{\phi}$, which can be efficiently minimized. In the context of regularized LRA, such majorizers are typically designed to be separable, enabling cost-efficient update rules. The originality in our approach lies in the presence of the regularizer $g$.

We now detail how to derive a separable majorizer $\bar{\phi}(x|\tilde{x})$ tight at a given point $\tilde{x}$ for the objective function $\phi(x) := f(x) + \mu g(x)$. The first term $f(x)$ is majorized following the methodology introduced in Févotte and Idier [2011], where the beta-divergence is decomposed as the sum of a convex and a concave function. The convex part is majorized using Jensen's inequality, while the concave part is upper-bounded by its first-order Taylor approximation. The resulting separable majorizer for $f(x)$ built at the current iterate $\tilde{x}$ is denoted $\tilde{f}(x|\tilde{x})$. For the second term, we make the following assumptions:

<footnote n="6">⁶Note that, with convention $a \log 0 = -\infty$ for $a > 0$ and $0 \log 0 = 0$, the KL-divergence is well-defined.</footnote>

<footnote n="7">⁷their union is the initial and full set of variables</footnote>
**Assumption 1.** Leplat et al. [2021] For the smooth and positive-definite regularization function $g: x \in \mathbb{R}_+^m \to \mathbb{R}_+$, and any $\tilde{x} \in \mathbb{R}_+^m$, there exists nonnegative constants $C_k$ with $1 \le k \le m$ such that the inequality

$$g(x) \le \bar{g}(x|\tilde{x}) := g(\tilde{x}) + \langle \nabla g(\tilde{x}), x - \tilde{x} \rangle + \sum_{k=1}^{m} \frac{C_k}{2} (x[k] - \tilde{x}[k])^2 \quad (21)$$

holds for all $x \in \mathbb{R}_+^m$.

In other words, the regularization function can be upper-bound by a tight separable majorizer given in Equation (21). This property holds for all $\ell_p^p$ norms employed in this work as regularizers. Consequently, a separable majorizer $\bar{\phi}(x|\tilde{x})$ for $\phi(x)$ can be constructed at any $\tilde{x} \in \mathcal{X}$ by combining $\tilde{f}(x|\tilde{x})$ and $\bar{g}(x|\tilde{x})$. This majorizer is then minimized to compute the update $\hat{x}$. The existence and uniqueness of the minimizer $\hat{x}$ can be established under mild conditions Leplat et al. [2021].

**Computation of the update rules for beta divergence and $\ell_p^p$ regularizations** The MM updates for $\beta$ in $[0, 2]$ and regularization functions that are $\ell_p^p$ norms, both of which are important particular cases, are derived in Appendix F and Appendix G. One of our contributions is to derive update rules for $\beta = 1$ and $\beta = 0$ with squared $\ell_2$ regularizations leading to convergence guarantees for the iterates. These cases have not been addressed in the existing literature. As an illustrative example, for the specific pair $\beta = 1$ and $p = 2$, the update for each entry of the block of variable under consideration, denoted $x$, has the general form

$$x[k] \leftarrow \frac{\sqrt{c^2 + 8\mu x[k]t} - c}{4\mu}, \quad (22)$$

where $\mu > 0$ and $c$ and $t$ are nonnegative real-valued functions of the input data and the other blocks of variables held fixed.

More generally, a closed form expression for the update $\hat{x}$ can be derived if $\beta$ belongs to $\{0, 1, 3/2, 2\}$ according to Leplat et al. [2021]. Outside this set of values for $\beta$, an iterative scheme such as Newton-Raphson is required to compute the minimizer. Note that while MM updates for $\beta = 2$ can be derived within the proposed framework, in this case, the literature suggests the superiority of an alternating least squares approach called HALS Gillis and Glineur [2012], Marmoret et al. [2020].

## 3.3 Balancing improves convergence speed

One of the contributions of this work is the introduction of a numerical balancing of the columns of low-rank components, building upon the results presented in Section 2. While the balancing procedure has been established, its practical utility remains to be clarified. To this end, we present an example where the existence of “scaling swamps” Papalexakis and Bro [2013] for the Alternating Least Squares (ALS) algorithm can be proven.

For a given real value $y$ and a regularization parameter $0 < \lambda \ll y$, consider the scalar optimization problem

$$\min_{x_1, x_2} f(x_1, x_2) \quad \text{where} \quad f(x_1, x_2) = (y - x_1 x_2)^2 + \lambda(x_1^2 + x_2^2). \quad (23)$$

Using the results from Section 2, we can immediately infer that the minimizers $x_1^*$ and $x_2^*$ must be balanced. This yields optimal solutions $x_1^* = x_2^* = \sqrt{y - \lambda}$ (for $\lambda \le y$).

Setting aside the fact that we know the optimal solution, a classical iterative BCD method from the regularized LRA literature for solving this problem would be the ALS algorithm. Denoting the iterates at iteration $k$ as $x_1^{(k)}$ and $x_2^{(k)}$, the update rules are straightforwardly derived as

$$x_1^{(k+1)} = \frac{x_2^{(k)} y}{x_2^{(k)} + \lambda} \quad \text{and} \quad x_2^{(k+1)} = \frac{x_1^{(k+1)} y}{x_1^{(k+1)} + \lambda}. \quad (24)$$

We show in Appendix C that when $\lambda$ is small relative to $y$, and for sufficiently large $k$,

$$\frac{x_1^{(k+1)} - \sqrt{y - \lambda}}{x_1^{(k)} - \sqrt{y - \lambda}} \approx 1 - 4 \frac{\lambda}{y} \quad (25)$$

with a similar expression holding for $x_2$. Thus, ALS converges almost sub-linearly when $x_1$ and $x_2$ are near the optimum, particularly when $\frac{\lambda}{y} \ll 1$. The smaller the ratio $\frac{\lambda}{y}$, the slower the convergence. A numerical simulation of this toy problem solved with the ALS algorithm experimentally confirmed these observations, see Figure 1.

Similar slow convergence behaviors are observed in more realistic regularized LRA problems, as shown in Section 4. To address this, the Meta-Algorithm below includes a step to optimally balance the iterates after each outer iteration.
<figure><img src="image_5.png" /><figcaption>Figure 1: Study of the convergence of ALS for the toy problem (23). We choose $y = 10$ and $\lambda = 10^{-3}$. The number of iterations of ALS is set to 20000. Top left: cost function $f$ with respect to iteration index. Top right: values of $x_1$ (orange), $x_2$ (blue) and $\sqrt{y-\lambda}$ (magenta). Bottom left: value of $x_1x_2$ (blue) and $y-\lambda$ (magenta). Bottom right: empirical value of the cost function variation (blue) and theoretical value $16\frac{\lambda^2}{y}e_k^2$ (magenta). Several observations can be made: The value of $x_1x_2$ rapidly converges towards $y-\lambda$ (bottom left graph). However the cost converged after at least 5000 iterations, and the individual values of $x_1$ and $x_2$ are slow to converge to the optimum $\sqrt{y-\lambda}$ with respect to how fast $x_1x_2$ converged to $y-\lambda$. We can observe a linear convergence rate with a multiplicative constant very close to 1 in the bottom right graph, and a close match between the theoretical approximate cost decrease and the practical observation of that decrease.</figcaption></figure>

## 3.4 Meta-Algorithm

The proposed method for solving problem (1) is summarized in Algorithm 2. The algorithm is iterative, alternating, and consists of two main steps within each block update. First the block $\{X_i\}_{i \le n}$ is updated using MM. Second, optimal balancing of the block of variables is performed. This balancing step is computationally inexpensive and necessarily decreases the objective function, see Proposition 4. Thus, it does not compromise the convergence guarantee offered by the MM-based updates, ensuring the convergence of the sequence of objective function values. However, the convergence guarantees for the sequence of iterates will depend on the choice of the value $\beta$ in Problem (77). These points are discussed in more detail in Section 3.6 below.

**Initialization** In Algorithm 2, we propose an initialization stage where the blocks of variables are scaled (lines 1-2). This scaling is essential to prevent the "zero-locking" phenomenon. To illustrate this issue, let us consider a matrix case $n = 2$. If the regularization functions $g_1$ and $g_2$ are chosen as $\ell_1$ norms, the updates for the columns of parameter matrix $X_1$ reduce to solving a series of LASSO problems, with the mixing matrix given by $X_2^T$. For such problems, it is well known that if $\|MX_2\|_\infty \le \mu_1$, the solution is necessarily $X_1 = 0$. In this scenario, $(0, 0)$ being a saddle point it is very likely that Algorithm 2 will return $X_1 = X_2 = 0$. The scaling step mitigates this issue by ensuring that the initial values of $X_1$ and $X_2$ are not excessively small.
**Algorithm 2 Meta-Algorithm**

**Require:** An input tensor $\mathcal{T} \in \mathbb{R}^{m_1 \times \cdots \times m_n}$, an initialization $X_i^{(0)}$ and regularization hyperparameters $\mu_i > 0$, the factorization rank $r$, a maximum number of iterations, maxiter.

1: $\eta \in \text{Argmin}_{\eta \ge 0} f(\{\eta X_i^{(0)}\}_{i \le n}) + \sum_{i=1}^n \mu_i \sum_{q=1}^r \eta^{p_i} g_i(X_i^{(0)}[:; q])$ % Rescaling

2: $X_i^{(1)} \leftarrow \eta X_i^{(0)}$ for all $i = 1, \dots, n$.

3: % Main optimization loop:

4: **for** $k = 1$ : maxiter **do**

5: &nbsp;&nbsp;&nbsp;&nbsp;**for** $i = 1$ : $n$ and $q = 1$ : $r$ **do**

6: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$X_i[;, q]^{(k+1)} \leftarrow \text{argmin}_{x \in \mathcal{X}_i} \bar{\phi}(x|X_i[;, q]^{(k)})$ % Update of parameter matrices

7: &nbsp;&nbsp;&nbsp;&nbsp;**end for**

8: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$\{X_i^{(k+1)}[;, q]\}_{i \le n}$ balanced using Algorithm 1 for all $q \le r$ % Optimal balancing

9: **end for**

10: **return** $\{X_i[;, q]\}_{i \le n}$

## 3.5 Convergence guarantees

We establish two convergence results for the proposed Meta-Algorithm: (1) the convergence of the cost function, and (2) under further assumptions, the convergence of the iterates to a stationary point. While these results are not unexpected, they are not straightforward corollaries of existing literature. Detailed proofs, technical details, and precise formulations of the assumptions are provided in Appendix D.

Before we give the main convergence results, it is convenient to write down a general formulation for our optimization problems. Let $z = (\{X_i[;, q]\}_{i \le n, q \le r})$ denote the concatenation of the variables $X_i[;, q]$ and let $\phi$ represent the objective function from Problem (1). We consider the following class of multiblock optimization problems:

$$ \min_{z_i \in \mathcal{Z}_i} \phi(z_1, \dots, z_s), \qquad (26) $$

where $z$ is decomposed into $s$ blocks $z = (z_1, \dots, z_s)$, $\mathcal{Z}_i \subseteq \mathbb{R}^{m_i}$ is a closed convex set, and $\mathcal{Z} = \mathcal{Z}_1 \times \dots \mathcal{Z}_s \subseteq \text{dom}(\phi)$. The function $\phi: \mathbb{R}^m \to \mathbb{R} \cup \{+\infty\}$, with $m = \sum_{i=1}^s m_i$, is continuously differentiable over its domain, lower bounded and has bounded level sets.

**Theorem 1.** Let $z^{(k)} = (\{X_i^{(k)}[;, q]\}_{i \le n, q \le r}) \ge 0$, and let $\bar{\phi}_i$ be block-wise tight majorizers<fnref n="8" /> for $\phi: \mathcal{Z} \to \mathbb{R}_+$ from Problem (26), with $\mathcal{Z} := \mathcal{Z}_1 \times \dots \mathcal{Z}_s \subseteq \text{dom}(\phi)$ and $1 \le i \le nr$. Then $\phi$ is monotone non-increasing under the updates of Algorithm 2, and the sequence of cost values converges.

**Theorem 2.** Let block-wise tight majorizers $\bar{\phi}_i(y_i, z)$ be quasi-convex in $y_i$ for $i = 1, \dots, nr$. Additionally, assume that the directional derivatives of the majorizers and the objective coincide locally for each block and that the majorizer is continuous. Further, assume that the MM updates and the optimal scaling have unique solutions. Under these conditions, every limit point $z^\infty$ of the iterates $\{z^{(k)}\}_{k \in \mathbb{N}}$ generated by Algorithm 2 is a coordinate-wise minimum of Problem (26). Furthermore, if $\phi$ is regular at $z^\infty$, then $z^\infty$ is a stationary point of Problem (26).

The proofs of these theorems are adapted from the convergence results to stationary points obtained in [Razaviyayn et al., 2013, Theorem 2] for the Block Successive Upperbound Minimization (BSUM) framework. They leverage the fact that the balancing procedure always decreases the cost $\phi$.

All the hypotheses of Theorem 1 and Theorem 2, except for directional differentiability (see Section 3.6), are easily satisfied by many regularized LRA models, in particular those discussed in Section 4.

## 3.6 Handling near-zero entries with the $\beta$-divergence

To ensure the convergence of the iterates generated by Algorithm 2 to critical points of Problem (26), one of the required assumptions is the directional differentiability of the function $\phi$, see Appendix D. However, this property does not hold at zero as soon as $\beta < 2$, see for instance the discussions in Hien and Gillis [2021]. To guarantee convergence for the iterates, a common solution is to impose constraints $X_i \ge \epsilon$ for all $i \in [1, n]$, where $\epsilon > 0$ Hien and Gillis [2021]. Under these constraints, the updates of parameter matrices in Algorithm 2 (Step 6) become

$$ X_i[;, q]^{(k+1)} \leftarrow \max_x \left( \epsilon, \operatorname*{argmin}_x \bar{\phi}(x|X_i[;, q]^{(k)}) \right), \qquad (27) $$

<footnote n="8">⁸See Definition 1</footnote>
see Takahashi and Ryota [2014] where this approach has been successfully employed to derive multiplicative updates algorithms for NMF with global convergence guarantees.

Given the modified updates, ensuring the convergence guarantees stated in Theorem 2 is not straightforward. The reason lies in the optimal balancing step of Algorithm 2. Indeed, if one entry of $X_i[\cdot, \cdot]$ reached $\epsilon$, the balancing procedure may decrease this entry below $\epsilon$, making it non-feasible. We discuss two strategies to address this issue:

* **Suboptimal scaling with projection** If a column of a parameter matrix becomes elementwise smaller than $\epsilon$, all the corresponding columns in other parameter matrices are set to $\epsilon$ as well. When at least one entry in each column is strictly larger than $\epsilon$, we first set all values smaller than $\epsilon$ to zero, optimally scale the resulting projected columns, and then project back the outcome such that all entries are larger than $\epsilon$. This approach leads to a suboptimal scaling and may therefore break convergence, but we anticipate that these theoretical considerations will not affect the algorithm behavior noticeably. This is the recommended approach in practice, used in Section 4.

* **Stop balancing to guarantee convergence** A second approach is to stop balancing as soon as any entry in the parameter matrices is smaller than $\epsilon$. This strategy ensures the convergence of the iterates to critical points of the $\epsilon$-perturbed Problem (26) through BSUM, with the balancing steps primarily serving to accelerate the convergence of the algorithms during the initial iterations. Balancing may even only be applied upon initialization. This approach appears to be the most suitable from a theoretical perspective, but in practice, the balancing may not be performed for many iterations which is undesirable. Indeed, anticipating on Section 4, performing optimal balancing along all iterations is practically profitable.

Finally, for the case $\beta \ge 2$, the $\beta$-divergence exhibits continuous differentiability over the feasible set, thereby ensuring the convergence of Algorithm 2 as outlined in Theorem 2. However, thresholding the entries in parameter matrices $X_i$ to $\epsilon$ when $\beta \ge 2$ may still be of practical interest. This is because, with multiplicative updates, any entry in a parameter matrix that reaches zero cannot be further modified by the algorithm. This can prevent convergence in practice, especially with poor initialization or insufficient machine precision. In the proposed implementations of the meta-algorithm, whenever nonnegativity is enforced, a small $\epsilon$ is used to lower-bound the entries of all parameter matrices.

# 4 Showcases

In this section, we study three specific instances of HRSI: sNMF, rNCPD, and sNTD. We begin by relating each model to HRSI, outlining the balancing update and the implicit equivalent formulation, and detailing the implementation of the proposed alternating MM algorithm. Next, we evaluate the performance of the proposed method with and without balancing.

Additionally, we provide an implementation of the proposed algorithm for each showcase problem, along with reproducible experiments, in a dedicated repository<fnref n="9" />. For multilinear algebra, we rely on Tensorly Kossaifi et al. [2019], and we use the package shootout<fnref n="10" /> for running the experiments.

## 4.1 Explicit formulations and relation to HRSI

We consider sNMF, rNCPD, and sNTD respectively, explicitly formulated as follows:

$$ \min_{X_1 \ge 0, X_2 \ge 0} \mathrm{KL}(M, X_1 X_2^T) + \mu_1 \|X_1\|_1 + \mu_2 \|X_2\|_1, \qquad (\mathrm{sNMF}) $$

$$ \argmin_{X_i \ge 0} \|T - I_r \times_1 X_1 \times_2 X_2 \times_3 X_3\|_F^2 + \mu (\|X_1\|_F^2 + \|X_2\|_F^2 + \|X_3\|_F^2), \qquad (\mathrm{rNCPD}) $$

$$ \argmin_{X_i \ge 0, G \ge 0} \mathrm{KL}(T, G \times_1 X_1 \times_2 X_2 \times_3 X_3) + \mu (\|G\|_1 + \|X_1\|_F^2 + \|X_2\|_F^2 + \|X_3\|_F^2), \qquad (\mathrm{sNTD}) $$

where sNTD has sparsity enforced on the core tensor $G$. For sNMF, two regularization hyperparameters $\mu_1 \ne \mu_2$ are introduced to test the observations made in Section 2.2.

As discussed in Section 1.3, sNMF, and similarly rNCPD, are instances of HRSI. sNTD may also be formulated as a rank-one HRSI model, by considering a vectorized formulation:

$$ \argmin_{x_1 \ge 0, x_2 \ge 0, x_3 \ge 0, g \ge 0} \mathrm{KL}(T, \{x_1, x_2, x_3, g\}) + \mu (\|g\|_1 + \|x_1\|_F^2 + \|x_2\|_F^2 + \|x_3\|_F^2), \qquad (28) $$

<footnote n="9">⁹https://github.com/vleplat/NTD-Algorithms</footnote>

<footnote n="10">¹⁰https://github.com/cohenjer/shootout</footnote>
where $x_i = \text{vec}(X_i)$ and $g = \text{vec}(\mathcal{G})$ for a chosen arbitrary vectorization procedure, and where KL is the KL divergence against tensor $\mathcal{T}$ with vectorized inputs.

### 4.1.1 Optimal balancing and implicit formulations

We can apply the results of Section 2 straightforwardly for each model to obtain the following balancing updates for each column indexed by $q \le r$ and each parameter matrix $X_i$:

$$X_i[\cdot, q] \leftarrow \frac{\sqrt{\mu_1 \mu_2 \|X_1[\cdot, q]\|_1 \|X_2[\cdot, q]\|_1}}{\mu_i \|X_i[\cdot, q]\|_1} X_i[\cdot, q], \qquad (\text{sNMFb})$$

$$X_i[\cdot, q] \leftarrow \left( \sqrt{2} \|X_i[\cdot, q]\|_2^{-\frac{1}{3}} \prod_{j \ne i} \|X_j[\cdot, q]\|_2^{\frac{2}{3}} \right) X_i[\cdot, q]. \qquad (\text{rNCPDb})$$

Balancing for sNTD is done with respect to the vectorized inputs and therefore there is virtually only one column to scale for each parameter matrix:

$$X_i \leftarrow \sqrt{\frac{2^{\frac{3}{2}} (\|\mathcal{G}\|_1 \prod_{j=1}^3 \|X_j\|_F)^{2/5}}{2 \|X_i\|_F^2}} X_i \quad \text{and} \quad \mathcal{G} \leftarrow \frac{2^{\frac{3}{2}} (\|\mathcal{G}\|_1 \prod_{j=1}^3 \|X_j\|_F)^{2/5}}{\|\mathcal{G}\|_1} \mathcal{G}. \qquad (\text{sNTDb})$$

Furthermore, the essentially equivalent implicit HRSI problem for sNMF and rNCPD can respectively be cast as

$$\min_{L_q \in \mathbb{R}_+^{m_1 \times m_2}, \text{rank}(L_q) \le 1} \text{KL}(M, \sum_{q=1}^r L_q) + 2\sqrt{\mu_1 \mu_2} \sum_{q=1}^r \sqrt{\|L_q\|_1}. \qquad (29)$$

$$\inf_{\{\mathcal{L}_q\}_{q \le r}, \text{rank}(\mathcal{L}_q) \le 1} \|\mathcal{T} - \sum_{q=1}^{r} \mathcal{L}_q\|_F^2 + 3\mu \sum_{q=1}^{r} \|\mathcal{L}_q\|_F^{\frac{2}{3}}, \qquad (30)$$

with $\mathcal{L}_q = X_1[\cdot, q] \otimes X_2[\cdot, q] \otimes X_3[\cdot, q]$. In a nutshell, sNMF in its implicit formulation is regularized with a mixed norm that promotes sparsity at both the entry and the component levels, while rNCPD is regularized with a group sparse norm on the rank-one components.

The sNTD model exhibits rotational invariances that our approach does not address, as they do not match the assumptions of the HRSI model. Based on formulation (sNTD), we can derive an implicit formulation; however, balancing is applied to all components simultaneously, which obscures component-wise group sparse regularization. In Appendix H, we demonstrate that a Sinkhorn-like procedure can be designed to balance sNTD algorithmically. However, this approach did not yield improvements over the simpler vectorization-based balancing and deserves further exploration.

### 4.1.2 Majorization-Minimization update rules

For sNMF, the MM updates take the form of “multiplicative” updates:

$$\hat{X}_1 = \max \left( X_1 \odot \frac{X_2 \frac{M}{X_2^T X_1}}{\mu_1 e_{m_1 \times r} + e_{m_1} \otimes X_2^T e_{m_2}}, \epsilon \right), \qquad (31)$$

where $e_n$ is a vector or matrix of ones of size $n$ and the maximum is taken elementwise. A similar update is obtained for $X_2$, detailed computations are found in Appendix F. For sNTD, consider the unfolded tensor along the $i$th mode $\mathcal{T}_{[i]} \in \mathbb{R}_+^{m_i \times \prod_{j \ne i} m_j}$ and the Kronecker products $U_i := \mathcal{G}_{[i]} (\otimes_{j \ne i} X_j)^T \in \mathbb{R}_+^{r_1 \times \prod_{j \ne i} m_j}$. The MM updates then read

$$\hat{X}_i = \max \left( \frac{[C_i^2 + S_i]^2}{4\mu} - C_i, \epsilon \right), \qquad (32)$$

$$\hat{\mathcal{G}} = \max \left( \mathcal{G} \odot \frac{\mathcal{N} \times_1 X_1^T \times_2 X_2^T \times_3 X_3^T}{\mu \mathcal{E}_{r_1 \times r_2 \times r_3} + \mathcal{E}_{m_1 \times m_2 \times m_3} \times_1 X_1^T \times_2 X_2^T \times_3 X_3^T}, \epsilon \right), \qquad (33)$$

where $C_i = EU_i^T$ with $E$ and $\mathcal{E}$ all-one matrices and tensors, $S_i = 8\mu X_i \odot (\frac{\mathcal{T}_{[i]}}{X_i U_i} U_i^T)$ and $\mathcal{N} = \frac{\mathcal{T}}{\mathcal{G} \times_1 X_1 \times_2 X_2 \times_3 X_3}$. The matrices and tensors of ones are not contracted as such in practice; instead, we use summations of columns in the parameter matrices.

Regarding rNCPD, we solve the problem using HALS Gillis and Glineur [2012] as explained in Section 3.
4.1.3 Synthetic experiments

**Experiments description** We study the impact of the balancing procedure in the meta-algorithm, compared to the same algorithm without explicit balancing. We observe the performance of both strategies in terms of loss value at the last iteration. In addition, we explore a third option in which balancing is performed only at initialization. Some experiments conducted here are adapted from our previous work Cohen [2023].

The process for generating synthetic data is as follows:

1. We construct the ground truth low-rank matrix/tensor $M$ with dimensions $m_i = 30$ using low-rank models sNMF, rCPD, or sNTD.

2. The rank $r$ of the model is fixed at 4, and the estimated rank $r_e$ for the approximation is also set to 4.

3. For the generation of noisy data:

• For sNMF and sNTD, the noisy data matrix/tensor is sampled from the Poisson distribution $P(\alpha M)$.

• For rNCPD, it is sampled from a Gaussian independent and identically distributed (i.i.d.) distribution with variance $\alpha$ and mean $M$.

4. The parameter $\alpha$ is selected to ensure that the average Signal-to-Noise Ratio (SNR) is 40 dB for sNMF and sNTD, and 200 dB for rNCPD.

5. Finally, the noisy data is normalized using its Frobenius norm.

The ground-truth parameter matrices $X_i$ are sampled elementwise from i.i.d. Uniform distributions over the interval $[0, 1]$. For sNMF, these matrices are sparsified using hard thresholding, resulting in 30% of their entries being set to zero. In the case of sNTD, the core tensor is sampled similarly and sparsified such that 70% of its entries are null.

The hyperparameters are selected as follows: We set $\epsilon = 10^{-16}$. The initial parameter matrices are drawn like the true ones; however, the first matrix $X_1$ is scaled by 100. This results in unbalanced parameter matrices, which serves to emphasize the effect of the balancing procedure. We optimally scale the initialization with respect to the $\ell_2$ loss, as described in Algorithm 2. The number of outer iterations is set to 500 for both sNMF and sNTD, while it is set to 50 for rNCPD. The number of inner iterations is set to 10. Each experiment is repeated 50 times. The rNCPD problem is addressed using the HALS algorithm, which is known for its fast convergence, thereby requiring fewer outer iterations.

We conduct several additional tests specific to each model. For sNMF, to validate that the individual values of $\mu_1$ and $\mu_2$ do not influence the relative sparsity levels of $X_1$ and $X_2$, we examine the ratio $\frac{\|X_1\|_0}{\|X_2\|_0}$ in relation to the regularization hyperparameters. An entry in a matrix is considered null if it is below the threshold $2\epsilon$. Moreover, we fix the regularization hyperparameter $\mu_1 = 1$, which allows us to study the relative sparsity of the parameter matrices concerning the individual values of the regularization hyperparameters. For rNCPD, to investigate the number of zero components, we overestimate the number of components for the approximation $r_e$ (estimated rank) to 6. We anticipate that the ridge penalization will function as a rank penalization, effectively pruning two components. A component is deemed null if the weight associated with it (the product of the $\ell_2$ norms of each parameter matrix along each column) is lower than $1000\epsilon$ (there are $30 \times 6 = 180$ entries). For sNTD, the first estimated rank $r_1$ is also set to 6 instead of 4, and we analyze both the number of nonzero first-mode slices and the overall sparsity in the core tensor $\mathcal{G}$.

**Experiments Results** Results are shown in Figures 2 to 4. We validate several observations from the theoretical derivations:

• For small values of $\mu$, as predicted in Section 3.3, the meta-algorithm without explicit balancing struggles to converge in terms of loss value. However, as shown in the supplementary material Figure 6, the parameter matrices are still well estimated without balancing. This indicates that while the unbalanced algorithm identifies good candidate solutions for small $\mu$, it requires many iterations to balance these estimates. Explicit balancing reduces the regularization terms, resulting in significantly lower costs after running the algorithm. Balancing the initial guess offers less improvement compared to balancing at each iteration, except for sNMF.

• The observations derived from the implicit formulations are validated numerically. The relative sparsity levels in the sNMF parameter matrices remain close to 1, regardless of the ratio $\frac{\mu_2}{\mu_1}$. For a wide range of regularization hyperparameter values, the estimated number of components in rNCPD matches the true count of four, indicating that ridge regularization effectively performs automatic rank selection. In sNTD, since all algorithms should converge to similar solutions, the unbalanced method produces solutions denser than expected. Thus, balancing aids in identifying zero entries. The decomposition rank for the first mode was set larger than the true rank, and the balanced algorithm more reliably prunes the extra components.
<figure><img src="image_6.png" /><figcaption>Figure 2: Results for the simulated experiment on sNMF. The left plot shows the loss after the stopping criterion is reached, plotted against the regularization parameter $\mu_2$ while $\mu_1 = 1$. The right plot displays the ratio of the sparsity of the estimated matrix $X_1$ to that of the estimated matrix $X_2$. Results are presented for three scenarios: balancing used only at initialization (blue), balancing applied at each iteration (red), and no balancing (green). Notably, the left plot indicates that balancing reduces the loss function value upon reaching the convergence criterion, particularly for the weakest regularizations. In the right plot, it is observed that tuning the regularization hyperparameter $\mu_2$ while fixing $\mu_1 = 1$ does not affect the sparsity ratio of the parameter matrices $X_2$ and $X_1$, confirming that the sparsity of matrix $X_2$ cannot be adjusted independently of matrix $X_1$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 3: Results for the simulated experiment on ridge Nonnegative CPD. The normalized loss after the stopping criterion is reached is shown with respect to the regularization parameter $\mu$ on the left plot, while the number of components within the estimated rNCPD model with respect to the regularization parameter is shown on the right plot. Results when balancing is used at initialization only, at each iteration, and not used are shown respectively in blue, red, and green. One may observe on the left plot that balancing helps reduce the loss function value when the convergence criterion is reached. This reduction is significant in particular for small regularization parameter $\mu$ and when balancing is performed at each iteration, although balancing only the initialization already improves over no balancing. On the right plot, one may observe that for a wide range of ridge regularization hyperparameter values, the number of estimated components in rNCPD matches true rank $r = 4$. Therefore ridge regularization has an implicit group sparse action on the rank-one components as predicted in theory. This phenomenon happens regardless of the balancing procedure.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 4: Results for the simulated experiment on sparse NTD. The normalized loss after 500 outer iterations is presented in the top left plot, showing its relationship with the regularization parameter $\mu$. The top right plot displays the scaled sparsity of the estimated core tensor against the regularization parameters, while the bottom plot illustrates the number of nonzero components along the first mode of the estimated core. Results are shown in blue for balancing used only at initialization, in red for balancing at each iteration (using scalar balancing), and in green for no balancing. The top left plot indicates that balancing effectively reduces the loss function value upon reaching the convergence criterion, with a notable reduction for small regularization parameter $\mu$, especially when balancing is applied at each iteration. Even balancing only at initialization improves results compared to no balancing. In the top right and bottom plots, $\ell_1$ regularization on the core tensor promotes sparsity and helps select the appropriate number of components, $r_1 = 4$. Balancing—whether at initialization or during iterations—broadens the range of regularization values for which components are effectively pruned.</figcaption></figure>

### 4.1.4 sNTD for Music Segmentation

We showcase sNTD with balancing in a music redundancy detection task. The use of NTD for analyzing music was proposed in Smith and Goto [2018] and further developed in Marmoret et al. [2020]. For an audio recording of a song segmented into bars, a tensor spectrogram is constructed by stacking time-frequency spectrograms computed from each bar of the audio signal. Each column of the time-frequency spectrogram matrix captures frequency information about the one-dimensional audio signal within a small time window. The resulting data forms a tensor spectrogram $\mathcal{T}$ of size $F \times T \times B$, where $F$ is the number of frequency bins, $T$ is the number of time bins, and $B$ is the number of bars in the song. By applying NTD to the tensor spectrogram, redundancies across the bars can be detected, allowing for the factorization of spectral and short-term temporal content into a collection of patterns. Notably, the redundancy information is contained in the last mode parameter matrix $X_3$, which can be utilized for downstream tasks such as automatic segmentation. Each column $X_3[\cdot, :q]$ of matrix $X_3$ corresponds to a specific pattern, defined by the product of the other parameters $X_1, X_2$, and the slice $G[\cdot, \cdot; q]$ of tensor $G$. Each row of matrix $X_3$ maps these patterns to the bar indices. Matrix $X_3$ visualizes the relationship between patterns and bars, with the expectation that bars containing similar content will utilize the same patterns. Consequently, matrix $X_3$ should exhibit sparsity, as only a limited number of patterns are employed for each bar.

We consider sNTD with sparsity imposed on the parameter matrix $X_3$ instead of $G$. The derivations for this variant are similar to those of sNTD with a sparse core tensor and are therefore omitted. The song selected for this experiment is "Come Together" by the Beatles, which has been used previously to demonstrate NTD's ability to detect patterns and redundancies in music. After processing the song as described in Marmoret et al. [2020], we obtain a tensor of size
$F = 80, T = 96$, and $B = 89$. The core dimensions are set to $32 \times 12 \times 10$, indicating that the song should be factorized into at most 10 patterns. The algorithm runs for 50 iterations, with 10 inner iterations per run. The regularization parameter is chosen from the grid $[0, 10, 100, 300, 500]$. Initialization is performed using HOSVD De Lathauwer and Vandewalle [2004]. While the choice of initialization method is important, we opted for the method proposed in Marmoret et al. [2020] without further exploration.

Figure 5 displays the loss function for two variants of the sNTD algorithm: one without balancing and one with scalar balancing, along with the estimated parameter matrices $X_3^T$ for each setup. The expert, ground truth segmentation of the song is indicated by green vertical lines. Notably, the balanced algorithms converge faster than the unbalanced variant. Additionally, the sparsity imposed on $X_3$ significantly affects the quality of redundancy detection. Ideally, the pattern indices used for each bar should align closely with the expert segmentation, resulting in group sparsity between the green vertical lines in the rows of the matrices shown in Figure 5. This alignment is particularly evident with a regularization parameter of three hundred. At this level, balancing visibly influences the final sparsity pattern of matrix $X_3$, although it is arguably difficult to assess the superiority of the balanced algorithm from visual inspection only.

# 5 Conclusions and perspectives

The design of efficient algorithms for computing low-rank approximations is an active research field, in particular in the presence of constraints and penalizations. In this work, we have shed light on the hidden effects of regularizations in nonnegative low-rank approximations, by studying a broader class of problems we term Homogeneous Regularized Scale-Invariant problems. We have shown that the scale-ambiguity inherent to these problems induces a balancing effect of penalizations imposed on several parameter matrices. This balancing effect has both practical and theoretical implications. In theory, implicit balancing induces implicit regularizations effects. For instance, imposing ridge penalizations on all parameter matrices in a low-rank approximation is equivalent to using a sparsity-inducing penalization on rank-one components of the model. In practice, enforcing this implicit balancing helps alternating algorithms efficiently minimize the cost function without hurting convergence guarantees. This is demonstrated formally on a toy problem, and practically in three different showcases, namely sparse Nonnegative Matrix Factorization, ridge Nonnegative Canonical Polyadic Decomposition, and sparse Nonnegative Tucker Decomposition.

Our findings open several avenues for future exploration. First, the implicit equivalent formulation of regularized low-rank approximations could be leveraged to derive guarantees regarding the solution space of these models. For instance, it remains unclear under what conditions $\ell_1$-regularized Nonnegative Matrix Factorization yields zero, sparse, or dense solutions. Additionally, our framework is currently limited to positive homogeneous and positive-definite regularizations applied column-wise to parameter matrices of scale-invariant models. While positive homogeneity and positive definiteness allow for a wide range of regularizations, including $\ell_p^p$ norms, some regularizations and constraints commonly used in the literature do not satisfy these criteria. Thus, it would be valuable to investigate how relaxing these assumptions affects our findings.

# Acknowledgements

The authors would like to thank Rémi Gribonval for an early assessment of this work and several pointers toward relevant literature about neural network training. The authors also thank Cédric Févotte and Henrique De Morais Goulart for helpful discussions around the implicit and explicit formulations, as well as pointers toward relevant literature in sNMF.
<figure><img src="image_9.png" /><figcaption>Figure 5: Results for the sNTD experiment on music redundancy detection. **Top:** the loss function for each tested algorithm is shown across iterations for various values of the regularization hyperparameters. **Bottom:** The reconstructed matrix $X_3$ is shown, transposed, for each tested algorithm and various values of the regularization hyperparameter. Green vertical bars indicate an expert segmentation of the song (Come-Together). The number of patterns (rows) is hard to estimate a priori, therefore ideally the estimated matrix $X_3$ should have sparse rows to prune unnecessary patterns. Moreover, it should be group sparse between the green bars to allow for easier segmentation by downstream methods such as dynamic programming or K-means.</figcaption></figure>
A HRSI with unregularized modes is ill-posed

First, let us consider without loss of generality the case $n = 2, r = 1$ and $\mu_2 = 0$ for Problem (5):

$$ \phi(X_1, X_2) = f(X_1, X_2) + \mu_1 g(X_1) \quad (34) $$

Function $\phi$ is bounded below by zero because $f$ is nonnegative, therefore it must admit an infimum over $\mathbb{R}^{m_1 \times r} \times \mathbb{R}^{m_2 \times r}$. Let us denote by $p$ this infimum.

We are going to show that, given $l := \inf_{X_1, X_2} f(X_1, X_2)$, it holds that $p = l$. We know that $p \ge l$ since both $g$ and $\mu_1$ are nonnegative. To prove the converse inequality, note that from Equation (34), for any $\lambda < 1$, it holds that $\phi(X_1, X_2) > \phi(\lambda X_1, \frac{1}{\lambda} X_2)$. Therefore for any admissible $(X_1, X_2)$,

$$ p \le \lim_{\lambda \to 0} \phi(\lambda X_1, \frac{1}{\lambda} X_2). \quad (35) $$

Now we may notice that $\phi(\lambda X_1, \frac{1}{\lambda} X_2) = f(X_1, X_2) + o(\lambda)$, and therefore prolonging by continuity the right hand side, we get that $\lim_{\lambda \to 0} \phi(\lambda X_1, \frac{1}{\lambda} X_2) = f(X_1, X_2)$. This shows that $p \le l$, which concludes our proof that $p = l$.

One can easily note that for the infimum $p$ to be attained, say at $(X_1^*, X_2^*)$, one must have $\phi(X_1^*, X_2^*) = f(X_1^*, X_2^*)$ and therefore $g(X_1^*) = 0$, i.e. $X_1^* = 0$ and $p = f(0, X_2^*)$ since $g_1$ is positive definite. By contraposition, as soon as $f(0, X_2^*)$ is not the infimum of $f$ on the set of matrices with rank at most $r$, the infimum $p$ cannot be attained. Finally, note that $f(0, X_2^*) = f(0, 0)$ by scale-invariance.

B Closed-form expression to solutions of the balancing procedure

We are interested in computing the solution to the following problem:

$$ \min_{\forall i \le n, \lambda_i} \sum_{i=1}^{n} \lambda_i^{p_i} a_i \quad \text{such that} \quad \prod_{i=1}^{n} \lambda_i = 1 \quad (36) $$

Note that we dropped the nonnegativity constraints which are not required to reach our closed-form solution. A particular solution will be naturally nonnegative as long as all products $p_i a_i$ are nonnegative. We also suppose that each constant $a_i$ is strictly positive (see Section 2.1 for a justification).

We use the KKT conditions to compute the solutions, as is traditionally done for the related variational formulation of the geometric mean. Indeed,

$$ \argmin_{\{\lambda_i\}_{i \le n}} \sum_{i=1}^{n} \lambda_i \quad \text{such that} \quad \prod_{i=1}^{n} \lambda_i = \prod_{i=1}^{n} x_i \quad (37) $$

is an alternative definition of the geometric mean $\beta = (\prod_{i \le n} x_i)^{\frac{1}{n}}$.

Let us define the Lagrangian for our problem,

$$ \mathcal{L}(\{\lambda_i\}_{i \le n}, \nu) = \sum_{i \le n} \lambda_i^{p_i} a_i + \nu \left(1 - \prod_{i \le n} \lambda_i\right). \quad (38) $$

Optimality conditions can be written as follows:

$$ \left\{ \begin{array}{l} \nabla_{\lambda_i} \mathcal{L}(\{\lambda_i^*\}_{i \le n}, \nu^*) = p_i a_i \lambda_i^{p_i - 1} - \nu^* \prod_{j \ne i} \lambda_j^* = 0 \\ 1 = \prod_{i \le n} \lambda_i^* \end{array} \right. \quad (39) $$

Using the second equality, we can inject $\lambda_i^*$ in the right hand side of the first equation, which yields the following necessary condition

$$ \lambda_i^* = \left( \frac{\nu^*}{p_i a_i} \right)^{\frac{1}{p_i}}. \quad (40) $$

There may be a sign ambiguity depending on the power value $p_i$, in which case we can chose the positive solution since the original problem involved nonnegativity constraints on $\lambda_i$. Moreover, we have used the facts that $p_i a_i$ are nonzero, and that $\lambda_i^*$ are also nonzero. This last condition can easily be checked since the constraint $\prod_{i \le n} \lambda_i^* = 1$ prevents any optimal $\lambda_i^*$ to be null.
To compute the optimal dual parameter, we can use again the primal feasibility condition:

$$ \prod_{i \le n} \left( \frac{\nu^*}{p_i a_i} \right)^{\frac{1}{p_i}} = 1 \qquad (41) $$

which, after some simple computation, shows that $\nu^*$ is exactly the weighted geometric mean of products $p_i a_i$ with weights $\frac{1}{p_i}$,

$$ \nu^* = \left( \prod_{i \le n} (p_i a_i)^{\frac{1}{p_i}} \right)^{\frac{1}{\sum_{i \le n} \frac{1}{p_i}}} \qquad (42) $$

Note that when applying this balancing operation to an optimal solution $X_i^*$ of HRSI, the scales $\lambda_i$ must equal one, therefore at optimality,

$$ p_i a_i = \nu^* \qquad (43) $$

If the solution $X_i^*$ is null, then $\nu^*$ is also null and the result holds.

## C Proof of sublinear convergence rate for ALS on the scaling problem

We will study the variation of loss function when $x_1$ is updated, $x_2$ fixed, across one iteration. The problem is symmetric therefore we will generalize for all updates of either $x_1$ or $x_2$.

Let us introduce a few notations. Let $k$ the iteration index, we set $x_i := x_i^{(k)}, \tilde{x}_1 = x_1^{(k+1)}$. We denote the residuals $r_k = y - x_1 x_2$ and $r_{k+1} = y - \tilde{x}_1 x_2$, and the regularizations $g_k = \lambda(x_1^2 + x_2^2)$ and $g_{k+1} = \lambda(\tilde{x}_1^2 + x_2^2)$. The cost function value at iteration $k$ is $f_k = r_k^2 + g_k$. Finally, we denote $\delta r_k = r_k - r_{k+1}$ the variation in residuals. Throughout the following derivations, we suppose that $\lambda$ is much smaller than $y$, but also $x_1^2$ and $x_2^2$. Since both $x_1^2$ and $x_2^2$ will eventually converge towards $y - \lambda$, this is equivalent to assuming that $\lambda$ is much smaller than $y$, and $k$ is large enough so that $x_1$ and $x_2$ are relatively close to $\sqrt{y - \lambda}$.

We study the variation of the cost $\delta f_k = f_k - f_{k+1}$. We observe that

$$ f_{k+1} = (-\delta r_k + r_k)^2 + g_{k+1} - g_k + g_k \qquad (44) $$

$$ = f_k - 2r_k \delta r_k + (\delta r_k)^2 + (g_{k+1} - g_k), \qquad (45) $$

and therefore $\delta f_k = 2r_k \delta r_k - (\delta r_k)^2 + (g_k - g_{k+1})$. We can expend the update of the regularization terms:

$$ g_k - g_{k+1} = \lambda (x_1^2 - \tilde{x}_1^2) \qquad (46) $$

$$ = \frac{\lambda}{x_2^2} ((y^2 - \tilde{x}_1^2 x_2^2) - (y^2 - x_1^2 x_2^2)) \qquad (47) $$

$$ = \frac{\lambda}{x_2^2} (r_{k+1}(y + \tilde{x}_1 x_2) - r_k(y + x_1 x_2)) \qquad (48) $$

Recall that $\tilde{x}_1 = \frac{y}{x_2} \frac{1}{1+\frac{\lambda}{x_2^2}}$. We perform a first approximation here of $\tilde{x}_1(\lambda)$ by using its Taylor expansion around $\lambda = 0$, which hold approximately when $\frac{\lambda}{x_2^2} \ll 1$. Then we get

$$ \tilde{x}_1 = \frac{y}{x_2} \left(1 - \frac{\lambda}{x_2^2} + \mathcal{O}(\lambda^2)\right) = \frac{y}{x_2} - \frac{y\lambda}{x_2^3} + \mathcal{O}(\lambda^2), \qquad (49) $$

This means that

$$ y + \tilde{x}_1 x_2 = 2y - \frac{y\lambda}{x_2^2} + \mathcal{O}(\lambda^2). \qquad (50) $$

Furthermore, under the same hypothesis, we can easily see that the residuals is of order of magnitude $\lambda$:

$$ r_{k+1} = y - \tilde{x}_1 x_2 \qquad (51) $$

$$ = \frac{y\lambda}{x_2^2} + \mathcal{O}(\lambda^2) \qquad (52) $$

<footnote n="11">¹¹ALS algorithm updates each variable as follows $\tilde{x}_1 = \argmin_{x_1} (y - x_1 x_2)^2 + \lambda(x_1^2 + x_2^2)$, solving the KKT condition leads to $\tilde{x}_1 = \frac{y x_2}{\lambda + x_2} = \frac{y}{x_2} \frac{1}{1 + \frac{\lambda}{x_2^2}}$.</footnote>
and therefore, we can relate the sums $y + x_1x_2$ and $y + \tilde{x}_1x_2$ to the residuals as

$$y + \tilde{x}_1 x_2 = 2y - r_{k+1} + \mathcal{O}(\lambda^2) \qquad (53)$$

$$y + x_1 x_2 = 2y - r_k + \mathcal{O}(\lambda^2). \qquad (54)$$

This yields the following formula for the variation of regularization terms

$$g_k - g_{k+1} = \frac{\lambda}{x_2^2} (2y(r_{k+1} - r_k) + (r_k^2 - r_{k+1}^2) + \mathcal{O}(\lambda^3)) \qquad (55)$$

$$= -2r_{k+1}\delta r_k + \mathcal{O}(\lambda^3). \qquad (56)$$

We have assumed that $\delta r_k$ is also of the order of magnitude $\lambda$, which is shown in the derivations below. Finally, the variation in loss function can be written as

$$\delta f_k = -2r_{k+1}\delta r_k + 2r_k\delta r_k - (\delta r_k)^2 + \mathcal{O}(\lambda^3) \qquad (57)$$

$$= 2\delta r_k(r_k - r_{k+1}) - (\delta r_k)^2 + \mathcal{O}(\lambda^3) \qquad (58)$$

$$= (\delta r_k)^2 + \mathcal{O}(\lambda^3) \qquad (59)$$

To further detail the expression of $\delta r_k$, we may first relate $\tilde{x}_1$ with $x_1$ since

$$\tilde{x}_1 = \frac{y}{x_2} \frac{x_2^2}{x_2^2 + \lambda} = \frac{y}{x_2} \left( 1 - \frac{\lambda}{x_2^2 + \lambda} \right). \qquad (60)$$

Viewing $x_2$ as the least squares update from $x_1$, we further obtain

$$\tilde{x}_1 = x_1 \frac{1 - \frac{\lambda}{x_2^2 + \lambda}}{1 - \frac{\lambda}{x_1^2 + \lambda}} \qquad (61)$$

Therefore we get

$$\delta r_k^2 = x_2^2(x_1 - \tilde{x}_1)^2 \qquad (62)$$

$$= x_1^2 x_2^2 \left( 1 - \frac{1 - \frac{\lambda}{x_1^2 - \lambda}}{1 - \frac{\lambda}{x_2^2 + \lambda}} \right)^2 \qquad (63)$$

$$= x_1^2 x_2^2 \lambda^2 \left( \frac{\frac{1}{x_2^2 + \lambda} - \frac{1}{x_1^2 + \lambda}}{1 - \frac{\lambda}{x_2^2 + \lambda}} \right)^2 \qquad (64)$$

The denominator can be neglected since it will generate terms in the order of $\lambda^4$ as soon as $\lambda$ is small with respect to $x_2^2$, which is our working hypothesis. This yields

$$\delta r_k^2 = x_1^2 x_2^2 \lambda^2 \left( \frac{x_1^2 - x_2^2}{(x_1^2 + \lambda)(x_2^2 + \lambda)} \right)^2 + \mathcal{O}(\lambda^3) \qquad (65)$$

$$= \frac{x_1^2 x_2^2 (x_1 + x_2)^2}{(x_1^2 + \lambda)^2 (x_2^2 + \lambda)^2} \lambda^2 (x_1 - x_2)^2 + \mathcal{O}(\lambda^3) \qquad (66)$$

Let us introduce the errors $e_1, e_2$ from $x_1, x_2$ to their target $y - \lambda$:

$$x_1 = \sqrt{y - \lambda} + e_1 \qquad (67)$$

$$x_2 = \sqrt{y - \lambda} + e_2 \qquad (68)$$

We may note that $e_1$ and $e_2$ are of opposite sign as soon as $k > 1$, and for simplicity of presentation, we will assume that $e_1 \approx -e_2 = e_k$. We further assume that $k$ is large enough so that $e_1, e_2$ are small enough (of order of magnitude $\mathcal{O}(\lambda)$) to approximate $(x_1 + x_2)^2$ as $4y$ and $(x_1 + \lambda)^2, (x_2 + \lambda)^2$ as $y$ up to an error $\mathcal{O}(\lambda)$. The term $x_1 - x_2$ is of order $\mathcal{O}(\lambda)$. Therefore for large enough $k$,

$$\delta f_k = y^2 \frac{4y}{y^2 y^2} \lambda^2 (x_1 - x_2)^2 + \mathcal{O}(\lambda^3) \qquad (69)$$

$$= 4 \frac{\lambda^2}{y} (x_1 - x_2)^2 + \mathcal{O}(\lambda^3) \qquad (70)$$

$$= 16 \frac{\lambda^2}{y} e_k^2 + \mathcal{O}(\lambda^3) \qquad (71)$$
We are left to characterize the convergence speed of $e_k$ towards zero. We can see that $e_k$ converges linearly with however a multiplicative constant that goes to one as $\frac{\lambda}{y}$ grows smaller. Indeed,

$$e_{k+1} = \tilde{x}_1 - x^* \qquad (72)$$

$$= x_1 - x^* - x_1 \left( 1 - \frac{1 - \frac{\lambda}{x_2^2 + \lambda}}{1 - \frac{\lambda}{x_1^2 + \lambda}} \right) \qquad (73)$$

$$= x_1 - x^* - x_1 \lambda \left( \frac{x_1^2 - x_2^2}{(x_1^2 + \lambda)(x_2^2 + \lambda)} \right) + \mathcal{O}(\lambda^2) \qquad (74)$$

$$= x_1 - x^* - \frac{2\lambda}{y}(x_1 - x_2) + \mathcal{O}(\lambda^2) \qquad (75)$$

where the last two lines are obtained by reusing the previous derivations for $\delta r_k$. We find that

$$e_{k+1} = e_k - e_k 4 \frac{\lambda}{y} + \mathcal{O}(\lambda^2) \qquad (76)$$

We see that $\frac{e_{k+1}}{e_k}$ is approximately constant, equal to $1 - 4\frac{\lambda}{y}$, which proves that ALS is a sublinear algorithm when $\frac{\lambda}{y}$ goes to zero.

## D Convergence of iterates for Algorithm 2

In this supplementary section we discuss in details the convergence guarantees offered by Algorithm 2 w.r.t the sequence of objective functions and the sequence of iterates.

For convenience, let us briefly recall the formulation of our optimization problems as a general class of multiblock optimization problems; denoting $z = (\{X_i[\cdot, q]\}_{i \le n, q \le r})$ (concatenation of the variables) and $\phi$ the objective function from Problem (1), we consider the problems

$$\min_{z_i \in \mathcal{Z}_i} \phi(z_1, \dots, z_s), \qquad (77)$$

where $z$ can be decomposed into $s$ blocks $z = (z_1, \dots, z_s)$, $\mathcal{Z}_i \subseteq \mathbb{R}^{m_i}$ is a closed convex set, $\mathcal{Z} = \mathcal{Z}_1 \times \dots \mathcal{Z}_s \subseteq \text{dom}(\phi)$, $\phi: \mathbb{R}^m \to \mathbb{R} \cup \{+\infty\}$, $m = \sum_{i=1}^s m_i$, is a continuous differentiable function over its domain, and $\phi$ is lower bounded and has bounded level sets.

In Theorem 3 we show that Algorithm 2 generates a sequence of iterates $\{z^i\}_{i=1}^\infty$ such that the sequence of objective values $\{\phi(z^i)\}$ is monotone non-increasing. First, let us recall the formal definition of a tight majorizer for convenience.

**Definition 1.** A function $\bar{\phi}_i: \mathcal{Z}_i \times \mathcal{Z} \to \mathbb{R}$ is said to be a tight majorizer of $\phi(z)$ if

$$1. \quad \bar{\phi}_i(z_i, z) = \phi(z), \forall z \in \mathcal{Z},$$

$$2. \quad \bar{\phi}_i(y_i, z) \ge \phi(z_1, \dots, z_{i-1}, y_i, z_{i+1}, \dots, z_s), \forall y_i \in \mathcal{Z}_i, \forall z \in \mathcal{Z},$$

where $z_i$ denotes the $i$-th block of variables of $z$ with $1 \le i \le s$.

At iteration $k+1$, Algorithm 2 fixes the last values of block of variables $j \ne i$ and updates block $z_i$ as follows:

$$z_i^{(k+1)} = \underset{z_i \in \mathcal{Z}_i}{\operatorname{argmin}} \bar{\phi}_i(z_i, z_1^{(k+1)}, \dots, z_{i-1}^{(k+1)}, z_i^{(k,*)}, \dots, z_s^{(k,*)})$$

where $z_s^{(k,*)}$ denotes the $s$-th block of variable computed at iteration $k$ and optimally scaled using Equations (9). After updating all the blocks as described above, Algorithm 2 performs optimal scaling on each block. We first present Proposition 4 which will be useful for the proof of Theorem 3.

**Proposition 4.** Let $\{X_i\}_{i \le n}$ be the current decomposition, the balanced columns $X_i[\cdot, q]^*$ computed according to Algorithm 1 necessarily decrease the objective function value of Problem (77).

*Proof.* The proof directly follows the results presented in Section 2.2.2, indeed we have:

$$\begin{align}
\phi(\{X_i\}_{i \le n}) &\stackrel{\Delta}{=} \sum_{\forall i \le n, \Lambda_i = I_r \times r} f(\{X_i \Lambda_i\}_{i \le n}) + \sum_{i=1}^n \mu_i \sum_{q=1}^r g_i(\Lambda_i[q, q] X_i[\cdot, q]) \tag{78} \\
&\ge f(\{X_i \Lambda_i^*\}_{i \le n}) + \sum_{i=1}^n \mu_i \sum_{q=1}^r g_i(\Lambda_i^*[q, q] X_i[\cdot, q]) = \phi(\{X_i^*\}_{i \le n}). \nonumber
\end{align}$$
where $\Lambda_i^*[q, q] = \left( \frac{\beta_q}{p_i \mu_i g_i(X_i[\cdot, q])} \right)^{1/p_i}$ from Equations (9) with $\beta_q$ as defined in Proposition 2. The second inequality holds by definition of $\Lambda_i^*[q, q]$, that is the optimal and unique solution of the scaling Problem (10). $\square$

**Theorem 3** (also Theorem 1). Let $z^{(k)} = (\{X_i^{(k)}[\cdot, q]\}_{i \le n, q \le r}) \ge 0$, and let $\bar{\phi}_i$ be tight majorizers, as per Definition 1, for $\phi: \mathcal{Z} \to \mathbb{R}_+$ from Problem (77), with $\mathcal{Z} := \mathcal{Z}_1 \times \dots \mathcal{Z}_s \subseteq \text{dom}(\phi)$ and $1 \le i \le s = nr$. Then $\phi$ is monotone non-increasing under the updates of Algorithm 2. Moreover, the sequences $\{\phi(z^{(k)},*)\}$ and $\{\phi(z^{(k)})\}$ converge.

*Proof.* From Definition 1, we have:

$$ \begin{aligned} \phi(z^{(k),*}) &= \bar{\phi}_1(z_1^{(k),*}, z_1^{(k),*}, \dots, z_s^{(k),*}) \ge \bar{\phi}_1(z_1^{(k+1)}, z_1^{(k),*}, \dots, z_s^{(k),*}) \\ &\ge \phi(z_1^{(k+1)}, z_2^{(k),*}, \dots, z_s^{(k),*}) = \bar{\phi}_2(z_2^{(k),*}, z_1^{(k+1)}, z_2^{(k),*}, \dots, z_s^{(k),*}) \\ &\ge \bar{\phi}_2(z_2^{(k+1)}, z_1^{(k+1)}, z_2^{(k),*}, \dots, z_s^{(k),*}) \ge \phi(z_1^{(k+1)}, z_2^{(k+1)}, \dots, z_s^{(k),*}) \ge \dots \\ &\ge \phi(z_1^{(k+1)}, z_2^{(k+1)}, \dots, z_s^{(k+1)}) \ge \phi(z_1^{(k+1),*}, z_2^{(k+1),*}, \dots, z_s^{(k+1),*}) \end{aligned} $$

where $z_i^{(k+1),*}$ are computed using Equations (9) for $1 \le i \le s = nr$.

By assumption, $\phi$ from Problem (77) is a finite nonnegative sum of functions that are all bounded below onto the feasible set $\mathcal{Z}$, hence $\phi$ is also bounded below onto $\mathcal{Z}$. Therefore, both the non-increasing monotone sequences $\{\phi(z^{(k),*})\}$ and $\{\phi(z^{(k)})\}$ converge. $\square$

**Convergence of the sequence of iterates** In the next Theorem, we establish the convergence of the iterates generated by Algorithm 2 towards coordinate-wise minimum, or critical points of Problem (77) under mild assumptions.

**Theorem 4** (also Theorem 2). Let majorizers $\bar{\phi}_i(y_i, z)$, as defined in Definition 1, be quasi-convex in $y_i$ for $i = 1, \dots, s$. Additionally, assume that Assumption 2 are satisfied. Consider the subproblems in Equation (81), and assume that both these subproblems and the optimal scaling, obtained using Equations (9), have unique solutions for any point

$$ (z_1^{(k+1)}, \dots, z_{i-1}^{(k+1)}, z_i^{(k),*}, \dots, z_s^{(k),*}) \in \mathcal{Z}. $$

Under these conditions, every limit point $z^\infty$ of the iterates $\{z^{(k)}\}_{k \in \mathbb{N}}$ generated by Algorithm 2 is a coordinate-wise minimum of Problem (77). Furthermore, if $\phi(\cdot)$ is regular at $z^\infty$, then $z^\infty$ is a stationary point of Problem (77).

The proof is handled in Appendix E below, and is an adaptation of the convergence to stationary points obtained in [Razaviyayn et al., 2013, Theorem 2] for block successive upperbound minimization (BSUM) framework.

## E Convergence proof for the iterates of the Meta-Algorithm

To better analyze the convergence of the iterates generated by Algorithm 2, it is important to review some key concepts presented in Razaviyayn et al. [2013]:

* **Directional derivative:** For a function $\phi: \mathcal{Z} \to \mathbb{R}$, where $\mathcal{Z} \subseteq \mathbb{R}^m$ is a convex set, the directional derivative of $\phi$ at point $z$ along a vector $d$ is defined as:

$$ \phi'(z; d) := \lim_{\lambda \to 0} \frac{\phi(z + \lambda d) - \phi(z)}{\lambda} \qquad (79) $$

* **Stationary points of a function:** For a function $\phi: \mathcal{Z} \to \mathbb{R}$, where $\mathcal{Z} \subseteq \mathbb{R}^m$ is a convex set, the point $z^*$ is considered a stationary point of $\phi(\cdot)$ if $\phi'(z^*; d) \ge 0$ for all $d$ such that $z^* + d \in \mathcal{Z}$.

* **Quasi-convex function:** A function $\phi$ is said to be quasi-convex if it satisfies the following inequality:

$$ \phi(\theta x + (1 - \theta)y) \le \max(\phi(x), \phi(y)) \quad \forall \theta \in (0, 1), \quad \forall x, y \in \text{dom}\phi \qquad (80) $$

* **Regularity of a function at a point:** The function $\phi: \mathcal{Z} \to \mathbb{R}$ is regular at the point $z \in \text{dom}\phi$ w.r.t. coordinates $m_1, \dots, m_s$ if $\phi'(z; d) \ge 0$ for all $d = (d_1, \dots, d_s)$ with $\phi'(z; d_i^0) \ge 0$, where $d_i^0 := (0, \dots, d_i, \dots, 0)$ and $d_i \in \mathbb{R}^{m_i}$ for all $i$.
* Coordinate-wise minimum of a function $z \in \text{dom}\phi \in \mathbb{R}^m$ is the coordinate-wise minimum of $\phi$ w.r.t. the coordinates in $\mathbb{R}^{m_1}, ..., \mathbb{R}^{m_s}$ if

$$ \phi(z + d_i^0) \geq \phi(z) \quad \forall d_i \in \mathbb{R}^{m_i} \quad \text{with } z + d_i^0 \in \text{dom}\phi \quad \forall i $$

or equivalently:

$$ \phi'(z; d_i^0) \geq 0 \quad \forall d_i \in \mathbb{R}^{m_i} \quad \text{with } z + d_i^0 \in \text{dom}\phi \quad \forall i $$

* Block Successive MM: Algorithm 2 updates only one block of variables at a time. Specifically, in iteration $k+1$, the selected block (referred to as block $i$) is updated by solving the subproblem:

$$ z_i^{(k+1)} = \operatorname*{argmin}_{z_i \in \mathcal{Z}_i} \bar{\phi}_i(z_i, z_1^{(k+1)}, \dots, z_{i-1}^{(k+1)}, z_i^{(k),*}, \dots, z_s^{(k),*}) \quad (81) $$

where $\bar{\phi}_i(\cdot, z)$ is a tight majorizer of $\phi(z)$ as defined in Definition 1.

To discuss the convergence of Algorithm 2, additional assumptions about the majorizers $\bar{\phi}_i$ are necessary. These assumptions are as follows:

**Assumption 2.**

1. $\bar{\phi}'_i(y_i, z; d_i)|_{y_i=z_i} = \phi'(z; d_i^0) \quad \forall d_i^0 = (0, \dots, d_i, \dots, 0)$ such that $z_i + d_i \in \mathcal{Z}_i, \forall i$. This assumption implies that the directional derivatives of the majorizers and the objective function coincide at the current iterate for each block of variables.

2. $\bar{\phi}_i(y_i, z)$ is continuous in $(y_i, z) \quad \forall i$.

Following the approach outlined in Razaviyayn et al. [2013], the convergence outcomes of Algorithm 2 are applicable in the asymptotic regime as $k \to \infty$. This convergence relies on the aforementioned assumptions, which include the directional differentiability and the continuity of the function $\phi$. In addition to this, the proof leverages two critical assumptions: firstly, the quasi-convexity of the majorizers $\bar{\phi}_i(y_i, z)$, and secondly, the uniqueness of their minimizers. Here-under, we give the proof of Theorem 4.

*Proof.* We follow the main of steps of the general proof of BSUM Razaviyayn et al. [2013], and slightly adapt them to the specificity of our Algorithm. First, by Theorem 3, we have:

$$ \phi(z^{(0)}) \geq \phi(z^{(1)}) \geq \phi(z^{(1),*}) \geq \phi(z^{(2)}) \geq \phi(z^{(2),*}) \geq \dots \quad (82) $$

Consider a limit point $z^\infty$, by combining Equation (82) and the continuity of $\phi$, we have:

$$ \lim_{k \to \infty} \phi(z^{(k)}) = \phi(z^\infty) = \lim_{k \to \infty} \phi(z^{(k),*}) \quad (83) $$

In the ensuing discussion, we focus on the subsequence denoted as $\{z^{(k_j),*}\}$, which converges to the limit point $z^\infty$. Moreover, given the finite number of variable blocks, there must exist at least one block that undergoes infinite updates within the subsequences $\{k_j\}$. Without loss of generality, we assume block $s$ is updated infinitely often. Consequently, by constraining our attention to these subsequences, we can express this as:

$$ z_s^{(k_j+1)} = \operatorname*{argmin}_{z_s \in \mathcal{Z}_s} \bar{\phi}_s(z_s, z_1^{(k_j+1)}, z_2^{(k_j+1)}, \dots, z_s^{(k_j),*}) \quad (84) $$

To establish the convergence $z^{(k_j+1)} \to z^\infty$ as $j \to \infty$, or equivalently, to demonstrate $z_i^{(k_j+1)} \to z_i^\infty$ for all $i$ as $j \to \infty$, we consider the block $i=1$. Assume the contrary, supposing $z_1^{(k_j+1)}$ does not converge to $z_1^\infty$. Consequently, by narrowing our focus to a subsequence, there exists $\tilde{\gamma} > 0$ such that:

$$ \tilde{\gamma} \leq \gamma^{(k_j)} = \|z_1^{(k_j+1)} - z_1^{(k_j),*}\| \quad \forall k_j $$

This implies that as $j \to \infty$, $z_1^{(k_j+1)}$ consistently deviates from $z_1^{(k_j),*}$, which converges to $z_1^\infty$. Consequently, by transitivity, $z_1^{(k_j+1)}$ does not converge to $z_1^\infty$.

We introduce the sequence $s^{(k_j)}$, defined as the normalized difference between $z_1^{(k_j+1)}$ and $z_1^{(k_j),*}$:
$$s(k_j) := \frac{z_1^{(k_j+1)} - z_1^{(k_j),*}}{\gamma^{(k_j)}} \qquad (85)$$

Notably, we observe that $||s^{(k_j)}|| = 1$ holds for all $\{k_j\}$; thus, $s^{(k_j)}$ is a member of a compact set and consequently possesses a limit point denoted as $\bar{s}$. By further narrowing our focus to a subsequence that converges to $\bar{s}$, we can express the ensuing set of inequalities:

$$\phi(z^{(k_j+1),*}) \le \phi(z^{(k_j+1)}) \le \phi(z^{(k_j+1),*}) \le \phi(z^{(k_j+1)}) \qquad (86)$$

which hold thanks to Equation (82) (consequence of Theorem 3) and given that $k_{j+1} \ge k_j + 1$. Moreover:

$$
\begin{align*}
\phi(z^{(k_j+1)}) &\le \bar{\phi}_1(z_1^{(k_j+1)}, z^{(k_j),*}) \quad \text{(Definition 1)} \\
&= \bar{\phi}_1(z_1^{(k_j),*} + \gamma^{(k_j)} s^{(k_j)}, z^{(k_j),*}) \quad \text{(Using (85))} \\
&\le \bar{\phi}_1(z_1^{(k_j),*} + \epsilon \bar{\gamma} s^{(k_j)}, z^{(k_j),*}) \quad (\epsilon \in [0, 1] \text{ and } z_1^{(k_j+1)} = \operatorname*{argmin}_{\gamma^{(k_j)}} \bar{\phi}_1(\cdot, z^{(k_j),*})) \\
&:= \bar{\phi}_1(z_1^{(k_j),*} + \theta(z_1^{(k_j+1)} - z_1^{(k_j),*}), z^{(k_j),*}) \quad \text{(With } \theta = \epsilon \frac{\bar{\gamma}}{\gamma^{(k_j)}} \in [0, 1]\text{)} \tag{87} \\
&\le \max \left( \bar{\phi}_1(z_1^{(k_j+1)}, z^{(k_j),*}), \bar{\phi}_1(z_1^{(k_j),*}, z^{(k_j),*}) \right) \quad \text{(quasi-convexity of } \bar{\phi}_1(\cdot, z^{(k_j),*})\text{)} \\
&= \bar{\phi}_1(z^{(k_j),*}, z^{(k_j),*}) \quad \text{(by definition of the arg min } z_1^{(k_j+1)}\text{)} \\
&= \phi(z^{(k_j),*}) \quad \text{(Majorizer is tight at } z^{(k_j),*}\text{)}
\end{align*}
$$

Letting $j \to \infty$ and combining Equations (86) and (87), we have:

$$\phi(z^{\infty}) \le \bar{\phi}_1(z_1^{\infty} + \epsilon \bar{\gamma} \bar{s}) \le \phi(z^{\infty}) \qquad (88)$$

or equivalently:

$$\phi(z^{\infty}) = \bar{\phi}_1(z_1^{\infty} + \epsilon \bar{\gamma} \bar{s}) \qquad (89)$$

Furthermore:

$$
\begin{align}
\bar{\phi}_1(z_1^{(k_j+1),*}, z^{(k_j+1),*}) &= \phi(z^{(k_j+1),*}) \le \phi(z^{(k_j+1),*}) \le \phi(z^{(k_j+1)}) \nonumber \\
&\le \bar{\phi}_1(z_1^{(k_j+1)}, z^{(k_j),*}) \tag{90} \\
&\le \bar{\phi}_1(z_1, z^{(k_j),*}) \quad \forall z_1 \in \mathbb{Z}_1 \nonumber
\end{align}
$$

Thanks to the continuity of majorizers, and letting $j \to \infty$, we obtain:

$$\bar{\phi}_1(z_1^{\infty}, z^{\infty}) \le \bar{\phi}_1(z_1, z^{\infty}) \quad \forall z_1 \in \mathbb{Z}_1$$

which implies that $z_1^{\infty}$ is the minimizer of $\bar{\phi}_1(\cdot, z^{\infty})$. By hypothesis, we assumed that the minimizer of each majorizer is unique, which is a contradiction with Equation (89). Therefore, the contrary assumption is not true, i.e., $z^{(k_j+1)} \to z^{\infty}$ for $j \to \infty$.

Since $z_1^{(k_j+1)} = \operatorname*{argmin}_{z_1 \in \mathbb{Z}_1} \bar{\phi}_1(z_1, z^{(k_j),*})$, then we have:

$$\bar{\phi}_1(z_1^{(k_j+1)}, z^{(k_j),*}) \le \bar{\phi}_1(z_1, z^{(k_j),*}) \quad \forall z_1 \in \mathbb{Z}_1$$

The computation of the limit $j \to \infty$ gives:

$$\bar{\phi}_1(z_1^{\infty}, z^{\infty}) \le \bar{\phi}_1(z_1, z^{\infty}) \quad \forall z_1 \in \mathbb{Z}_1$$

The rest of the proof is identical to Razaviyayn et al. [2013], we detail the remaining steps for the sake of completeness. The above inequality implies that all the directional derivatives of $\phi_1(\cdot, z^{\infty})$ w.r.t. block 1 of variables, and evaluated at $z_1^{\infty}$ are nonnegative, that is:

$$\bar{\phi}'(z_1, z^{\infty}; d_1) \Big|_{z_1=z_1^{\infty}} \ge 0 \quad \forall d_1 \in \mathbb{R}^{m_1} \text{ such that } z_1^{\infty} + d_1 \in \mathbb{Z}_1$$
Following the same rationale for the other blocks of variables, we obtain:

$$ \bar{\phi}'(z_i, z^\infty; d_i) \Big|_{z_i=z_i^\infty} \ge 0 \quad \forall d_i \in \mathbb{R}^{m_i} \text{ such that } z_i^\infty + d_i \in \mathcal{Z}_i, \forall i = 1, \dots, s \quad (91) $$

Using Assumption 2 and Equations (91), we obtain:

$$ \phi'(z^\infty; d_i^0) \ge 0 \quad \forall d_i^0 = (0, \dots, d_i, \dots, 0) \text{ such that } z_i^\infty + d_i \in \mathcal{Z}_i, \forall i = 1, \dots, s \quad (92) $$

Equations (92) implies that $z^\infty$ is a coordinate-wise minimum of $\phi$. Finally, assuming that $\phi$ is regular on its domain, including $z^\infty$, then we can write:

$$ \phi'(z^\infty; d) \ge 0 \quad \forall d = (d_1, \dots, d_s) \text{ such that } z^\infty + d \in \mathcal{Z} \quad (93) $$

Under this additional regularity assumption of $\phi$, Equation (93) implies that the limit point $z^\infty$ is a stationary point of $\phi$. This concludes the proof. □

## F MM update rules for the meta-algorithm

We consider a generic minimization problem of the form

$$ \underset{X_1 \ge 0}{\operatorname{argmin}} D_{\beta}(T|X_1 U) + \mu \sum_{q \le r} \|X_1[:q]\|_p^p \quad (94) $$

where $T = T_{[1]} \in \mathbb{R}_+^{m_1 \times m_2 \times m_3}$, $U = G_{[1]}(X_2 \boxtimes X_3)^T \in \mathbb{R}_+^{r_1 \times m_2 m_3}$, $D_\beta$ is the beta-divergence and $\mu > 0$. Notably, we considered the decomposition of a third-order tensor $T$, but the approach can be generalized to any order.

In Problem (94), the regularization term is separable w.r.t. each entry of $X_1[:q]$ for all $q \le r$. For the data fitting term, we use the majorizer proposed in Févotte and Idier [2011]. For the sake of completeness and to ease the understanding of the further developments, we briefly recall it in the following. It consists in first decomposing the $\beta$-divergence between two nonnegative scalars $a$ and $b$ as a sum of a convex, concave and constant terms as follows:

$$ d_{\beta}(a|b) = \check{d}_{\beta}(a|b) + \hat{d}_{\beta}(a|b) + \bar{d}_{\beta}(a|b), \quad (95) $$

where $\check{d}$ is a convex function of $b$, $\hat{d}$ is a concave function of $b$ and $\bar{d}$ is a constant w.r.t. $b$; see Table 1, and then in majorizing the convex part of the $\beta$-divergence using Jensen's inequality and majorizing the concave part by its tangent (first-order Taylor approximation).

<table><caption>Table 1: Differentiable convex-concave-constant decomposition of the $\beta$-divergence under the form (95) Févotte and Idier [2011].</caption><thead><tr><th></th><th>\tilde{d}(a|b)</th><th>\hat{d}(a|b)</th><th>\bar{d}(a)</th></tr></thead><tbody><tr><td>β = (-∞, 1) \ {0}</td><td>1/(1-β)ab<sup>β-1</sup></td><td>1/βb<sup>β</sup></td><td>1/(β(β-1))a<sup>β</sup></td></tr><tr><td>β = 0</td><td>ab<sup>-1</sup></td><td>log(b)</td><td>a(log(a) - 1)</td></tr><tr><td>β ∈ [1, 2]</td><td>d<sub>β</sub>(a|b)</td><td>0</td><td>0</td></tr></tbody></table>

Since $\mu \sum_{q \le r} \|X_1[:q]\|_p^p = \mu \|X_1\|_p^p$ is convex, the subproblem obtained after majorization (step 1 of MM, see Section 3.2) to solve has the form:

$$ \underset{X_1 \ge 0}{\operatorname{argmin}} G(X_1|\tilde{X}_1) + \mu \|X_1\|_p^p \quad (96) $$

where $G(X_1|\tilde{X}_1)$ is the separable, tight majorizer for $D_\beta(T|X_1 U)$ at $\tilde{X}_1 \ge 0$ of the following form

$$ \begin{align} G(X_1|\tilde{X}_1) &= \sum_{i \le m_1, k \le r} l(X_1[i,k]|\tilde{X}_1) \nonumber \\ \text{such that } l(X_1[i,k]|\tilde{X}_1) &= \sum_{j \le m_2 m_3} \frac{\tilde{X}_1[i,k]U[k,j]}{\tilde{T}[i,j]} \check{d}(T[i,j]|\tilde{T}[i,j]) \frac{X_1[i,k]}{\tilde{X}_1[i,k]} \tag{97} \\ &\quad + \sum_{j \le m_2 m_3} [U[k,j]\hat{d}(T[i,j]|\tilde{T}[i,j])] X_1[i,k] + \text{cst} \nonumber \end{align} $$
according to [Févotte and Idier, 2011, Lemma 1] where $\tilde{T} = \tilde{X}_1 U$.

As explained in Section 3.2, the optimization problem (96) is a particular instance of the family of problems covered by the framework proposed in Leplat et al. [2021] in the case no additional equality constraints are required. Assuming $\mu > 0$ and the regularization terms satisfy Assumption 1, the minimizer of Problem (96) exists and is unique in $(0, \infty)^{m_1 \times r_1}$ as soon as $\beta < 2$.

Moreover, according to Leplat et al. [2021], it is possible to derive closed-form expressions for the minimizer of Problem (96) for the following values of $\beta$: $\beta \in \{0, 1, 3/2\}$ (if constants $C > 0$ in Assumption 1) or $\beta \in (-\infty, 1] \cup \{5/4, 4/3, 3/2\}$ (if constants $C = 0$ in Assumption 1 when $p = 1$ for instance). Outside these values for $\beta$, an iterative scheme is required to numerically compute the minimizer. We can cite the Newton's method, the Muller's method and the procedure developed in Rechtschaffen [2008] which is based on the explicit calculation of the intermediary root of a canonical form of cubic. This procedure is suited for providing highly accurate numerical results in the case of badly conditioned polynomials.

In this section, we will consider the interval $\beta \in [1, 2)$. For this setting, we can show that computing the minimizer of Problem (96) corresponds to solving, for a single element $X_1[i, k]$:

$$ \nabla_{X_1[i,k]} [G(X_1|\tilde{X}_1) + \mu \|X_1\|_p^p] = 0 \qquad (98) \\ \Leftrightarrow aX_1[i, k]^{p-1+2-\beta} + bX_1[i, k] - c = 0 \text{ since } X_1[i, k] \in (0, \infty). $$

where the second equation has been obtained by multiplying both sides by $X_1[i, k]^{\beta-2}$, and:

* $a = p\mu > 0$ by hypothesis,

* $b = \sum_j U[k, j] \left( \frac{\tilde{T}[i, j]}{\tilde{X}_1[i, k]} \right)^{\beta-1} \ge 0$

* $c = \sum_j U[k, j] T[i, j] \left( \frac{\tilde{T}[i, j]}{\tilde{X}_1[i, k]} \right)^{\beta-2} \ge 0$.

At this stage, one may observe that finding the solution $X_1$ boils down to solving several polynomial equations in parallel. Nevertheless, depending on the values of $\beta$ and $p$, we can further simplify the update.

**Case $\beta=1, p=1$** Equation (98) becomes $aX_1[i, k] + bX_1[i, k] - c = 0$. The positive (real) root is computed as follows:

$$ \begin{aligned} \hat{X}[i, k] &= \frac{c|_{\beta=1}}{a + b|_{\beta=1}} \\ &= \tilde{X}_1[i, k] \frac{\sum_j U[k, j] \frac{T[i, j]}{T[i, j]}}{\mu + \sum_j U[k, j]} \end{aligned} \qquad (99) $$

which is nonnegative. In matrix format,

$$ \hat{X}_1 = \tilde{X}_1 \odot \frac{\left[ \frac{T}{[\tilde{X}_1 U]} U^\top \right]}{[\mu e_{r1 \times m_1} + e_{m_1 \times m_2 m_3} U^\top]} \qquad (100) $$

where $\odot$ denotes the Hadamard product, $[A]/[B]$ denotes the elementwise division between matrices $A$ and $B$. Note that in Hoyer [2002], the regularization parameters $\mu$ is mistakenly written in the numerator. This error has been corrected in subsequent works without an explicit mention of the error Le Roux and Hershey [2015].

**Case $\beta=1, p=2$** Equation (98) becomes $aX_1[i, k]^2 + bX_1[i, k] - c = 0$. The positive (real) root is computed as follows:

$$ \hat{X}_1[i, k] = \frac{\sqrt{\left(\sum_j U[k, j]\right)^2 + 8\mu \tilde{X}_1[i, k] \sum_j U[k, j] \frac{T[i, j]}{T[i, j]} - \sum_j U[k, j]}}{4\mu} \qquad (101) $$

with $\mu > 0$. Note that although the closed-form expression in Equation (101) has a negative term in the numerator of the right-hand side, it can be easily checked that it always remains nonnegative given $T[i, j], U[k, j]$ and $\tilde{X}[i, k]$ are nonnegative for all $i, j, k$.

Equation (101) can be expressed in matrix form as follows:
$$ \hat{X}_1 = \frac{[C \cdot 2 + S]^{\frac{1}{2}} - C}{4\mu} \qquad (102) $$

where $[A]^{\alpha}$ denotes the elementwise $\alpha$-exponent, $C = eU^{\top}$ with $e$ is a all-one matrix of size $m_1 \times m_2 m_3$ and $S = 8\mu \tilde{X}_1 \odot \left( \frac{[T]}{[\tilde{X}_1 U]} U^T \right)$.

Some interesting insights are discussed here-under:

* Computing the limit $\lim_{\mu \to 0} \hat{X}_1(\mu)$ makes Equation (102) tends to the original Multiplicative Updates introduced by Lee and Seung Lee and Seung [1999]. Indeed let us compute this limit in the scalar case from Equation (101) and pose $\alpha = \sum_n u[j, k]$ and $\eta = \tilde{w}[j] \sum_j U[k, j] \frac{T[i, j]}{T[i, j]}$ for convenience:

$$ \begin{aligned} \lim_{\mu \to 0} \frac{\sqrt{\alpha^2 + 8\mu\eta - \alpha}}{4\mu} &= \frac{0}{0} \\ &= \lim_{H \to 0} \frac{\frac{\partial}{\partial_{\mu}} \sqrt{\alpha^2 + 8\mu\eta - \alpha}}{\frac{\partial}{\partial_{\mu}} 4\mu} \qquad (103) \\ &= \lim_{\mu \to 0} (\alpha^2 + 8\mu\eta)^{-\frac{1}{2}} \eta = \frac{\eta}{\alpha} = \tilde{X}_1[i, k] \frac{\sum_j U[k, j] \frac{T[i, j]}{T[i, j]}}{\sum_j U[k, j]} \end{aligned} $$

In matrix form we have:

$$ \lim_{\mu \to 0} \hat{X}_1(\mu) = \tilde{X}_1 \odot \frac{\left[ \frac{T}{[\tilde{X}_1 U]} U^T \right]}{[eU^T]} \qquad (104) $$

which are the MU proposed by Lee and Seung for $\beta = 1$.

* Regarding the analysis of univariate polynomial equation given in (98), interestingly the existence of a unique positive real root could have been also established by using Descartes rules of sign. Indeed we can count the number of real positive roots that $p(X_1[i, k]) = aX_1[i, k]^{p-1+2-\beta} + bX_1[i, k] - c$ has (for $\beta \in [1, 2]$). More specifically, let $v$ be the number of variations in the sign of the coefficients $a, b, c$ (ignoring coefficients that are zero). Let $n_p$ be the number of real positive roots. Then:

1. $n_p \le v$,

2. $v - n_p$ is an even integer.

Let us consider the case $\beta = 1$ and $p = 2$ as an example: given the polynomial $p(X_1[i, k]) = aX_1[i, k]^2 + bX_1[i, k] - c$, assuming that $c$ is positive. Then $v = 1$, so $n_p$ is either 0 or 1 by rule 1. But by rule 2, $v - n_p$ must be even, hence $n_p = 1$. Similar conclusion can be made for any $\beta \in [1, 2]$.

Finally, similar rationale can be followed to derive MM updates in closed-form for $\beta = 3/2$. The case $\beta = 2$ is less interesting to tackle with MM updates rules since each subproblem is L-smooth, allowing the use of highly efficient methods such as the Projected Accelerated Gradient Descent algorithm from Nesterov [2013], or the famous HALS algorithm. The derivation of the MM updates in the case $\beta = 0$ with $\ell_2^2$ regularizations on the factors is below in Appendix G.

## G Factor updates for $\beta = 0$ with $\ell_2^2$ regularizations

In this section, we present the steps followed to derive MM update when the data fitting term corresponds to the well-known "Itakura-Saito" divergence, that is for $\beta = 0$, and with $\ell_2^2$ regularization on $\{X_i[i, q]\}_{i \le n}$ for all $q \le r$. Since the update of all factors can be derived in a similar fashion, we only detail the update for $X_1 \in \mathbb{R}_+^{m_1 \times r_1}$, that is for $i = 1$, without loss of generality. The MM update for $X_1$ boils down to tackle the following subproblem:

$$ \underset{\tilde{X}_1 \ge 0}{\operatorname{argmin}} D_0(T|X_1 U) + \mu \sum_{q \le r} \|X_1[i, q]\|_2^2 \qquad (105) $$

where $T = T_{[1]} \in \mathbb{R}^{m_1 \times m_2 \times m_3}$, $U = G_{[1]}(X_2 \boxtimes X_3)^T \in \mathbb{R}_+^{r_1 \times m_2 m_3}$ and $\mu > 0$. We follow the methodology detailed in Section 3.2. In Problem (105), the regularization terms are separable w.r.t. each entry of $X_1[i, q]$ for all $q \le r$. For
the data fitting term, we use the results from Févotte and Idier [2011]. Since $\sum_{q \le r} \|X_1[:q]\|_2^2 = \|X_1\|_F^2$, this leads to solve the following problem:

$$ \underset{X_1 \ge 0}{\operatorname{argmin}} \quad G(X_1 | \tilde{X}_1) + \mu \|X_1\|_F^2 \qquad (106) $$

where $G(X_1|\tilde{X}_1)$ is the separable, tight majorizer for $D_0(T|X_1U)$ at $\tilde{X}_1$ as proposed in [Févotte and Idier, 2011, Lemma 1].

Given that the objective function in Problem (106) is separable with respect to each entry $X_1[i, k]$ of $X_1$ where $1 \le i \le m_1$ and $1 \le k \le r_1$, we can update each entry separately. To find the value of $X_1[i, k]$, we solve

$$ \nabla_{X_1[i,k]} [G(X_1|\tilde{X}_1) + \mu \|X_1\|_F^2] = 0 \qquad (107) $$

and check if the solution is unique and satisfies the nonnegativity constraints. Considering the discussion in Section 3.2 and the findings of [Leplat et al., 2021, Proposition 1], we are aware that this minimizer exists and is unique over the range of $(0, +\infty)$ due to the condition of $\mu > 0$ and the regularization function squared $\ell_2$ norm satisfying Assumption 1 with constants $C = 2$. For now, we put this result aside and simply compute the solution of Equation (107) and then check if it exists and is unique over the positive real line.

From Lemma [Févotte and Idier, 2011, Lemma 1] and Table 1 for $\beta = 0$, one can show that solving (107) reduces to finding $X_1[i, k] > 0$ of the following scalar equation:

$$ -\bar{a}(X_1[i, k])^{-2} + 2\mu X_1[i, k] + \bar{c} = 0 \qquad (108) $$

where $\bar{a} = (\tilde{X}_1[i, k])^2 \sum_j U[k, j] \frac{T[i, j]}{(\bar{T}[i, j])^2} \ge 0$, $\tilde{T} = \tilde{X}_1 U$ and $\bar{c} = \sum_j \frac{U[k, j]}{\bar{T}[i, j]}$ $^{12}$. Posing $x = X_1[i, k]$, Equation (108) can be equivalently rewritten as follows:

$$ x^3 + px^2 + qx + r = 0, \qquad (109) $$

where $p = \frac{\bar{c}}{2\mu} (\ge 0)$, $q = 0$ and $r = -\frac{\bar{a}}{2\mu} (\le 0)$. Equation (109) can be rewritten in the so-called normal form by posing $x = z - p/3$ as follows:

$$ z^3 + az + b = 0, \qquad (110) $$

where $a = \frac{1}{3}(3q-p^2) = -\frac{p^2}{3}$ and $b = \frac{1}{27}(2p^3-9pq+27r) = \frac{1}{27}(2p^3+27r)$. Under the condition $\frac{b^2}{4} + \frac{a^3}{27} > 0$, Equation (110) has one real root and two conjugate imaginary roots. This condition boils down to $\frac{1}{4}(2p^3+27r)^2 - p^6 > 0$. By definition of $p$ (non-negative) and $r$ (non-positive) above, this condition holds if $p > 0$ and $r < -\frac{4p^3}{27}$. The first condition on $p$ is satisfied assuming that $\bar{c}$ is positive. The second condition on $r$ is satisfied by further assuming $\mu > \sqrt{\frac{1}{27} \frac{(\bar{c})^3}{\bar{a}}}$ $^{13}$. Therefore, the positive real root has the following closed-form expression:

$$ z = \sqrt[3]{-\frac{b}{2} + \sqrt{\frac{b^2}{4} + \frac{a^3}{27}}} + \sqrt[3]{-\frac{b}{2} - \sqrt{\frac{b^2}{4} + \frac{a^3}{27}}}, \qquad (111) $$

and hence

$$ X_1[i, k] = z - \frac{\bar{c}}{6\mu}, \qquad (112) $$

With a bit of technicalities, one can show this root is positive, under the assumptions $\bar{c} > 0$ and $\mu > \sqrt{\frac{1}{27} \frac{(\bar{c})^3}{\bar{a}}}$. In matrix form, we obtain the following update

$$ X_1 = Z - \frac{[\bar{C}]}{[6\mu]}, $$

where $Z = [\frac{-B}{2} + [\frac{[B]^2}{4} + \frac{[A]^3}{27}]^{(1/2)}]^{(1/3)} + [\frac{-B}{2} - [\frac{[B]^2}{4} + \frac{[A]^3}{27}]^{(1/2)}]^{(1/3)}$, $A = -\frac{1}{3} [\frac{[\bar{C}]}{[2\mu]}]^2$, $B = \frac{1}{27} (2 [\frac{[\bar{C}]}{[2\mu]}]^3 + 27 [\frac{-\bar{A}}{[2\mu]}])$, $\bar{A} = [\tilde{X}_1]^2 \odot (\frac{[T]}{[\tilde{X}_1 U]^2} U^\top)$ and $\bar{C} = \frac{[E]}{[\tilde{X}_1 U]} U^\top$, with $E$ denoting a matrix full of ones of size $m_1 \times m_2 m_3$.

Note that we skipped in the discussion above the following cases:

<footnote n="12">¹²The solution $X_1[i, k] = 0$ only exists if $\bar{c} = 0$, that is, by definition of $\bar{c}$, when the column $U[k, :]$ is zero.</footnote>

<footnote n="13">¹³In practice, this condition means that $\mu$ should not be chosen too small.</footnote>
1. $\frac{b^2}{4} + \frac{a^3}{27} = 0$, that is when $r = -\frac{4p^3}{27}$, or equivalently $\mu = \sqrt{\frac{1}{27} \frac{(\bar{c})^3}{\bar{a}}}$: here, Equation (110) has one nonnegative real root given by $z = \sqrt{-\frac{\bar{a}}{3}} = \frac{1}{3} \frac{\bar{c}}{\mu}$ (if $b < 0$) and hence $X_1[i, k] = z - \frac{\bar{c}}{6\mu} = \frac{\bar{c}}{6\mu} \ge 0$. If $b > 0$, the root is given by $z = \frac{\bar{c}}{6\mu}$, and then $X_1[i, k] = 0$.

2. $\frac{b^2}{4} + \frac{a^3}{27} < 0$, that is when $\mu < \sqrt{\frac{1}{27} \frac{(\bar{c})^3}{\bar{a}}}$: then Equation (110) has one nonnegative real root given among $z = \frac{1}{3} \frac{\bar{c}}{\mu} \cos(\Phi + \frac{2k\pi}{3})$ with $k = 0, 1, 2$ and where $\cos \Phi = \frac{1}{2} \frac{|2p^3+27r|}{p^3}$.

Those cases are treated in the implementation of the algorithm.

## H Scaling slices and the core with tensor Sinkhorn

At first glance the sNTD does not fit in the HRSI framework. While there is a form of scale ambiguity in NTD, it is not of the form we introduced in Assumption A1 in Section 1.3. The scaling ambiguity works pairwise between the core and each factor, for instance on the first mode,

$$ (\mathcal{G} \times_1 \Lambda) \times_1 (X_1 \Lambda^{-1}) \times_2 X_2 \times_3 X_3 = \mathcal{G} \times_1 X_1 \times_2 X_2 \times_3 X_3 \quad (113) $$

for a diagonal matrix $\Lambda$ with positive entries. This means that the core matrix has scaling ambiguities affecting its rows, columns, and fibers. It then makes little sense to talk about scaling the columns of $\mathcal{G}$ within the HRSI framework.

The fact that the scaling invariances do not fit the HRSI framework does not change the fact that the regularizations will enforce implicitly a particular choice of the parameters $\mathcal{G}$ and $X_1, X_2, X_3$ at optimality. Therefore it should still be valuable to design a balancing procedure to enhance convergence speed, as we have proposed for sNMF and rCPD.

The scalar scale invariance detailed in Section 4 does not encompass all scaling invariances that occur in NTD. Again, in general, for each mode, the core and factor can be rescaled by inserting a diagonal matrix in-between, e.g. a matrix $\Lambda_1$ for the first mode: $\mathcal{G} \times_1 X_1 = (\mathcal{G} \times_1 \Lambda_1) \times_1 (X_1 \Lambda_1^{-1})$, where $\Lambda_1$ has strictly positive diagonal elements. The scaling on each mode affects all the slices of the core, therefore the HRSI framework does not directly apply. Nevertheless, it is possible to derive an algorithmic procedure similar to the Sinkhorn algorithm to scale the factors and the core. We provide detailed derivations for the special case of sNTD, but the proposed method extends to any regularized NTD with homogeneous regularization functions.

To understand the connection between the Sinkhorn algorithm, which aims at scaling a matrix or tensor so that the slices are on the unit simplex Sinkhorn and Knopp [1967], Michael [1970], Lamond and Stewart [1981], and the proposed procedure, consider the following: the results of Section 2 applied to scale invariance on the first mode inform that at optimality, $\beta_{X_1}[q] := 2\|X_1^*[\cdot, q]\|_2^2 = \|\mathcal{G}_{[1]}^*[\cdot, q]\|_1$, where $\mathcal{G}_{[i]}$ is the mode $i$ unfolding of the core tensor $\mathcal{G}$. This observation implies that solutions of sparse NTD must in particular be solutions to a system of equations of the form:

$$ \beta_{X_1}[q] = \|\mathcal{G}_{[1]}[\cdot, q]\|_1 \quad (114) $$

$$ \beta_{X_2}[q] = \|\mathcal{G}_{[2]}[\cdot, q]\|_1 \quad (115) $$

$$ \beta_{X_3}[q] = \|\mathcal{G}_{[3]}[\cdot, q]\|_1 \quad (116) $$

Optimally balancing sparse NTD can therefore be reduced to two operations: first, compute the optimal regularizations $\beta_{X_i}$, then transform the core and the factors so that(114)-(116) holds. Supposing we know the $\beta_{X_i}$ values, computing the scalings of the slices of the core tensor in each mode is exactly what the tensor Sinkhorn algorithm is designed for Csiszar [1975], Sugiyama et al. [2017]. However, since the $\beta_{X_i}$ values are not known, we cannot directly apply a tensor Sinkhorn algorithm. Tensor Sinkhorn usually involves an alternating normalization procedure, which we shall mimic in the following.

We propose a scaling heuristic which, starting from initial diagonal matrices $\Lambda_1, \Lambda_2, \Lambda_3$ set to identity matrices, alternatively balances each mode optimally:

$$ \underset{\Lambda_1 \text{ diagonal}}{\operatorname{argmin}} \quad \|\mathcal{G} \times_1 \Lambda_1 \times_2 \Lambda_2 \times_3 \Lambda_3\|_1 + \|X_1 \Lambda_1^{-1}\|_F^2 \quad (117) $$

$$ \underset{\Lambda_3 \text{ diagonal}}{\operatorname{argmin}} \quad \|\mathcal{G} \times_1 \Lambda_1 \times_2 \Lambda_2 \times_3 \Lambda_3\|_1 + \|X_2 \Lambda_2^{-1}\|_F^2 \quad (118) $$

$$ \underset{\Lambda_3 \text{ diagonal}}{\operatorname{argmin}} \quad \|\mathcal{G} \times_1 \Lambda_1 \times_2 \Lambda_2 \times_3 \Lambda_3\|_1 + \|X_3 \Lambda_3^{-1}\|_F^2. \quad (119) $$

This procedure computes estimates of the $\beta_{X_i}$ marginals iteratively for each mode, and scales the core tensor and the factors accordingly on the fly. It can be seen as a modified tensor Sinkhorn algorithm where the $\beta_{X_i}$ marginals
<figure><img src="image_10.png" /><figcaption>Figure 6: Left: sNMF experiment, center: rNCPD experiment, right: sNTD experiment.</figcaption></figure>

are learned iteratively from the input factors $X_1, X_2, X_3$ and core $\mathcal{G}$ in order to satisfy Equations (114)-(116). Each minimization step boils down to a balancing with two factors in the spirit of what is showcased for sNMF in Section 4. Balancing the first mode yields

$$ \mathcal{G}[q, : ,:] \leftarrow \frac{\beta_{X_1}[q]}{\|\Lambda_2 \mathcal{G}[q, : ,:] \Lambda_3\|_1} \mathcal{G}[q, : ,:] \quad \text{and} \quad X_1[\cdot, : ,:] \leftarrow \sqrt{\frac{\beta_{X_1}[q]}{2\|X_1[\cdot, : ,:]\|_2^2}} X_1[\cdot, : ,:] \quad (120) $$

with $\beta_{X_1} = (\sqrt{2}\|X_1\|_F\|\mathcal{G} \times_2 \Lambda_2 \times_3 \Lambda_3\|_1)^{2/3}$. The other two modes are updated in a similar way. Note that in the alternating scaling algorithm, the scales are not explicitly stored. Rather they are embedded in the core at each iteration because of the scaling update.

As the cost decreases at each iteration and is bounded below, the algorithm is guaranteed to converge. The update rule is well-defined if the factors and the core have non-zero columns/slices. Additionally, since each subproblem has a unique minimizer that is attained, and given that the objective function to minimize in Equation (117) is continuously differentiable, every limit point of this block-coordinate descent algorithm is a critical point (see [Bertsekas, 1999, Proposition 2.7.1]). If the objective function is a proper closed function and is continuous over its domain (not necessarily continuously differentiable), then any limit point of the sequence is a coordinate-wise minimum (see [Beck, 2017, Theorem 14.3]).

While this guarantees that the scales of the factors and core will be improved, unlike for the tensor Sinkhorn, we are not able to guarantee that the proposed heuristic computes an optimal scaling although unreported experiments hint towards this direction. The proposed adaptive tensor Sinkhorn algorithm is summarized for the sNTD problem in Algorithm 3. Finally, we handle the null components and the entries threshold $\epsilon$ for the scaling as explained in Section 3.6.

We did not report the results of the Sinkhorn scaling since they were overall similar to what is obtained with the scalar scaling. The scalar scaling is moreover much easier to implement.

# I Additional experiments: Factor match score evaluation

Figure 6 shows the Factor Match Score (FMS) Acar et al. [2011] for each experiments conducted in Section 4. There is significant variability in the results, but overall balancing does not generally improve FMS. Therefore when scaling-swamps occur, typically the loss is not well minimized but the estimated components are well estimated up to scaling.

# References

* Richard A. Harshman. Foundations of the PARAFAC procedure: Models and conditions for an "explanatory" multimodal factor analysis. *UCLA working papers in phonetics*, 16, 1970.

* Pierre Comon and Christian Jutten. *Handbook of Blind Source Separation: Independent component analysis and applications*. Academic press, 2010.

* Daniel D. Lee and H. Sebastien Seung. Learning the parts of objects by non-negative matrix factorization. *Nature*, 401 (6755):788-791, 1999.

* Nicolas Gillis. *Nonnegative Matrix Factorization*. SIAM, Philadelphia, 2020.

* Morten Mørup, Lars Kai Hansen, and Sidse M Arnfred. Algorithms for sparse nonnegative Tucker decompositions. *Neural computation*, 20(8):2112-2131, 2008.

* Tamara G. Kolda and Brett W. Bader. Tensor decompositions and applications. *SIAM Review*, 51(3):455-500, 2009.
<caption>**Algorithm 3 Balancing algorithm for sNTD**</caption>

**Require:** $X_1, X_2, X_3, \mathcal{G}$, maximal number of iterations, itermax.
1: Initialize $\beta_{X_1}[q], \beta_{X_2}[q]$ and $\beta_{X_3}[q]$ respectively with $\|X_1[\cdot, q]\|_2^2, \|X_2[\cdot, q]\|_2^2, \|X_3[\cdot, q]\|_2^2$ for all $q$.
2: **for** $i = 1$ : itermax **do**
3: &nbsp;&nbsp;&nbsp;&nbsp;**for** $q = 1: r_1$ **do**
4: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $\beta_{X_1}[q] = (\sqrt{2\beta_{X_1}[q]}\|\mathcal{G}[q, \cdot, \cdot]\|_1)^{\frac{2}{3}}$
5: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $\mathcal{G}[q, \cdot, \cdot] = \frac{\beta_{X_1}[q]}{\|\mathcal{G}[q, \cdot, \cdot]\|_1} \mathcal{G}[q, \cdot, \cdot]$ or $\mathcal{G}[q, \cdot, \cdot] = 0$ if $\beta_{X_1}[q] = 0$
6: &nbsp;&nbsp;&nbsp;&nbsp;**end for**
7: &nbsp;&nbsp;&nbsp;&nbsp;**for** $q = 1: r_2$ **do**
8: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $\beta_{X_2}[q] = (\sqrt{2\beta_{X_2}[q]}\|\mathcal{G}[\cdot, q, \cdot]\|_1)^{\frac{2}{3}}$
9: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $\mathcal{G}[\cdot, q, \cdot] = \frac{\beta_{X_2}[q]}{\|\mathcal{G}[\cdot, q, \cdot]\|_1} \mathcal{G}[\cdot, q, \cdot]$ or $\mathcal{G}[\cdot, q, \cdot] = 0$ if $\beta_{X_2}[q] = 0$
10: **end for**
11: **for** $q = 1: r_3$ **do**
12: &nbsp;&nbsp;&nbsp;&nbsp; $\beta_{X_3}[q] = (\sqrt{2\beta_{X_3}[q]}\|\mathcal{G}[\cdot, \cdot, q]\|_1)^{\frac{2}{3}}$
13: &nbsp;&nbsp;&nbsp;&nbsp; $\mathcal{G}[\cdot, \cdot, q] = \frac{\beta_{X_3}[q]}{\|\mathcal{G}[\cdot, \cdot, q]\|_1} \mathcal{G}[\cdot, \cdot, q]$ or $\mathcal{G}[\cdot, \cdot, q] = 0$ if $\beta_{X_3}[q] = 0$
14: **end for**
15: **end for**
16: **for** $q = 1: r_1$ **do**
17: &nbsp;&nbsp;&nbsp;&nbsp; $X_1[\cdot, q] = \sqrt{\frac{\beta_{X_1}[q]}{2\|X_1[\cdot, q]\|_F^2}} X_1[\cdot, q]$ or $X_1[\cdot, q] = 0$ if $\|X_1[\cdot, q]\|_2^2 = 0$
18: **end for**
19: **for** $q = 1: r_2$ **do**
20: &nbsp;&nbsp;&nbsp;&nbsp; $X_2[\cdot, q] = \sqrt{\frac{\beta_{X_2}[q]}{2\|X_2[\cdot, q]\|_F^2}} X_2[\cdot, q]$ or $X_2[\cdot, q] = 0$ if $\|X_2[\cdot, q]\|_2^2 = 0$
21: **end for**
22: **for** $q = 1: r_3$ **do**
23: &nbsp;&nbsp;&nbsp;&nbsp; $X_3[\cdot, q] = \sqrt{\frac{\beta_{X_3}[q]}{2\|X_3[\cdot, q]\|_F^2}} X_3[\cdot, q]$ or $X_3[\cdot, q] = 0$ if $\|X_3[\cdot, q]\|_2^2 = 0$
24: **end for**
25: **return** $X_1, X_2, X_3, \mathcal{G}$

Robert Tibshirani. Regression Shrinkage and Selection via The Lasso: A Retrospective. *Journal of the Royal Statistical Society Series B: Statistical Methodology*, 73(3):273–282, 2011. ISSN 1369-7412, 1467-9868. doi:10.1111/j.1467-9868.2011.00771.x. URL https://academic.oup.com/jrsssb/article/73/3/273/7034363.

Małgorzata Bogdan, Ewout van den Berg, Chiara Sabatti, Weijjie Su, and Emmanuel J. Candès. SLOPE—adaptive variable selection via convex optimization. *The annals of applied statistics*, 9(3):1103–1140, 2015. ISSN 1932-6157. doi:10.1214/15-AOAS842. URL https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4689150/.

Claire Boyer, Antonin Chambolle, Yohann De Castro, Vincent Duval, Frédéric De Gournay, and Pierre Weiss. On Representer Theorems and Convex Regularization. *SIAM Journal on Optimization*, 29(2):1260–1281, 2019. ISSN 1052-6234, 1095-7189. doi:10.1137/18M1200750. URL https://epubs.siam.org/doi/10.1137/18M1200750.

Michael Unser. A Unifying Representer Theorem for Inverse Problems and Machine Learning. *Foundations of Computational Mathematics*, 21(4):941–960, 2021. ISSN 1615-3383. doi:10.1007/s10208-020-09472-x. URL https://doi.org/10.1007/s10208-020-09472-x.

Ryan J. Tibshirani. The lasso problem and uniqueness. *Electronic Journal of Statistics*, 7:1456–1490, 2013. ISSN 1935-7524, 1935-7524. doi:10.1214/13-EJS815. URL https://projecteuclid.org/journals/ electronic-journal-of-statistics/volume-7/issue-none/The-lasso-problem-and-uniqueness/ 10.1214/13-EJS815.full.

Bradley Efron, Trevor Hastie, Iain Johnstone, and Robert Tibshirani. Least angle regression. *The Annals of Statistics*, 32(2):407–499, 2004. ISSN 0090-5364, 2168-8966. doi:10.1214/009053604000000067. URL https://projecteuclid.org/journals/annals-of-statistics/volume-32/issue-2/ Least-angle-regression/10.1214/009053604000000067.full.
Andrzej Cichocki, Rafal Zdunek, and Shun-ichi Amari. *Csiszár's Divergences for Non-negative Matrix Factorization: Family of New Algorithms*, volume 3889. Springer, LNCS-3889, 2006. doi:10.1007/11679363_5.

Cédric Févotte and Jérôme Idier. Algorithms for nonnegative matrix factorization with the $\beta$-divergence. *Neural computation*, 23(9):2421–2456, 2011.

Felix J. Weninger Le Roux, Jonathan and John R. Hershey. Sparse NMF–half-baked or well done? *Mitsubishi Electric Research Labs (MERL), Cambridge, MA, USA, Tech. Rep., no. TR2015-023*, 11:13–15, 2015.

J.M. Bioucas-Dias, A. Plaza, N. Dobigeon, M. Parente, Q. Du, P. Gader, and J. Chanussot. Hyperspectral unmixing overview: Geometrical, statistical, and sparse regression-based approaches. *IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing*, 5(2):354–379, 2012.

Charly Caredda, Jérémy E. Cohen, Laurent Mahieu-Williame, Raphaël Sablong, Michäel Sdika, Jacques Guyotat, and Bruno Montcel. Separable spectral unmixing based on the learning of periodic absorbance changes: Application to functional brain mapping using RGB imaging. In *Translational Biophotonics: Diagnostics and Therapeutics III* (2023), Paper 126271D, page 126271D. Optica Publishing Group, June 2023. doi:10.1117/12.2670506.

Axel Marmoret, Jérémy E. Cohen, Nancy Bertin, and Frédéric Bimbot. Uncovering audio patterns in music with nonnegative Tucker decomposition for structural segmentation. In *ISMIR 2020-21st International Society for Music Information Retrieval*, 2020.

R. Bro. *Multi-way Analysis in the Food Industry: Models, Algorithms, and Applications*. PhD thesis, University of Amsterdam, The Netherlands, 1998.

Marie Roald, Carla Schenker, Vince D. Calhoun, Tülay Adalı, Rasmus Bro, Jeremy E. Cohen, and Evrim Acar. An AO-ADMM approach to constraining PARAFAC2 on all modes. *SIAM Journal on Mathematics of Data Science*, 4 (3):1191–1222, 2022. ISSN 2577-0187. doi:10.1137/21M1450033. URL http://arxiv.org/abs/2110.01278.

N. D. Sidiropoulos, R. Bro, and G. B. Giannakis. Parallel factor analysis in sensor array processing. *IEEE Trans. Sig. Proc.*, 48(8):2377–2388, aug 2000.

Joseph B Kruskal. Three-way arrays: rank and uniqueness of trilinear decompositions, with application to arithmetic complexity and statistics. *Linear algebra and its applications*, 18(2):95–138, 1977.

Andrzej Cichocki, Rafal Zdunek, Anh-Huy Phan, and Shun-ichi Amari. *Nonnegative Matrix and Tensor Factorization*. Wiley, 2009.

Lek-Heng Lim and Pierre Comon. Nonnegative approximations of nonnegative tensors. *Journal of chemometrics*, 23 (7-8):432–441, 2009.

Martin J. Mohlenkamp. The Dynamics of Swamps in the Canonical Tensor Approximation Problem. *SIAM Journal on Applied Dynamical Systems*, 18(3):1293–1333, January 2019. ISSN 1536-0040. doi:10.1137/18M1181389.

Arthur Marmin, José Henrique de Morais Goulart, and Cédric Févotte. Majorization-Minimization for Sparse Nonnegative Matrix Factorization With the $\$/\beta$$\$-Divergence. *IEEE Transactions on Signal Processing*, 71:1435–1447, January 2023. ISSN 1053-587X. doi:10.1109/TSP.2023.3266939.

Nathan Srebro, Jason Rennie, and Tommi Jaakkola. Maximum-Margin Matrix Factorization. In *Advances in Neural Information Processing Systems*, volume 17. MIT Press, 2004. URL https://proceedings.neurips.cc/paper/2004/hash/e0688d13958a19e087e123148555e4b4-Abstract.html.

Ming Yuan and Yi Lin. Model selection and estimation in regression with grouped variables. *Journal of the Royal Statistical Society: Series B (Statistical Methodology)*, 68(1):49–67, 2006. ISSN 1467-9868. doi:10.1111/j.1467-9868.2005.00532.x. URL https://onlinelibrary.wiley.com/doi/abs/10.1111/j.1467-9868.2005.00532.x.

Nicholas D. Sidiropoulos Papalexakis, Evangelos E. and Rasmus Bro. From k-means to higher-way co-clustering: Multilinear decomposition with sparse latent factors. *IEEE transactions on signal processing*, 61(2):493–506, 2013.

Alexis Benichoux, Emmanuel Vincent, and Remi Gribonval. A fundamental pitfall in blind deconvolution with sparse and shift-invariant priors. In *IEEE International Conference on Acoustics, Speech and Signal Processing*, pages 6108–6112, 2013. ISBN 978-1-4799-0356-6. doi:10.1109/ICASSP.2013.6638838. URL http://ieeexplore.ieee.org/document/6638838/.

Morten Mørup and Lars Kai Hansen. Automatic relevance determination for multi-way models. *Journal of Chemometrics*, 23(7-8):352–363, 2009. ISSN 1099-128X. doi:10.1002/cem.1223.

Vincent Y.F. Tan and Cédric Févotte. Automatic Relevance Determination in Nonnegative Matrix Factorization with the /spl beta/-Divergence. *IEEE Transactions on Pattern Analysis and Machine Intelligence*, 35(7):1592–1605, July 2013. ISSN 1939-3539. doi:10.1109/TPAMI.2012.240.
André Uschmajew. Local convergence of the alternating least squares algorithm for canonical tensor approximation. *SIAM J. matrix Analysis*, 33(2):639–652, 2012.

Behnam Neyshabur, Ryota Tomioka, and Nathan Srebro. In search of the real inductive bias: On the role of implicit regularization in deep learning. In *3rd International Conference on Learning Representations, ICLR 2015, San Diego, CA, USA, Workshop Track Proceedings*, 2015. URL http://arxiv.org/abs/1412.6614.

Ryan J Tibshirani. Equivalences between sparse models and neural networks. *Working Notes*. URL https://www.stat.cmu.edu/ryantibs/papers/sparsitynn.pdf, 2021.

Tolga Ergen and Mert Pilanci. Implicit convex regularizers of {cnn} architectures: Convex optimization of two- and three-layer networks in polynomial time. In *International Conference on Learning Representations*, 2021. URL https://openreview.net/forum?id=0N8jUH4JMv6.

Pierre Stock, Benjamin Graham, Rémi Gribonval, and Hervé Jégou. Equi-normalization of neural networks. In *International Conference on Learning Representations*, 2019. URL https://openreview.net/forum?id=r1gEqiC9FX.

Rémi Gribonval, Theo Mary, and Elisa Riccietti. Scaling is all you need: quantization of butterfly matrix products via optimal rank-one quantization. In *29ème Colloque sur le traitement du signal et des images (GRETSI), Aug 2023, Grenoble, France.*, pages pp. 497–500, 2023.

Rahul Parhi and Robert D. Nowak. Deep Learning Meets Sparse Regularization: A signal processing perspective. *IEEE Signal Processing Magazine*, 40(6):63–74, September 2023. ISSN 1558-0792. doi:10.1109/MSP.2023.3286988.

Francis Bach, Rodolphe Jenatton, Julien Mairal, and Guillaume Obozinski. Optimization with Sparsity-Inducing Penalties. *Foundations and Trends® in Machine Learning*, 4(1):1–106, 2011.

C. Févotte, N. Bertin, and J-L. Durrieu. Nonnegative matrix factorization with the Itakura-Saito divergence: With application to music analysis. *Neural computation*, 21(3):793–830, 2009.

Ayanendranath Basu, Ian R Harris, Nils L Hjort, and MC Jones. Robust and efficient estimation by minimising a density power divergence. *Biometrika*, 85(3):549–559, 1998.

Shinto Eguchi and Yutaka Kano. Robustifying maximum likelihood estimation. *Tokyo Institute of Statistical Mathematics*, *Tokyo, Japan, Tech. Rep.*, 2001.

A. Lefèvre. *Méthode d'apprentissage de dictionnaire pour la séparation de sources audio avec un seul capteur*. PhD thesis, Ecole Normale Supérieure de Cachan, 2012.

Valentin Leplat, Nicolas Gillis, and Jérôme Idier. Multiplicative updates for nmf with beta-divergences under disjoint equality constraints. *SIAM Journal on Matrix Analysis and Applications*, 42(2):730–752, 2021. doi:10.1137/20M1377278.

Nicolas Gillis and Francois Glineur. Accelerated multiplicative updates and hierarchical ALS algorithms for nonnegative matrix factorization. *Neural computation*, 24(4):1085–1105, 2012.

Meisam Razaviyayn, Mingyi Hong, and Zhi-Quan Luo. A unified convergence analysis of block successive minimization methods for nonsmooth optimization. *SIAM Journal on Optimization*, 23(2):1126–1153, 2013.

L. T. K. Hien and N. Gillis. Algorithms for nonnegative matrix factorization with the kullback-leibler divergence. *Journal of Scientific Computing*, (87):93, 2021.

T. Takahashi and H. Ryota. Global convergence of modified multiplicative updates for nonnegative matrix factorization. *Computational Optimization and Applications*, (57):417–440, 2014.

Jean Kossaifi, Yannis Panagakis, Anima Anandkumar, and Maja Pantic. Tensorly: Tensor learning in python. *The Journal of Machine Learning Research*, 20(1):925–930, 2019.

Jérémy E Cohen. Régularisation implicite des factorisations de faible rang pénalisées. *GRETSI 2023, Grenoble*, 2023.

Jordan B. L. Smith and Masataka Goto. Nonnegative Tensor Factorization for Source Separation of Loops in Audio. In *2018 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)*, pages 171–175, 2018. doi:10.1109/ICASSP.2018.8461876. URL https://ieeexplore.ieee.org/abstract/document/8461876?casa_token=SzYH8C--TmQAAAAA:8kHygCQIvSDP8U0pWsE9NTwBigBV4_3IJ0hmEq4rSc9reVFZt1LBvkELXeWaLIbgItiqAZnRBN8Z.

Lieven De Lathauwer and Joos Vandewalle. Dimensionality reduction in higher-order signal processing and rank-($R_1, R_2, ..., R_n$) reduction in multilinear algebra. *Linear Algebra and its Applications*, 391:31–55, 2004.

Edgar Rechtschaffen. Real roots of cubics: explicit formula for quasi-solutions. *The Mathematical Gazette*, (524): 268–276, 2008.
Patrik O. Hoyer. Non-negative sparse coding. In *IEEE Workshop on Neural Networks for Signal Processing*, pages 557–565. IEEE, 2002.

Yurii Nesterov. *Introductory lectures on convex optimization: A basic course*. Springer Science & Business Media, 2013.

Richard Dennis Sinkhorn and Paul Joseph Knopp. Concerning nonnegative matrices and doubly stochastic matrices. *Pacific Journal of Mathematics*, 27(2):343–348, 1967.

Bacharach Michael. *Biproportional Matrices and Input-Output Change*. Cambridge University Press, 1970.

B. Lamond and N.F. Stewart. Bregman’s balancing method. *Transportation Research Part B: Methodological*, 15(4): 239–248, 1981. ISSN 0191-2615.

I. Csiszar. *i*-divergence geometry of probability distributions and minimization problems. *The Annals of Probability*, 3 (1):146–158, 1975.

Mahito Sugiyama, Hiroyuki Nakahara, and Koji Tsuda. Tensor balancing on statistical manifold. In Doina Precup and Yee Whye Teh, editors, *Proceedings of the 34th International Conference on Machine Learning*, volume 70 of *Proceedings of Machine Learning Research*, pages 3270–3279. PMLR, 06–11 Aug 2017.

Dimitri P. Bertsekas. *Nonlinear Programming*. Athena Scientific, Belmont, MA, 1999.

Amir Beck. *First-Order Methods in Optimization*. Society for Industrial and Applied Mathematics, Philadelphia, PA, 2017.

Evrim Acar, Daniel M. Dunlavy, Tamara G. Kolda, and Morten Mørup. Scalable tensor factorizations for incomplete data. *Chemometrics and Intelligent Laboratory Systems*, 2011.