# Rapture of the deep: highs and lows of sparsity in a world of depths

Rémi Gribonval, Elisa Riccietti, Quoc-Tung Le, Léon Zheng

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

## To cite this version:

Rémi Gribonval, Elisa Riccietti, Quoc-Tung Le, Léon Zheng. Rapture of the deep: highs and lows of sparsity in a world of depths. IEEE Signal Processing Magazine, 2025, 22 p. hal-04954574v4

HAL Id: hal-04954574

https://inria.hal.science/hal-04954574v4

Submitted on 9 Sep 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
Rapture of the deep:
highs and lows of sparsity in a world of depths

Rémi Gribonval* Elisa Riccietti * Quoc-Tung Le† Léon Zheng‡

September 9, 2025

Abstract

Promoting sparsity in deep networks is a natural way to control their complexity, and a timely endeavor, since practical neural model sizes have grown to unprecedented levels. The lessons from sparsity in linear inverse problems also bear the promise of many other benefits beyond such computational aspects, from statistical significance to explainability. Can these promises be fulfilled? Can we safely leverage the know-how of sparsity-promoting regularizers for inverse problems to harness sparsity in deeper contexts, linear or not? This article surveys the curses and blessings of deep sparsity. After a reminder on the main lessons from inverse problems, we tour a number of results that challenge their immediate deep extensions, both from a mathematical and a computational perspective. In particular, we highlight that $l^1$ regularization does not always lead to sparsity, and that optimization with a prescribed set of allowed nonzero coefficients can be NP-hard. We emphasize the role of rescaling-invariances in these phenomena, and the need to favor structured sparsity to keep sparse network training problems under control, ensure their stability, and actually enable efficient network implementations on GPUs. We finally outline the promises and challenges of a flexible family of so-called Kronecker sparsity structures, which extend the classical butterfly structure and appear in many classical scientific computing applications and that have recently emerged also in deep learning.

Introduction

The notion of sparsity is ubiquitous in signal and image processing. Sparse vectors – i.e., vectors with only a small number of nonzero entries, measured by their $l^0$ pseudo-norm – can be represented with much fewer bits than dense vectors of the same dimension. This motivates the quest for sparse representations to compress signals and images, culminating in the MP3 and JPEG/JPEG2000 compression formats. The success of these approaches is largely due to some sort of a *gift from nature*: luckily, many natural signals and images can be sparsely represented in familiar transformed domains such as the (time)-frequency domain, the wavelet domain, or in learned dictionaries [39], as illustrated on Figure 1.

In many settings, sparsity indeed comes as a *natural objective* for efficient data processing: not only for *data compression*, but also for *efficient computations*, as the existence of fast transforms such as the fast Fourier transform (FFT) or the fast Wavelet transform correspond to somewhat deep sparse matrix products (see Figure 2).

Much beyond this natural role for compression, the last two decades have also seen tremendous advances in the context of linear inverse problems that have uncovered a less obvious but fundamental role of sparsity: when used as a *prior knowledge*, it is a key enabler to allow

<footnote>*ENS de Lyon, CNRS, Inria, Université Claude Bernard Lyon 1, LIP, UMR 5668, 69342, Lyon cedex 07, France</footnote>

<footnote>†Toulouse School of Economics, ANITI – Artificial and Natural Intelligence Toulouse Institute</footnote>

<footnote>‡Huawei, Lagrange Research Center</footnote>
<figure><img src="image_3.png" /><figcaption>Figure 1: Natural signals and their sparse transforms. The waveform of a sound signal is less sparse than its time-frequency representation (heatmap colors: black corresponds to negligible time-frequency coefficients). Similarly, an image is generally not sparse in the pixel-domain, but its wavelet transform is sparser (white corresponds to negligible wavelet coefficients)</figcaption></figure>

the identification of sparse vectors from their low-dimensional projection, using algorithms of bounded complexity with provably good performance. This is notably behind the development of compressive sensing [1].

# 1 Sparsity: main lessons from inverse problems

The vast body of signal processing literature on inverse problems and sparsity (see e.g., the book [10]) has led to a few simple lessons that we may summarize as follows:

* **Lesson 1:** $\ell^1$ minimization induces sparsity, while $\ell^2$ minimization rather induces “flatness”.

* **Lesson 2:** greed is good, (iterative) thresholding allows to identify “significant” coefficients.

* **Lesson 3:** recovering a sparse vector is easy if its support is known.

In slightly more details, the objective of a linear inverse problem is to recover an unknown high-dimensional signal $x \in \mathbb{R}^n$ from a (noisy) linear observation

$$y \approx \mathbf{M}x \in \mathbb{R}^m. \quad (1)$$

This covers many classical problems such as denoising, deblurring, inpainting, or source separation, where a key difficulty is that the dimension $m$ of the observation is (much) lower than the dimension $n$ of the signal. This is where sparsity (together with adequate properties of the matrix $\mathbf{M}$ [10, Chapter 2]) enters the game: if we observe $y:=\mathbf{M}x^*$ where $x^*$ is sufficiently sparse ($\|x^*\|_0$ is small enough), then $x^*$ can be recovered as the unique sparsest solution to the linear problem $y = \mathbf{M}x$. that is to say $x^* = \hat{x}_0$ where for $0 \le p \le \infty$ we denote

$$\hat{x}_p := \arg \min \|x\|_p \text{ s.t. } y = \mathbf{M}x. \quad (2)$$

While (2) with $p = 0$ is computationally intractable, the problem becomes convex when $p = 1$ and it is endowed with efficient algorithms to obtain estimates reaching the minimum $\ell^1$ norm up to a tolerance<fnref n="1" />. Moreover, this $\ell^1$ minimization problem also leads to the guarantee that $\hat{x}_1 = x^*$, under slightly stronger assumptions on the sparsity of $x^*$ and the properties of $\mathbf{M}$ [10, Chapter 4] than those ensuring $\hat{x}_0 = x^*$. This is then even stable to the presence of noise and numerical inaccuracies.

Irrespective of the sparsity of $x_0$ or of the properties of $\mathbf{M}$, this $\ell^1$ optimization problem is also known to promote sparsity in the sense that its set of minimizers always contains at least

<footnote n="1">¹Accurately estimating a minimizer can however be much harder, see e.g. [5].</footnote>
<figure><img src="image_4.png" /><figcaption>Figure 2: Usual Fast Transforms as sparse matrix products. The $n \times n$ matrix of the Discrete Fourier Transform (DFT) is a product of $\log_2 n$ factors with the so-called butterfly structure, each with two nonzero entries per row and per column, leading to the complexity $O(n \log_2 n)$ of the FFT. The matrix of the Discrete Wavelet Transform (DWT) is the product of sparse factors corresponding to subsampled high-pass and low-pass filtering, leading to the $O(n)$ complexity of the DWT. In each case the number of flops is bounded by $\sum_{\ell} \|\mathbf{W}_{\ell}\|_0$ where $\mathbf{W}_{\ell}$ is the $\ell$-th factor.</figcaption></figure>

one vector $x$ such that $\|x\|_0 \le m$. In contrast, the $\ell^2$ optimization problem (2) with $p = 2$ (or its relaxed versions corresponding to ridge regression / Tikhonov regularization) have solutions that depend linearly on $y$ and are typically non sparse. This is the essence of **Lesson 1**.

A convenient notion to think about the sparsity of a vector $x$ or a matrix $X$ is its support

$$ \mathrm{supp}(x) := \{x_i \neq 0\}, \quad \mathrm{supp}(X) := \{(i, j) : X_{i,j} \neq 0\}. \qquad (3) $$

The size of this support indeed measures sparsity: $\#\mathrm{supp}(x) = \|x\|_0$ and $\#\mathrm{supp}(X) = \|X\|_0$. Interestingly, by coupling (partial) knowledge on the support of an unknown $x^*$ with the observation of $y = \mathbf{M}x^*$ we can easily recover $x^*$: if we know a set of indices $T \subseteq \llbracket 1, n \rrbracket := \{1, \dots, n\}$ such that $\mathrm{supp}(x^*) \subset T$ then (assuming linear independence of the columns of $\mathbf{M}$ indexed by $T$, i.e., injectivity of the matrix $\mathbf{M}_T = \mathbf{M}[\cdot, T]$) simple linear algebra guarantees that the entries of $x^*$ on $T$ can be computed via a pseudo-inverse: $x_T^* = \mathbf{M}_T^\dagger y$, while all other entries of $x^*$ are zero. Even in the noisy case, while sparse recovery is known to be NP-hard [28], finding the optimal coefficients for such a candidate support $T$

$$ \min_x \|y - \mathbf{M}x\|_2^2 \quad \text{s.t.} \quad \mathrm{supp}(x) \subseteq T \qquad (4) $$

is a simple least squares problem solved by $x_T = \mathbf{M}_T^\dagger y$ and $x_{T^c} = 0$. This yields **Lesson 3**.

Algorithmically speaking, all state-of-the-art sparse recovery algorithms rely on the magnitude of coefficients or inner products to detect “significant” ones, either via greedy approaches or via thresholding steps: they iteratively refine an estimate of the support of $x^*$ on the basis of how large are the coefficients of some correlation vector. This is notably the case in iterative thresholding algorithms addressing the $\ell^1$-penalized least-square regression variant of problem (2) with $p = 1$, as well as in greedy algorithms such as Orthogonal Matching Pursuit or Hard Thresholding Pursuit [10]. This is **Lesson 2**.

## 2 Surprises and pitfalls of shallow and deep sparsity

Arguably, the most classical model of deep neural network is the multilayer perceptron (MLP) with $L \ge 2$ affine layers parameterized by $L$ weight matrices $\mathbf{W}_1, \dots, \mathbf{W}_L$ (as well as bias vectors that we will omit for simplicity of exposition). For example, with the rectified linear unit (ReLU) activation function $\rho(t) = \mathrm{ReLU}(t) := \max(t, 0)$ and its entrywise extension $\rho(x) = (\rho(x_i))_{i=1}^n$ for vectors, $f_\theta = \mathbf{W}_1 \circ \rho \circ \dots \circ \rho \circ \mathbf{W}_L$ where $\theta := (\mathbf{W}_\ell)_{\ell=1}^L$. An even simpler related model appears when $\rho$ is the identity, leading to a simple multilayer linear model: $f_\theta(x) = \mathbf{A}_\theta x$ where $\mathbf{A}_\theta = \mathbf{W}_1 \dots \mathbf{W}_L$.
Can sparsity also be leveraged in the context of deep learning, where the need to control the complexity/performance tradeoffs of large neural models is becoming more evident everyday? Intuitively, a natural objective to reduce the memory required to store the parameters $\theta$ of a network and/or the computational cost of computing the network output $f_{\theta}(x)$ is to constrain the weight matrices $W_{\ell}$ of each layer to be sparse, as measured by their matrix $\ell^0$ (pseudo)norm $\|W\|_0$ which counts their number of nonzero entries (see Figure 3). As illustrated in this section the road to mathematically harness deep sparsity bears a few surprises, as essentially **each of the lessons from linear inverse problems breaks down in a multilayer context**, whether linear or not.

## New lesson 1 Deep $\ell^p$ regularization does not promote solutions with the expected properties

As seen in the lessons from linear inverse problems, the minimization of the $\ell^1$ norm promotes sparse solutions whereas minimizing the $\ell^2$ norm induces “flat” ones. This knowledge, however, does not hold anymore when the relation between the input and output is no longer linear, i.e., $Mx = y$.

<figure><img src="image_5.png" /><figcaption>Figure 3: Deep neural models involve several dense linear layers associated with dense matrices $W_{\ell}$. Their sparsified version corresponds to sparse matrices $W_{\ell}$, reminiscent of the factorized structure of fast transforms (cf Figure 2) but with much less structure in the sparsity pattern. Performing such a sparsification aims to achieve networks that can be encoded with *fewer bits*, and implemented with *fewer flops*, but are they up to their promises?</figcaption></figure>

As a warmup, consider the *sparse blind deconvolution problem*, where one seeks to recover an unknown filter $h_0$ and an unknown signal $s_0$, both assumed to be sparse, from the only observation of their convolution $y = h_0 \star s_0$. Remember that **Lesson 1** from inverse problems [10, 34] indicates that under appropriate assumptions on the filter $h_0$ (or equivalently on the matrix $M$ such that $Ms = h_0 \star s$), the sparsity of $s_0$ allows to recover it as the solution of an $\ell^1$ minimization problem

$$ \min_{s} \|s\|_{1} \text{ s.t. } y = h_{0} \star s, $$

and vice-versa $h_0$ can be guaranteed to be the minimum $\ell^1$ norm solution to a similar problem. However, it has been shown that addressing the joint optimization problem

$$ \min_{h,s} \|h\|_1 + \|s\|_1 \text{ s.t. } y = h \star s, \qquad (5) $$

*totally fails* to recover $h_0$ and $s_0$: it admits a *global optimum* at $h^* := \lambda\delta_0$ (the Dirac delta at zero) and $s^* := y/\lambda$ for some $\lambda \neq 0$ [2] (and another one at $h^* := y/\lambda$ and $s^* := \lambda\delta_0$).

In the same vein, but closer to deep learning consider the popular **weight decay regularization** approach during neural network learning, which simply corresponds to $\ell^2$ regularization on the parameters. To be concrete, let us focus specifically on the simple example of least-squares regression with a ReLU neural network containing one hidden layer and no bias. Each neuron $j$ is parameterized by a vector $v_j$ of input weights and a vector $u_j$ of outgoing weights, and implements a function $f_{\theta}(x) := \sum_j u_j \text{ReLU}(v_j^\top x)$ (where we gather all parameters
into $\theta = (u_j, v_j)_j$). Training such a network on (a batch of) samples $(x_i, y_i)_{i=1}^N$ with weight decay corresponds to gradient descent on the cost function

$$ \min_{\theta} \frac{1}{N} \sum_{i=1}^{N} (y_i - f_{\theta}(x_i))^2 + \lambda \sum_{j} (\|u_j\|_2^2 + \|v_j\|_2^2). \quad (6) $$

While **Lesson 1** would suggest that the $\ell^2$ regularization term promotes “flat” solutions rather than sparse ones, it has been documented that instead it tends to promote sparse/low-rank solutions [30, 31], and that every optimal parameter vector satisfies a so-called **Neural Balance** property $\|u_j\|_2 = \|v_j\|_2$.

Behind the scene of these two examples, an intrinsic scaling ambiguity plays a key role: for any positive scalar we have the identity $h \star s = (h/\lambda) \star (\lambda s)$. This is due to the bilinearity of the model $y = h \star s$, and this implies that the apparent regularization $\|h\|_1 + \|s\|_1$ in (5) actually corresponds to the implicit regularization

$$ R(h, s) := \min_{\lambda > 0} \|h/\lambda\|_1 + \|\lambda s\|_1 \propto \sqrt{\|h\|_1 \|s\|_1}. $$

### [Box 1] Rescaling equivalence

Consider a single neuron with no bias. In this simple case we have $f_{\theta}(x) = v \max(0, \langle u, x \rangle)$ with $\theta = (u, v)$. From the positive-homogeneity of the ReLU function, it is thus straightforward that for any $\lambda > 0$, the “rescaled” parameter $\tilde{\theta} = (\lambda u, \frac{v}{\lambda})$ implements the same function. A similar rescaling-invariance property holds for more general networks [14], leading to a notion of **rescaling-equivalent parameters**: $\theta' \sim \theta$ if $\theta' = D\theta$ where D belongs to the group of diagonal matrices generated by elementary rescaling operations. Rescaling-equivalence ($\theta' \sim \theta$) implies $f_{\theta'} = f_{\theta}$, but the converse is not always true [37].

The phenomenon illustrated in the case of blind deconvolution then simply follows from the inequality $\|y\|_1 \|\delta_0\|_1 = \|y\|_1 \le \|h\|_1 \|s\|_1$ that is valid [2] for any filter $h$ and source $s$ such that $y = h \star s$. In the case of weight-decay regularization, we similarly have $(u/\lambda)\text{ReLU}(\lambda v^\top x)$ for any positive $\lambda$ and any $u, v, x$ (see Box 1). This is due to the fact that the ReLU function, $\text{ReLU}(t) := t_+ = \max(t, 0)$, is positively homogeneous: it satisfies $\text{ReLU}(\lambda t) = \lambda \text{ReLU}(t)$ for every non-negative $\lambda$. As a consequence, the apparent regularization $\sum_j (\|u_j\|_2^2 + \|v_j\|_2^2)$ in (6) actually corresponds to the implicit regularization

$$ R(\theta) := \sum_j \|u_j v_j^\top\|_2 $$

reminiscent of Group Lasso (we refer to the recent survey [30] for more details).

## New lesson 2 Deep greed can be bad

Another consequence of rescaling invariances in ReLU networks is the lack of guarantees for greedy approaches and thresholding to select “significant” coefficients.

Let us illustrate this in a network compression scenario with two players:

* Alice openly releases a large-scale trained network, with parameters $\theta$;

* Bob wants to obtain a compressed network from it via pruning (i.e., replacing $\theta$ by an approximation $\hat{\theta}$ keeping a few entries of $\theta$ and setting the other entries to zero), cf Figure 3.

The first step of the popular (iterative) magnitude pruning (IMP) [19] approach is to keep a given proportion of the coefficients with the largest magnitude. Due to rescaling invariances, a malicious player (Malice) could replace $\theta$ by a rescaling-equivalent parameter $\theta'$ (see Box 1) corresponding to a network with exactly the same behavior: $f_{\theta'} = f_{\theta}$. The performance of
IMP can indeed be severely degraded when applied on $\theta'$ rather than $\theta$ [14], indicating that the magnitude of coefficients is not a direct indicator of their significance, unlike suggested by **Lesson 2**.

## New lesson 3 Support identification is not the only challenge

Most sparse recovery algorithms for linear inverse problems somehow combine two ingredients: the exploration of different possible supports for the candidate solution, and the optimization of the coefficients for a given support. This is, for example, explicitly done in greedy algorithms such as orthogonal matching pursuit or hard thresholding pursuit [10,27], which iteratively explore the support and compute the optimal coefficients for each new support. Despite the NP-hardness of sparse recovery, finding the optimal coefficients in (4) for a given candidate support $T$ is simple, as we discussed in **Lesson 3**. Thus, the main computational difficulty in classical sparse recovery is to find the best support.

In contrast, in the multilayer context, the problem remains challenging even with known support. Let us focus on sparse matrix factorization, where the goal is to find sparse factors $\mathbf{X}$ and $\mathbf{Y}$ such that $\mathbf{A} \approx \mathbf{X}\mathbf{Y}$. Assume that we know two subsets of matrix indices, $I$ and $J$, that contain the supports of $\mathbf{X}$ and $\mathbf{Y}$. We then seek a solution to the following fixed support matrix factorization (FSMF) problem:

$$ \min_{\mathbf{X},\mathbf{Y}} \|\mathbf{A} - \mathbf{X}\mathbf{Y}\|_F^2, \quad \text{s.t.} \quad \mathrm{supp}(\mathbf{X}) \subseteq I, \mathrm{supp}(\mathbf{Y}) \subseteq J. \qquad (\text{FSMF}) $$

This seemingly simple problem is severely more difficult that it may look [22]:

* **Hardness:** it is NP-hard to approximate$^2$, in sharp contrast to **Lesson 3**; this fact should not however be surprising as it is in-line with previously known NP-hardness results for variants of matrix factorization problems such as non-negative matrix factorization [40] and matrix completion [12].

* **Instability:** it may not even admit an optimal solution, leading to diverging optimization schemes, as illustrated and discussed in Boxes 2 and 3.

It is noteworthy that sparsity might not be the only cause of the absence of a global optimum. A perhaps better-known example (but not related to sparse models) of the non-existence of optimal solutions is low-rank tensor approximation. For tensors of order at least three, the best rank-$r$ ($r \ge 2$) tensor approximation problem does not always admit a minimizer [35]. Moreover, the set of tensors whose instances do not admit a minimizer can have positive Lebesgue measure [35]. Other examples of this phenomenon are ReLU neural networks training [25] and its sparse counterpart [23].

## [Box 2] A minimalistic instance of Problem (FSMF) without an optimum

For certain choices of target matrix $\mathbf{A}$ and of support constraints $(I, J)$, the problem (FSMF) does not admit a global optimum. We illustrate this with a minimalistic example related to LU factorization. Consider $\mathbf{A} = \begin{bmatrix} 0 & 1 \\ 1 & 0 \end{bmatrix}$, $I = \begin{bmatrix} 1 & 0 \\ 1 & 1 \end{bmatrix}$, $J = \begin{bmatrix} 1 & 1 \\ 0 & 1 \end{bmatrix}$. Denoting $\mathbf{X}_k = \begin{bmatrix} 1/k & 0 \\ 0 & -k \end{bmatrix}$, $\mathbf{Y}_k = \begin{bmatrix} 1/k & k \\ 0 & k \end{bmatrix}$ for each $k \ge 1$, we have $\lim_{k \to \infty} \|\mathbf{A} - \mathbf{X}_k \mathbf{Y}_k\|_F = \lim_{k \to \infty} \frac{1}{k^2} = 0$ hence $\mathbf{A}$ is approximated arbitrarily well by factors satisfying the support constraints $I, J$. Yet, there is no pair $(\mathbf{X}, \mathbf{Y})$ such that $\mathbf{A} = \mathbf{X}\mathbf{Y}$ and $\mathrm{supp}(\mathbf{X}) \subseteq I, \mathrm{supp}(\mathbf{Y}) \subseteq J$. Indeed, any such couple would need to satisfy: $\mathbf{X}_{1,1}\mathbf{Y}_{1,2} = 1, \mathbf{X}_{2,1}\mathbf{Y}_{1,1} = 1, \mathbf{X}_{1,1}\mathbf{Y}_{1,1} = 0$ (where $\mathbf{X}_{i,j}$ is the coefficient at the index $(i, j)$ of the matrix $\mathbf{X}$). However, the third equation implies that either $\mathbf{X}_{1,1} = 0$ or $\mathbf{Y}_{1,1} = 0$, which makes either $\mathbf{X}_{1,1}\mathbf{Y}_{1,2} = 0$ or $\mathbf{X}_{2,1}\mathbf{Y}_{1,1} = 0$. This leads to a contradiction.

<footnote>$^2$There exists $\epsilon > 0$ such that finding $\mathbf{X}, \mathbf{Y}$ reaching the optimal value of (FSMF) up to an $\epsilon$ margin is NP hard [22, Th. 2.4].</footnote>
(Brand) new lesson Sparsity is not equivalent to efficiency (in the world of depths)

Last, but not least, we highlight that in the context of deep learning where GPUs are a major hardware tool, naive sparsification of dense network layers does not bring the expected computational benefits [20, Section 7]. This is due to the irregular memory layout of unstructured sparse matrices, which forbids to leverage the processing power of GPUs due to high costs of memory accesses to non-contiguous data. Even with structured sparse matrices, performance might still suffer without dedicated implementations because the additional overhead of memory operations (e.g., data transfer and rewriting operations, both of which are typically higher for sparse models) may offset the benefits of existing efficient implementations for general matrix multiplication on the GPU [15].

In summary, the major lessons of sparse linear inverse problems do not hold anymore in a deeper context, even in the seemingly simplistic bilinear case, for at least two reasons: rescaling invariances of ReLU networks, which require revisiting the implicit regularizers behind apparently explicit regularization schemes; and perhaps more importantly, unfavorable properties of unstructured deep sparsity, which generate key theoretical and practical bottlenecks (from the hardness and instability of training problems to the inefficiency of GPU implementations).

## 3 Rescaling-invariances and non-closedness: how to deal with them

It may seem despairing that most lessons about sparsity learned from inverse problems break down with depth. Yet, having diagnosed underlying causes, one can envision cures to harness the main underlying phenomena: rescaling-invariances and non-closedness. We briefly review possible approaches here.

[Box 3] Spurious valleys, the absence of a global optimum, and non-closedness

The unpleasant phenomenon of the absence of global minimum from Box 2 corresponds to the presence of a “valley at infinity” of the optimization problem, so that any optimization algorithm to address it will either lead to diverging iterates or to a sub-optimal point. As an illustration, consider a sliced view on the landscape of the optimization problem (FSMF) via the function

$$g(\delta) := \min \|\mathbf{A} - \mathbf{X}\mathbf{Y}\|_F^2 \text{ s.t. } \mathrm{supp}(\mathbf{X}) \subseteq I, \mathrm{supp}(\mathbf{Y}) \subseteq J, \mathbf{X}_{2,2}\mathbf{Y}_{2,2} = \delta$$

with $(\mathbf{A}, I, J)$ from Box 2.

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

Note that $\inf_{\delta} g(\delta) = 0$ is also the infimum of (FSMF) and it can be approached only if $\delta = \mathbf{X}_{2,2}\mathbf{Y}_{2,2}$ diverges to infinity. Mathematically, the root of this phenomenon is the non-closedness of the set $\mathcal{B}_{I,J} := \{\mathbf{X}\mathbf{Y} | \mathrm{supp}(\mathbf{X}) \subseteq I, \mathrm{supp}(\mathbf{Y}) \subseteq J\}$. Indeed, one can reformu-
late (FSMF) as

$$ \min_{\mathbf{B}} \|\mathbf{A} - \mathbf{B}\|_F^2 \quad \text{such that} \quad \mathbf{B} \in \mathcal{B}_{I,J}. \qquad (7) $$

By standard arguments, (7) (and also (FSMF)) admits an optimum for every target matrix $\mathbf{A}$ if and only if the set $\mathcal{B}_{I,J}$ is closed (with the usual topology).

In the example of Box 2, the choice $I = \begin{bmatrix} 1 & 0 \\ 1 & 1 \end{bmatrix}$, $J = \begin{bmatrix} 1 & 1 \\ 0 & 1 \end{bmatrix}$ corresponds to the supports of lower/upper triangular matrices, and (FSMF) is equivalent to the LU decomposition of a square matrix. It is then known that $\mathcal{B}_{I,J}$ is dense in $\mathbb{R}^{2 \times 2}$ but $\mathcal{B}_{I,J} \neq \mathbb{R}^{2 \times 2}$. Hence, $\mathcal{B}_{I,J}$ is not closed.

We emphasize that matrix factorization without sparsity constraint or with some other types of constraints (e.g., non-negativity) does not suffer from this non-closedness. Indeed, a problem without sparsity constraint corresponds to index sets $I$ and $J$ both containing all possible indices. In such a case $\mathcal{B}_{I,J}$ is a set of low-rank matrices, which is closed. As a consequence, in the matrix factorization problem, whether there always exist an optimum solution depends on the sparsity (support) constraints.

## 3.1 Explicit handling of implicit regularizers

The relation between apparently explicit regularizers and actual implicit regularizers in the presence of rescaling invariances has a long history in signal processing, optimization and machine learning [16,38]. This has notably led to approaches that directly exploit this property to ease optimization. For example, the minimization of a criterion penalized by an $\ell^1$ regularizer $\lambda \|x\|_1$ leads to a non-smooth optimization problem, introducing surrogate variables $u, v$ such that the unknown $x$ is the Hadamard product $x = u \odot v$ (entrywise product) allows to recast the problem with a smooth $\ell^2$ penalty $\lambda (\|u\|_2^2 + \|v\|_2^2)$ [32]. Vice-versa, when starting from a problem known to satisfy rescaling-invariances, the explicit computation of optimal rescaling parameters during gradient descent can be beneficial [4,36].

## 3.2 Capturing rescaling-invariances with the path-lifting $\Phi$

The rescaling-invariance illustrated in Box 1 for single-hidden-layer ReLU networks still holds for deeper ones: with all parameters (weights and biases) gathered in a parameter vector $\theta$, it is always possible to replace the parameters $\theta$ of any ReLU neural network by rescaling-equivalent ones $\theta'$ that lead to a network with the exact same realization: $f_{\theta}(x) = f_{\theta'}(x)$ for all $x$. These rescaling-invariances can be captured by a suitable tool: the path-lifting (see Box 4), a finite-dimensional polynomial representation $\Phi(\theta)$ of the parameters of ReLU networks, which is invariant under such rescalings [13] as opposed to the networks parameters themselves. First introduced for “pure ReLU” DAG networks [29], this tool was recently extended to a very general “DAG ReLU” network model with max-pool, etc. It has lead to state-of-the-art generalization bounds [13, 29], regularizers for network training [21], network identifiability results [3, 37], or rescaling-invariant Lipschitz bounds on $f_{\theta}$ [13]. Another recently established property [14] is a rescaling invariant Lipschitz bound for the map $\theta \mapsto f_{\theta}$

$$ \|f_{\theta}(x) - f_{\theta'}(x)\|_{1} \leq \max(\|x\|_{\infty}, 1) \|\Phi(\theta) - \Phi(\theta')\|_{1}. \qquad (8) $$

This has led to a pruning algorithm that is both effective and invariant under symmetries. Instead of taking pruning decisions based on a criterion in the parameter space, e.g., based on the magnitude of the weights, it prunes based on the $\ell^1$ norm of the path-lifting, using a criterion that is efficiently computable using one forward pass and autodifferentiation. This pruning matches the accuracy of magnitude pruning when applied to ResNets trained on Imagenet in the lottery ticket context [11], while being rescaling-invariant [14] and thus allowing to avoid the malicious attacks of Malice from our New lesson 2.
[Box 4] Path-lifting

The path-lifting $\Phi(\theta) \in \mathbb{R}^P$ is a vector indexed by the set $\mathcal{P}$ of paths of the network, i.e., a sequence of connected nodes starting at some neuron and ending at an output neuron. For instance, in the simple neural network below with no bias $p = u \to v \to w$ is a path with $u$ an input neuron, $v$ a hidden neuron, and $w$ an output neuron.

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

The coordinate of $\Phi(\theta)$ associated with a path is the product of the weights along this path, ignoring the non-linearities: here for example $\Phi_p(\theta) = \theta^{u \to v} \theta^{v \to w}$. A twin vector to the path-lifting is the binary path-activation vector $a(\theta, x) \in \{0, 1\}^P$ which contains the information about nonlinearities: $a_p(\theta, x) = 1$ if and only if all neurons along path $p$ are activated for the input vector $x$. The definition of these two quantities is valid for so-called DAG-ReLU neural networks, with a directed acyclic graph (DAG) structure, that can mix ReLU and max-pooling activations and include biases [13].

Importantly, both $\Phi(\theta)$ and $a(\theta, x)$ are rescaling-invariant: if $\theta' \sim \theta$ then $\Phi(\theta') = \Phi(\theta)$ and $a(\theta', x) = a(\theta, x)$ for every $x$ [13]. Another key property is that the output of a scalar-valued DAG ReLU network $f_\theta: \mathbb{R}^n \to \mathbb{R}$ can be expressed as:

$$f_{\theta}(x) = \langle \Phi(\theta) \odot a(\theta, x), \mathbf{A}(\frac{x}{1}) \rangle \quad (9)$$

(a similar simple formula holds for vector-valued networks [13, Theorem A.1]) where the binary matrix $\mathbf{A}$, of size $|\mathcal{P}| \times (n+1)$, encodes from which neuron each path starts: $\mathbf{A}_{p,u} := 1$ if and only if $p$ starts at the input neuron $u$ (by convention $u = n+1$ corresponds to a “bias” neuron).

Last, but not least, even though $\Phi(\theta)$ is a vector of combinatorial dimension (the number of paths in the network) whose actual computation would be intractable, **many important quantities such as its $\ell^p$ norms can be efficiently computed**, typically using a single forward pass on a slight modification of the underlying network [13].

## 3.3 Avoiding instability by ensuring closedness

To avoid the cruel dilemma of having to choose between convergence to a local minimum and divergence to infinity of gradient descent in deep or shallow optimization problems such as (FSMF), a standard approach is regularization via the addition of a small coercive regularization term (e.g., any norm on the parameters). While this approach is natural, it imposes to choose tradeoffs between approximation quality and weight magnitudes, leading not only to the difficult problem of calibrating the penalty level but also possibly to the absence of any satisfying tradeoff [23].

However, as the root for the instability of (FSMF) is diagnosed as the non-closedness of the set $\mathcal{B}_{I,J}$ of products of support-constrained factors (see Box 3), another perhaps more satisfying alternative is to *restrict the search to “structured” support constraints* ensuring the closedness of $\mathcal{B}_{I,J}$. Related ideas already appear, e.g., in tensor approximation problems [35], where changing the notion of tensor rank to *border rank* ensures that the low-rank tensor approximation problem always admits minimizers.
### 3.3.1 Closedness conditions in deep sparse matrix factorization

How to characterize “good supports” (I, J) for Problem (FSMF), i.e., pairs of support constraints such that (FSMF) admits minimizers for every matrix A? This question can be extended to the multi-factor setting, i.e., the case in which there are L > 2 factors with support constraints (I₁, ..., I_L). In the vein of the previous analysis, finding “good supports” (I₁, ..., I_L) amounts to answering to the following question:

**Question** (Closedness decision problem). Given $\mathcal{I} := (I_1, \dots, I_L)$ a sequence of binary matrices representing support constraints, decide whether the following set of matrices is closed or not:

$$ \mathcal{B}_{\mathcal{I}} := \{ \mathbf{W}_1 \dots \mathbf{W}_L \mid \mathrm{supp}(\mathbf{W}_{\ell}) \subseteq I_{\ell}, \ell = 1, \dots, L \}. $$

There is indeed an algorithmic answer to this question (that is thus *decidable*) using the Quantifier Elimination Algorithm, a classical tool from real algebraic geometry. However, as discussed in [23], this remains unsatisfactory as: a) the algorithm complexity is doubly exponential, e.g., $O(2^{2^n})$ (where n grows linearly with the dimension of the product $\mathbf{W}_1 \dots \mathbf{W}_L$ and the cardinality of the support constraints $I_\ell$); and b) the algorithm does not provide any insight on how to generate good supports.

An alternative is to explicitly exhibit families of good supports. We discuss in Box 5 such a construction, given in [22], which is also endowed with a polynomial algorithm to find the optimal factors.

To generalize this construction one can express the product $XY$ as a sum of rank-one matrices:

$$ \mathbf{X}\mathbf{Y} = \sum_{i=1}^{r} x_i y_i^{\top} = \sum_{i=1}^{r} \mathbf{C}_i \qquad (10) $$

where $\mathbf{X} = [x_1 \dots x_r] \in \mathbb{R}^{m \times r}$, $\mathbf{Y} = \begin{bmatrix} y_1^{\top} \\ \vdots \\ y_r^{\top} \end{bmatrix} \in \mathbb{R}^{r \times n}$ and $\mathbf{C}_i = x_i y_i^{\top}$ has rank at most one. Observe that $\mathbf{X}$ and $\mathbf{Y}$ satisfy support constraints $(I, J)$, if and only if their columns (resp rows) satisfy corresponding constraints $\mathrm{supp}(x_i) \subseteq I_i$ and $\mathrm{supp}(y_i) \subseteq J_i$, so that $\mathrm{supp}(\mathbf{C}_i) \subseteq I_i \times J_i$ where $I_i$ (resp. $J_i$) is the set specifying the support constraint for the i th column of $\mathbf{X}$ (resp. row of $\mathbf{Y}$) according to $(I, J)$.

### [Box 5] Example of support constraints with closedness and tractability properties

We describe an example of $(I, J)$ for which (FSMF) is not only *guaranteed to admit minimizers for every choice of A*, but also *endowed with an efficient way to compute these minimizers*. Consider a pair of support constraints $(I, J)$ and a matrix $\mathbf{A} \in \mathbb{R}^{4 \times 4}$ written in block form as follows:

$$ I = \begin{pmatrix} 1 & 1 & 0 & 0 \\ 1 & 1 & 0 & 0 \\ 0 & 0 & 1 & 1 \\ 0 & 0 & 1 & 1 \end{pmatrix}, \quad J = \begin{pmatrix} 1 & 1 & 0 & 0 \\ 0 & 0 & 1 & 1 \\ 1 & 1 & 0 & 0 \\ 0 & 0 & 1 & 1 \end{pmatrix}, \quad \mathbf{A} = \begin{pmatrix} \mathbf{A}_1 & \mathbf{A}_2 \\ \mathbf{A}_3 & \mathbf{A}_4 \end{pmatrix} $$

where $\mathbf{A}_i \in \mathbb{R}^{2 \times 2}$. Matrices $\mathbf{X}, \mathbf{Y} \in \mathbb{R}^{4 \times 4}$ satisfying the support constraints write as

$$ \mathbf{X} = \left( \begin{array}{c|c} x_1 & x_2 \\ \hline \mathbf{0} & \mathbf{0} \end{array} \right), \quad \mathbf{Y} = \left( \begin{array}{c|c} y_1^{\top} & \mathbf{0} \\ \hline \mathbf{0} & y_2^{\top} \end{array} \right) $$

with $x_i, y_i \in \mathbb{R}^2$ so that the loss function $\|\mathbf{A} - \mathbf{X}\mathbf{Y}\|_F^2$ can be decomposed as:

$$ \| \mathbf{A}_1 - x_1 y_1^{\top} \|_F^2 + \| \mathbf{A}_2 - x_2 y_2^{\top} \|_F^2 + \| \mathbf{A}_3 - x_3 y_3^{\top} \|_F^2 + \| \mathbf{A}_4 - x_4 y_4^{\top} \|_F^2 $$

The optimization thus decouples into four independent best rank-one approximation problems solved via a (truncated) Singular Value Decomposition (SVD) on each block $\mathbf{A}_i$.
When the subsets $I_i \times J_i$ are **pairwise disjoint or identical**, the instance of (FSMF) is tractable by truncated SVDs on blocks extracted from **A**. This two-factor construction serves as a building block for another structured family of deep sparse factorizations (**the butterfly sparsity**), endowed with closedness guarantees and algorithms with provable performance, as discussed in the upcoming dedicated section.

### 3.3.2 Sparse neural network training: closed or not?

The closedness decision problem (Section 3.3.1) naturally extends to the problem of sparse neural network training, i.e.,

$$ \begin{array}{ll} \text{Minimize} & \displaystyle\sum_{i=1}^{N} (\mathbf{W}_L \text{ReLU}(\mathbf{W}_{L-1} \dots \text{ReLU}(\mathbf{W}_1 x_i + b_1) \dots b_{L-1}) + b_L - y_i)^2 \\ \text{such that} & \text{supp}(\mathbf{W}_\ell) \subseteq I_\ell, \forall \ell = 1, \dots, L. \end{array} \qquad \text{(Sparse-training)} $$

For which support constraints $\mathcal{I} := (I_\ell)_{\ell=1}^L$ is the optimization problem (Sparse-training) guaranteed to admit an optimum, regardless of the training data $\mathcal{D} := \{(x_i, y_i)\}_{i=1}^N$? Partial answers to this question [23] take the form of a necessary condition and a sufficient one:

* **Necessary condition:** If all instances of (Sparse-training) with support constraints $\mathcal{I}$ admit an optimum, then $\mathcal{I}$ is also a good support in the sense of deep matrix factorization, i.e., the matrix set $B_{\mathcal{I}}$ (cf. Section 3.3.1) is closed. This condition is however not sufficient: even without support constraint (in which case $B_{\mathcal{I}}$ is closed, see Box 3 for the two layer case) there are examples of simple networks where (Sparse-training) fails to have an optimum for some datasets [25].

* **Sufficient condition:** For any scalar-valued shallow network, i.e., $L = 2$ with a single output neuron, (Sparse-training) always admits an optimum, irrespective of the support constraint $\mathcal{I}$ or the dataset. This contrasts with the examples of Box 3 where $L = 2$ but the output dimension is at least 2.

It remains a challenge to find necessary and sufficient conditions (ideally associated to an efficient decision algorithm) to determine whether (Sparse-training) is guaranteed to have an optimum for a given support constraint $\mathcal{I}$. Yet, the above conditions can already be used to unveil whether there are many support constraints leading to problematic situations.

### [Box 6] Bad supports are not rare in shallow sparse networks

Consider support constraints for a shallow ReLU network: $L = 2$, and $\mathcal{I} = (I_1, I_2)$. It is known [23] that the presence of certain easy-to-detect LU patterns (cf. Box 2) in $\mathcal{I}$ prevents $B_{\mathcal{I}}$ from being closed, implying the existence of datasets for which (Sparse-training) with constraint $\mathcal{I}$ has no optimum.

How frequent is this phenomenon? We reproduce here an experiment [23] illustrating that it is all the more frequent as the support constraints are sparse. Consider a shallow ReLU network architecture where the dimension of the $L = 2$ weight matrices is $100 \times 100$. For sparsity levels $s_i \in \{0.1, 0.2, \dots, 1\}$, $i = 1, 2$, the support constraints $I_1, I_2$ are drawn uniformly at random among all constraints with cardinality $\#I_i = s_i \times (100 \times 100)$ ($s_i = 1$ thus corresponds to the absence of support constraint).
<figure><img src="image_8.png" /></figure>

For each pair $s_1, s_2$ we display above the empirical probability that at least one LU pattern is detected, implying that $B_I$ is not-closed and that (Sparse-training) can fail to have an optimum with the sparsity constraint $I$. Obviously, the phenomenon is not rare: in particular in the sparse regime $s_1, s_2 \le 0.2$, most of the randomly generated support constraints are bad. For less sparse regimes, it remains open whether the phenomenon persists because the displayed empirical probability just provides a lower bound (the presence of an LU pattern is only sufficient, but not necessary for the non-closedness of $B_I$; the latter is also only sufficient for (Sparse-training) to admit bad instances).

As illustrated in Box 6 on shallow ReLU networks, non-closedness is indeed a very common phenomenon when considering random sparse support patterns. This suggests that in order to avoid the problem of non-closedness, it is necessary to consider sparsity constraints with structures. This is the main subject of the upcoming section.

# 4 Focus on butterfly structured sparsity

In this section, we discuss in depth a structure known as “butterfly” in the literature. As we will see, the butterfly structure avoids most of the challenges of deep sparsity discussed previously.

## 4.1 Origins and expected benefits.

The butterfly structure is inspired by the computational structure of the Fast Fourier Transform (FFT) and more generally of tools appearing in Fast Multipole Methods (FMM) [8] and Hierarchical Matrices (H-matrices) [18]. Such matrices are widely used in signal processing and numerical linear algebra. More recently, the butterfly structure has acquired a quickly increasing popularity in deep learning [6,7,26] for neural network compression. This structure is indeed quite widespread, due to its powerful expressivity: a large number of classical linear operators $\mathbf{A}$ can be exactly expressed or approximated with good accuracy by a product $\mathbf{A} \approx \mathbf{W}_1 \dots \mathbf{W}_L$ of factors $\mathbf{W}_\ell$ having a prescribed structured sparse support (illustrated on Figure 4), typically with $L \approx \log_2(n)$ factors if $\mathbf{A} \in \mathbb{R}^{n \times n}$.

This structure appears for instance in the Discrete Fourier Transform (DFT) –leading to the FFT–, which in turn allows to factorize circulant matrices. Given the extreme sparsity of the corresponding support constraints (just two coefficients per row and per column in the most classical DFT when the dimension is a power of two, see Figure 2), the complexity of the corresponding matrix-vector product is reduced from $O(n^2)$ to $O(n \log_2 n)$ for $n \times n$ matrices ($2n$ operations for each factor, for $\log_2 n$ factors).
<figure><img src="image_9.png" /><figcaption>Figure 4: Illustration of the support of a factor with pattern $\pi = (a, b, c, d)$, cf. Definition 1 below. The colored squares indicate the indices belonging to the support. We illustrate respectively the concepts of factor, block and sub-block.</figcaption></figure>

## 4.2 Definition and illustration.

We discuss and illustrate a flexible parameterization of butterfly structures [24] that covers many existing variants (e.g., [6, 7, 26]) and allows for an easy mathematical and computational analysis and manipulation. In particular, the supports of sparse structured factors is encoded by four integers as follows:

**Definition 1** (Kronecker-sparse factors and their sparsity patterns [24]). For $a, b, c, d \in \mathbb{N}$, a Kronecker-sparse factor of pattern $\pi := (a, b, c, d)$ (or $\pi$-factor) is a matrix in $\mathbb{R}^{m \times n}$ or $\mathbb{C}^{m \times n}$, where $m := \text{abd}$, $n := \text{acd}$, with support included in that of the binary matrix $S_{\pi} := \mathbf{I}_a \otimes \mathbf{1}_{b \times c} \otimes \mathbf{I}_d \in \{0, 1\}^{m \times n}$ ($\otimes$ denotes the Kronecker product, see Box 7).

The tuple $\pi$ is called an elementary Kronecker-sparse pattern, or simply a pattern.

### [Box 7] The Kronecker product and its properties

Given two matrices $A \in \mathbb{R}^{n_A \times m_A}$, $B \in \mathbb{R}^{n_B \times m_B}$ where $A = (A_{i,j})_{1 \le i \le n_A, 1 \le j \le m_A}$, the Kronecker product between $A$ and $B$, denoted as $A \otimes B$, is defined as:

$$A \otimes B := \begin{pmatrix} A_{1,1}B & \dots & A_{1,m_A}B \\ \vdots & \ddots & \vdots \\ A_{n_A,1}B & \dots & A_{n_A,m_A}B \end{pmatrix} \in \mathbb{R}^{n_A n_B \times m_A m_B},$$

where $A_{i,j}B \in \mathbb{R}^{n_B \times m_B}$ is the matrix $B$ multiplying with the scalar $A_{i,j}$. An important property of the Kronecker product is the following identity:

$$(A \otimes B)(C \otimes D) = (AC) \otimes (BD), \qquad (11)$$

provided that the sizes of $A, B, C, D$ make the matrix products $AC$ and $BD$ well-defined.

Figure 4 illustrates the structure of a $\pi$-factor with $\pi = (a, b, c, d)$. Such a factor is block-diagonal, with $a$ diagonal blocks; each of these blocks it itself a block matrix of size $b \times c$, where each sub-block is a diagonal matrix of size $d \times d$. This is a flexible sparsity structure covering many examples:

* Dense matrices: a matrix of size $m \times n$ is a $(1, m, n, 1)$-factor.

* Block-diagonal matrices with $m$ sub-blocks of sizes $t \times t$ are a $(m, t, t, 1)$-factors.
* **Factors in a square dyadic butterfly factorization** [7, 24]: the pattern of the $\ell$-th factor is $\pi_{\ell} = (2^{\ell-1}, 2, 2, 2^{L-\ell})$ for $\ell = 1, \dots, L$, see Figure 2.

* **Factors in a Monarch factorization** [6]: this factorization of $m \times n$ matrices involves two factors, which patterns are $\pi_1 = (1, p, q, m/p)$, $\pi_2 = (q, m/p, n/q, 1)$ for $p, q$ such that $p \mid m$ and $q \mid n$.

A butterfly architecture $\beta$ is a sequence of patterns $\pi_{\ell} = (a_{\ell}, b_{\ell}, c_{\ell}, d_{\ell})$, $\ell = 1, \dots, L$ such that matrix multiplication between any successive $\pi_{\ell}$ and $\pi_{\ell+1}$-factors is well-defined. As an example, the factors in **square dyadic butterfly factorization** [7, 24] and **Monarch factorization** [6] admit a butterfly architecture. As described next, a key feature of butterfly architectures (and in particular of so-called *chainable* ones [24]) is that they overcome all difficulties of deep sparsity evoked so far.

## 4.3 Butterfly approximation problem

Given an architecture $\beta := (\pi_{\ell})_{\ell=1}^L$, a natural objective is to find the best approximation of any given target matrix $\mathbf{A}$ by factors admitting this architecture:

$$ \underset{\mathbf{W}_1, \dots, \mathbf{W}_L}{\text{Minimize}} \quad \|\mathbf{A} - \mathbf{W}_1 \dots \mathbf{W}_L\|_F \quad \text{s.t.} \quad \mathbf{W}_\ell \text{ is a } \pi_\ell\text{-factor, for every } \ell. \qquad (\text{BF}) $$

Is the minimum achieved? Can a (near) minimizer be computed efficiently? In the shallow case (with $L = 2$ factors), (BF) is a particular instance of (FSMF). For deeper cases ($L > 2$), this corresponds to a multi-factor version of (FSMF). Luckily, due to the special structure of Kronecker-sparse patterns, and *in contrast to the general hardness and instability results* for (FSMF) and its multi-factor version, the optimization problem (BF) is nicely behaved.

For $L = 2$, the Kronecker-sparse structure is aligned with the properties observed on the small example of Box 5. This gives rise to a polynomial algorithm to solve (BF) via truncated SVDs on blocks extracted from the target matrix $\mathbf{A}$. For the deep case, it is still an open question whether the optimal factors can be found efficiently. Nevertheless, for butterfly architectures $\beta$ satisfying a so-called *chainability* property (ensuring that the products of successive Kronecker-sparse factors in the architecture remains a Kronecker-sparse factor, see Box 8), it has been shown that a hierarchical application of the algorithm for $L = 2$ gives rise to a polynomial algorithm finding *quasi-optimal* factors $(\widehat{\mathbf{W}}_{\ell})_{\ell=1}^L$ for (BF), in the sense that:

$$ \|\mathbf{A} - \widehat{\mathbf{W}}_1 \dots \widehat{\mathbf{W}}_L\|_F \le \sqrt{(L-1)} \times \inf_{(\mathbf{W}_\ell)_{\ell=1}^L} \|\mathbf{A} - \mathbf{W}_1 \dots \mathbf{W}_L\|_F. \qquad (\text{Quasi-Opt}) $$

This applies to practical butterfly structures such as **square dyadic butterfly factorization** (the architecture behind the FFT) and **Monarch factorization** (proposed to speedup deep learning), which are indeed *chainable*. For any *chainable* architecture we may outline consequences of (Quasi-Opt):

1. When $L = 2$, the quasi-optimal factors are in fact the optimal ones since $\sqrt{L-1} = 1$. The algorithm for $L = 2$ is, in fact, the building block for the deep (i.e., $L > 2$) case.

2. When the target matrix $\mathbf{A}$ is indeed equal to the product of some (unknown) factors admitting the target (chainable) butterfly architecture $\beta$, the factors $(\widehat{\mathbf{W}}_{\ell})_{\ell=1}^L$ obtained via the algorithm indeed satisfy $\widehat{\mathbf{W}}_1 \dots \widehat{\mathbf{W}}_L = \mathbf{A}$. This is also known as the *exact recovery guarantee* [24].

3. From (Quasi-Opt) and the exact recovery guarantee, the matrix set $\mathcal{B}_I$ is closed where $I = (I_1, \dots, I_L)$, with $I_{\ell} = \text{supp}(S_{\pi_{\ell}})$. Hence (BF) admits optimal solutions for any target matrix $\mathbf{A}$.
### [Box 8] Stability of Kronecker-sparse factors under matrix multiplication – Chainability

In general, given two Kronecker-sparse patterns $\pi_1, \pi_2$, the product between $\pi_\ell$-factors is no longer a non-trivial$^a$ Kronecker-sparse factor. There are however conditions on $(\pi_1, \pi_2)$ that ensure that such a product remains a $\pi$-factor with non-trivial $\pi$. Here is example:

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

Mathematically, this can be shown by combining two following arguments:

1. If $\text{supp}(X) \subseteq \text{supp}(S_{\pi_1}), \text{supp}(Y) \subseteq \text{supp}(S_{\pi_2})$, then $\text{supp}(XY) \subseteq \text{supp}(S_{\pi_1}S_{\pi_2})$.

2. We have: $S_{\pi_1}S_{\pi_2} = S_{\pi}$. Indeed,

$$ \begin{aligned} S_{\pi_1}S_{\pi_2} &= [(\mathbf{I}_2 \otimes \mathbf{1}_{2 \times 2}) \otimes \mathbf{I}_2] [\mathbf{I}_4 \otimes \mathbf{1}_{2 \times 2}] \stackrel{(11)}{=} [(\mathbf{I}_2 \otimes \mathbf{1}_{2 \times 2})\mathbf{I}_4] \otimes [\mathbf{I}_2\mathbf{1}_{2 \times 2}] \\ &= \mathbf{I}_2 \otimes \mathbf{1}_{2 \times 2} \otimes \mathbf{1}_{2 \times 2} = \mathbf{I}_2 \otimes \mathbf{1}_{4 \times 4} = S_{\pi}. \end{aligned} $$

This example shows that for certain $(\pi_1, \pi_2)$, the products of the corresponding factors remain Kronecker-sparse, with a finer $\pi$ than dense matrices. More generally, **chainability** [24] ensures the stability of Kronecker-sparse factors under matrix multiplication.

**Definition 2 (Chainable pairs and architectures [24]).** A pair $(\pi_1, \pi_2)$ of Kronecker-sparse patterns, with $\pi_\ell = (a_\ell, b_\ell, c_\ell, d_\ell)$, $\ell = 1, 2$, is chainable if $\frac{a_1c_1}{a_2} = \frac{b_2d_2}{d_1} \in \mathbb{N}$, $a_1 | a_2$ and $d_2 | d_1$ ($a | b$ means a divides $b$). This implies the existence of $r \in \mathbb{N}$ and a Kronecker-sparse pattern $\pi$ such that $S_{\pi_1}S_{\pi_2} = rS_{\pi}$. An architecture $\beta = (\pi_\ell)_{\ell=1}^L$ is chainable if $(\pi_\ell, \pi_{\ell+1})$ is chainable for every $\ell$.

The time required by the algorithm of [24] to compute its provably good solution can be two orders of magnitude smaller for the same accuracy than gradient-based methods, with or without momentum, as illustrated in Figure 5.

## 4.4 Quantization-friendliness.

Butterfly factorization performs model compression by reducing the number of parameters to describe (or approximate) a matrix. A complementary tool is to rely on reduced-precision numerical formats via quantization. It turns out that the peculiar properties of butterfly architectures make them quantization-friendly, giving rise to specific quantization approaches that outperform the naive independent rounding of each entry of each factor. To illustrate this consider the problem of quantizing a given (square dyadic) butterfly factorization. Given $L := \log_2(n)$ Kronecker-sparse factors $\mathbf{W}_1, \dots, \mathbf{W}_L \in \mathbb{R}^{n \times n}$ (that may have been obtained using an existing algorithm to approximate some target dense matrix $\mathbf{A}$) we want to replace the

<footnote>$^a$in fact, the product can be interpreted as a dense $m \times n$ matrix and thus a $(1, m, n, 1)$-Kronecker-sparse factor</footnote>
<figure><img src="image_11.png" /><figcaption>Figure 5: Relative approximation errors $||\mathbf{A} - \widehat{\mathbf{W}}_1 \dots \widehat{\mathbf{W}}_L||_F/||\mathbf{A}||_F$ vs. running time of the different algorithms. The butterfly algorithm in the figure is from [24], and it is compared to a gradient-based method [7] and an alternating least square method [26] (see [24] for more details about this experimental comparison). The target matrix $\mathbf{A}$ is the Hadamard matrix of size $1024 \times 1024$, and $\widehat{\mathbf{W}}_1 \dots \widehat{\mathbf{W}}_L$ is the computed approximation for Problem (BF) associated with the square dyadic butterfly architecture.</figcaption></figure>

factors by quantized versions $\widehat{\mathbf{W}}_1, \dots, \widehat{\mathbf{W}}_L$, while minimizing the relative error

$$ \frac{\|\mathbf{W}_1 \dots \mathbf{W}_L - \widehat{\mathbf{W}}_1 \dots \widehat{\mathbf{W}}_L\|_F}{\|\mathbf{W}_1 \dots \mathbf{W}_L\|_F}. \qquad (12) $$

The entries of the factors $\widehat{\mathbf{W}}_\ell$ should belong to $\mathbb{F}_t$, a set of floating-point numbers with t-bit significand:

$$ \mathbb{F}_t := \{\pm k 2^{e-t}, k \in [2^{2^{t-1}}, 2^t - 1], e \in \mathbb{Z}\} \cup \{0\} \qquad (13) $$

and retain the same support as the unquantized factors: $\text{supp}(\widehat{\mathbf{W}}_\ell) \subseteq \text{supp}(\mathbf{W}_\ell)$, $\ell = 1: L$.

To achieve better quantization than naive independent rounding of each (nonzero) entry of each factor, it has been proposed [17] to exploit two key properties:

* thanks to chainability (see Box 8), the partial products $\mathbf{X} := \mathbf{W}_p \dots \mathbf{W}_q \in \mathbb{R}^{n \times n}$ and $\mathbf{Y} := \mathbf{W}_{q+1} \dots \mathbf{W}_r \in \mathbb{R}^{n \times n}$ remain Kronecker-sparse factors; one can thus exploit property (10): the product $\mathbf{X}\mathbf{Y}$ is a sum of $n$ rank-one matrices $\mathbf{C}_i := x_i y_i^\top \in \mathbb{R}^{n \times n}$ with disjoint supports.

* a recently proposed algorithm (Box 9) allows to optimally quantize rank-one matrices on $\mathbb{F}_t$.

Noticeably, this optimum quantization is generally achieved with $\lambda^* \neq 1$, meaning that it differs from the naive RTN (see definition in Box 9) strategy. Although the analysis leading to this algorithm leverages, once more, the underlying rescaling-invariance (see also Box 1), in general we also have $\mu^* \neq 1/\lambda^*$.

### [Box 9] Optimal quantization of rank-one matrices

Given high-precision vectors $x \in \mathbb{R}^n$, $y \in \mathbb{R}^m$, the rank-one matrix $xy^\top \in \mathbb{R}^{m \times n}$ can be optimally quantized in $\mathbb{F}_t$ [17] by finding the optimal scaling factor $\lambda^*$ such that, if we define $\hat{x}^* \in \mathbb{F}_t^n, \hat{y}^* \in \mathbb{F}_t^m$ as

$$ \hat{x}^* = \text{round}(\lambda^* x), \quad \hat{y}^* \in \text{round}(\mu^* y) \quad \text{with } \mu^* = \begin{cases} \frac{x^\top \hat{x}}{\|\hat{x}\|^2}, & \text{if } \hat{x} \neq 0 \\ 0, & \text{otherwise,} \end{cases} \qquad (14) $$

with $\text{round}(\cdot)$ the round-to-nearest (RTN) operator, then $\hat{x}^*$ and $\hat{y}^*$ are such that

$$ (\hat{x}^*, \hat{y}^*) \in \arg \min_{\hat{x} \in \mathbb{F}_t^m, \hat{y} \in \mathbb{F}_t^n} \|xy^\top - \hat{x}\hat{y}^\top\|^2. \qquad (15) $$
To find an optimal $\lambda^*$, an algorithm of complexity $O(mn2^t)$ has been proposed, making it tractable for the small values of $t$ that are of interest in the context of quantization to low-precision formats.

Let us see why these two properties make it almost straightforward to achieve optimal quantization in the two-factor case. We will then see how they can also be leveraged in a deeper setting.

### 4.4.1 Two-factor case

When $L = 2$ considering $X = W_1$ and $Y = W_2$ optimal quantization reads

$$ \min_{\hat{X} \in \mathcal{X}, \hat{Y} \in \mathcal{Y}} \|XY - \hat{X}\hat{Y}\|^2 \quad (16) $$

where $\mathcal{X}$ and $\mathcal{Y}$ are the sets of matrices with the same supports as $X$, $Y$ and coefficients in $\mathbb{F}_t$ (different bits per significand for $X$ and $Y$ can also be handled). Exploiting the decomposition (10) into rank-one components with disjoint supports, it is not difficult to check that this couples into $n$ independent optimal rank-one quantization problems. Its solution can thus be computed by $n$ independent applications of the optimal algorithm of Box 9 to quantize each of these rank-one matrices. This yields diagonal matrices $\Lambda, M \in \mathbb{R}^{n \times n}$ such that $\hat{X} = \text{round}(X\Lambda)$ and $\hat{Y} = \text{round}(MY)$ ($\Lambda$ and $M$ contain the optimal $\lambda^*, \mu^*$ on the diagonal).

### 4.4.2 Multi-factor case

The optimal two-factor quantization algorithm serves as a building block in a heuristic procedure to quantize a decomposition $W_1...W_L$ with $L > 2$. The key idea is to rely on two-factor quantizations, using chainability to consider certain blocks of consecutive Kronecker-sparse factors $W_\ell$ as a single Kronecker-sparse factor.

**Pairwise quantization heuristics.** For instance, we can form groups of two as

$$ (W_1W_2)(W_3W_4)...(W_{L-1}W_L) $$

(assuming that $L$ is even), and quantize each pair of consecutive factors $(W_\ell W_{\ell+1})$ separately. If the number of factors $L$ is odd, the last factor $W_L$ can simply be quantized by RTN, $\hat{W}_L = \text{round}(W_L)$. This is referred to as “pairwise” heuristic.

**Left-to-right heuristics.** Another option is to iterate and proceed from left to right. Typically at the first step the factors are parenthesized as

$$ (W_1) \underbrace{(W_2...W_L)}_{W'_2} $$

Chainability of the overall architecture is known to imply chainability of the Kronecker-sparse patterns of $W_1$ and $W'_2$, so that the optimal two-factor quantization can be run to obtain $\hat{W}_1 = \text{round}(W_1\Lambda)$, $\hat{W}'_2 = \text{round}(MW'_2) = \text{round}(MW_2...W_L)$. Keeping $\hat{W}_1$ as the quantized version of $W_1$, the process is repeated: at the generic step $\ell$ the first $\ell-1$ factors would already been quantized, and the remaining factors are regrouped as $(MW_\ell)(W_{\ell+1}...W_L)$, where the current diagonal scaling $M$ comes from the quantization at step $\ell-1$ ($M$ is the matrix of the values $\mu^*$ coming from the optimal rank-one quantization). The quantized factor $\hat{W}_\ell$ is obtained using the two-factor quantization with $X = MW_\ell$ and $Y = W_{\ell+1}...W_L$. This heuristic is referred to as left-to-right. In an analogous way a right-to-left heuristic can be defined.

**Empirical numerical efficiency.** Even if not optimal for $L > 2$, such quantization strategies show very good performance in practice, allowing for important memory savings compared
<figure><img src="image_12.png" /><figcaption>Figure 6: Relative quantization error (12) for butterfly factors with either the RTN baseline or pairwise or left-to-right algorithms, as a function of *t* and *n*. Left: *n* = 2<sup>16</sup>, *t* varies; right: *t* = 8, *n* varies.</figcaption></figure>

to the RTN strategy: Figure 6 highlights that they achieve an accuracy equivalent to the RTN
baseline using a lower precision of *t'* = *t*/1.4 for the left-to-right heuristic, which represents a
1 - 1/1.4 ≈ 30% reduction of storage for the same accuracy. Similarly, the pairwise heuristic
achieves a 1 - 1/1.3 ≈ 23% reduction.

The pairwise and the left-to-right procedures have time complexities in $O(2^t n \log n)$ and in $O(2^t n^2)$, respectively, and a space complexity in $O(2^t + n \log n)$ [17]. Notably, in practice both heuristics have much better time complexities than predicted by the theoretical bounds and both are tractable for values of $t$ of interest and for reasonably large values of $n$.

4.5 GPU-friendliness

When it comes to implementing sparse matrix multiplication on the GPU, the total running
time consists of the time to perform the computations *and the time to transfer data* between
different levels of memory on the GPU. In practice, it is observed that the time for data transfer
in the current GPU implementation of general sparse matrix multiplication that makes no
assumptions about the structure of the sparsity pattern (e.g., implementation based on the CSR
format) can become non-negligible [20], sometimes making it even slower than dense general
matrix multiplication implementations that do not assume any sparsity in the matrix (see, e.g.,
experiments in [15]). In this context, *structured sparsity*, such as N:M sparsity (where there are
at most *N* nonzero entries in every *M* consecutive entries in the same row of a matrix) or block
sparsity (where the matrix is assumed to be a block matrix, with few nonzero blocks), becomes
crucial for accelerating sparse matrix multiplication on the GPU. These structures typically
enable efficient reductions of data transfer costs in GPU implementations, thanks to the regular
memory layout of the input data.

It turns out that the butterfly sparsity also presents a form of structure that can be lever-
aged for an efficient matrix multiplication implementation on the GPU. Indeed, leveraging the
Kronecker-sparse structure (cf. Definition 1), it is possible to design a so-called *tiling strategy*
(i.e., a strategy to split the computational task into smaller subtasks that can be processed
independently and efficiently in parallel) that is well-suited for the sparse multiplication on the
GPU by a Kronecker-sparse factor [15]. Concretely, this well-adapted tiling strategy enables
a GPU implementation within a *single* kernel, eliminating the excessive memory transfer op-
erations that were a limitation of previous implementations based on batched general matrix
multiplication routines [6] (e.g., the benchmark in [15] revealed that earlier implementations
could spend up to 50% of their total runtime on read and write operations). Experimentally,
the well-adapted tiling strategy leads to a median speed-up of ×1.4 compared to previous imple-
mentations, while also cutting energy consumption by 15% [15]. In particular, it was identified
empirically that pattern $\pi = (a, b, c, d)$ with large ratio $(b+c)/bc$ leads to more time- and energy-efficient multiplication by a Kronecker-sparsity compared to its dense counterpart [15]. Such a heuristic lays the groundwork to design more efficient Kronecker-sparse neural networks in the future, for instance by selecting the most efficient pattern $(a, b, c, d)$ among those with the same number of non-zeros.

# 5 Perspectives and Open Challenges

The main lesson of this survey is that, while deep sparsity remains an attractive and natural objective to enable much needed flexible resource-performance tradeoffs in learning, promoting it in a deep context is far from straightforward. This is indeed not just a trivial transposition from the shallow context as, as we have seen, current sparse know-how severely breaks down in a deep context and unforeseen difficulties arise: ill-posed training problems, NP-hard, possibly with no optimum.

A first challenge is to fully harness the natural rescaling-invariances of deep networks, from quantization and pruning to optimization. The promise is to achieve principled and robust approaches to train and compress deep ReLU networks, fully leveraging and extending the successful schemes developed so far in the deep linear case for butterfly matrix factorizations.

It also becomes of crucial importance to clearly distinguish between “good” and “bad” problems, and to have a mean to individuate and construct the good ones, for instance through a tractable and efficient algorithm to check if a selected support is pathological or not. A related problem, in the context of general DAG ReLU network architectures, is to understand which DAG structures correspond to good optimization problems, possibly via an algebraic characterization involving the corresponding path-lifting. This would significantly consolidate algorithmic foundations of pruning, by ensuring that the sparsity they induce maintains stability of the subsequent re-training steps. From a more mathematical perspective, knowing that an optimization problem admits a global optimum is also an important step to enable the analysis of the properties of this optimum, as was done e.g., with success to characterize the success of $ℓ^1$ minimization in linear inverse problems.

Kronecker sparsity (extending classical butterfly sparsity) is a promising first step to build a comprehensive framework for *structured deep sparsity*. Yet, while butterflies ensure stable and tractable learning, such guarantees so far only hold in a matrix factorization setting and extending such architectures and the associated guarantees to nonlinear contexts, e.g., *ReLU networks*, is a promising perspective.

A particular challenge is to handle the natural permutation-invariance of standard network layers (at the neuron level of the channel level for convolutive layers). Indeed it is documented that even if a matrix admits an exact butterfly factorization, permuting its rows and/or columns can severely degrade its error of best approximation via a butterfly matrix [41]. Moreover, mixing butterfly structures for “dense” network layers with other types of structured sparse approximations of the tensors appearing in other deep learning modules such as convolutional layers is desirable, and the similarity between butterfly approximation algorithms [24] and near-optimal algorithms for tree-based tensor network approximation [9] offer promising avenues. Special architectures, such as neural operators, seem particularly suited for such approximations since butterfly structures are more likely to appear naturally in this context.

Overall, the main challenge is to design deep sparsity structures that combine the provably good computational properties of butterflies with the ability to closely match the performance of today’s best networks on real learning problems, while being endowed with statistical generalization guarantees.
# 6 acknowledgements

This work was supported in part by the AllegroAssai ANR-19-CHIA-0009 project of the French Agence Nationale de la Recherche and by the SHARP ANR project ANR-23-PEIA-0008 in the context of the France 2030 program. Tung Le was supported by AI Interdisciplinary Institute ANITI funding, through the French "Investments for the Future – PIA3" program under the grant agreement ANR-19-PI3A-0004, and Air Force Office of Scientific Research, Air Force Material Command, USAF, under grant numbers FA8655-22-1-7012. The authors thank the Blaise Pascal Center of Ecole Normale Supérieure de Lyon for the computational means. It uses the SIDUS solution [33] developed by Emmanuel Quemener.

# References

[1] R. G. Baraniuk, E. Candes, R. Nowak, and M. Vetterli. Compressive Sampling [From the Guest Editors]. *IEEE Signal Processing Magazine*, 25(2):12–13, Mar. 2008. Conference Name: IEEE Signal Processing Magazine.

[2] A. Benichoux, E. Vincent, and R. Gribonval. A fundamental pitfall in blind deconvolution with sparse and shift-invariant priors. In *ICASSP*, 2013.

[3] J. Bona-Pellissier, F. Malgouyres, and F. Bachoc. Local Identifiability of Deep ReLU Neural Networks: the Theory. In *NeurIPS*, 2022.

[4] J. E. Cohen and V. Leplat. Efficient Algorithms for Regularized Nonnegative Scale-invariant Low-rank Approximation Models, June 2024. arXiv:2403.18517.

[5] M. J. Colbrook, V. Antun, and A. C. Hansen. The difficulty of computing stable and accurate neural networks: On the barriers of deep learning and smale's 18th problem. *Proceedings of the National Academy of Sciences*, 119(12):e2107151119, 2022.

[6] T. Dao, B. Chen, N. Sohoni, A. Desai, M. Poli, J. Grogan, A. Liu, A. Rao, A. Rudra, and C. Ré. Monarch: Expressive structured matrices for efficient and accurate training. In *ICML*, 2022.

[7] T. Dao, A. Gu, M. Eichhorn, A. Rudra, and C. Ré. Learning fast algorithms for linear transforms using butterfly factorizations. In *ICML*, pages 1517–1527, 2019.

[8] N. Engheta, W. D. Murphy, V. Rokhlin, and M. S. Vassiliou. The fast multipole method (FMM) for electromagnetic scattering problems. *IEEE Transactions on Antennas and Propagation*, 40(6):634–641, 1992.

[9] A. Falco, W. Hackbusch, and A. Nouy. Tree-based tensor formats. *SeMA Journal*, 78(2):159–173, June 2021. arXiv:1810.01262.

[10] S. Foucart and H. Rauhut. *A Mathematical Introduction to Compressive Sensing*. Springer, New York, NY, 2013.

[11] J. Frankle, D. J. Schwab, and A. S. Morcos. The early phase of neural network training. In *ICLR*, 2020.

[12] N. Gillis and F. Glineur. Low-Rank Matrix Approximation with Weights or Missing Data Is NP-Hard. *SIAM Journal on Matrix Analysis and Applications*, 32(4):1149–1165, Oct. 2011.

[13] A. Gonon, N. Brisebarre, E. Riccietti, and R. Gribonval. A path-norm toolkit for modern networks: consequences, promises and challenges. In *ICLR*, 2024.
[14] A. Gonon, N. Brisebarre, E. Riccietti, and R. Gribonval. A rescaling-invariant Lipschitz bound based on path-metrics for modern ReLU network parameterizations, Jan. 2025. arXiv:2405.15006.

[15] A. Gonon, L. Zheng, P. Carrivain, and Q.-T. Le. Fast Inference with Kronecker-Sparse Matrices. In *ICML 2025*, Vancouver (BC), Canada, July 2025.

[16] Y. Grandvalet. Least Absolute Shrinkage is Equivalent to Quadratic Penalization. In *ICANN 98*, Perspectives in Neural Computing, pages 201–206, London, 1998. Springer.

[17] R. Gribonval, T. Mary, and E. Riccietti. Optimal quantization of rank-one matrices in floating-point arithmetic—with applications to butterfly factorizations. working paper or preprint, June 2023.

[18] W. Hackbusch. *Hierarchical Matrices: Algorithms and Analysis*, volume 49 of Springer Series in Computational Mathematics. Springer, Berlin, Heidelberg, 2015.

[19] S. Han, J. Pool, J. Tran, and W. J. Dally. Learning both weights and connections for efficient neural networks. In *NeurIPS*, 2015.

[20] T. Hoefler, D. Alistarh, T. Ben-Nun, N. Dryden, and A. Peste. Sparsity in deep learning: Pruning and growth for efficient inference and training in neural networks. *Journal of Machine Learning Research*, 22(241):1–124, 2021.

[21] F. Latorre, A. Bonnet, P. Rolland, N. Hallak, and V. Cevher. 1-path-norm regularization of deep neural networks. In *LatinX in AI Workshop at ICML 2023*, 2023.

[22] Q.-T. Le, E. Riccietti, and R. Gribonval. Spurious Valleys, NP-hardness, and Tractability of Sparse Matrix Factorization With Fixed Support. *SIAM Journal on Matrix Analysis and Applications*, 2022.

[23] Q.-T. Le, E. Riccietti, and R. Gribonval. Does a sparse ReLU network training problem always admit an optimum? In *NeurIPS*, 2023.

[24] Q.-T. Le, L. Zheng, E. Riccietti, and R. Gribonval. Butterfly factorization with error guarantees, Nov. 2024. arXiv:2411.04506.

[25] L.-H. Lim, M. Michalek, and Y. Qi. Best $k$-layer neural network approximations. *Con-structive Approximation*, June 2021.

[26] R. Lin, J. Ran, K. H. Chiu, G. Chesi, and N. Wong. Deformable butterfly: A highly structured and sparse linear transform. In *NeurIPS*, 2021.

[27] S. G. Mallat and Z. Zhang. Matching pursuits with time-frequency dictionaries. *IEEE Trans. Sig. Proc.*, 41(12):3397–3415, 1993.

[28] B. Natarajan. Sparse approximate solutions to linear systems. *SIAM J. Computing*, 25(2):227–234, 1995.

[29] B. Neyshabur, R. Salakhutdinov, and N. Srebro. Path-SGD - Path-Normalized Optimization in Deep Neural Networks. In *NeurIPS*, 2015.

[30] R. Parhi and R. D. Nowak. Deep Learning Meets Sparse Regularization: A signal processing perspective. *IEEE Signal Processing Magazine*, 40(6):63–74, Sept. 2023.

[31] M. Pilanci and T. Ergen. Neural Networks are Convex Regularizers: Exact Polynomial-time Convex Optimization Formulations for Two-layer Networks. In *ICML 2020*.
[32] C. Poon and G. Peyré. Smooth over-parameterized solvers for non-smooth structured optimization. *Mathematical Programming*, 201(1):897–952, Sept. 2023.

[33] E. Quemener and M. Corvellec. Sidus—the solution for extreme deduplication of an operating system. *Linux Journal*, 2013.

[34] J. Romberg. Compressive Sensing by Random Convolution. *SIAM Journal on Imaging Sciences*, 2(4):1098–1128, Jan. 2009.

[35] V. Silva and L.-H. Lim. Tensor rank and the ill-posedness of the best low-rank approximation problem. *SIAM Journal on Matrix Analysis and Applications*, 30:1084–1127, 08 2006.

[36] P. Stock, B. Graham, R. Gribonval, and H. Jégou. Equi-normalization of Neural Networks. In *ICLR*, May 2019.

[37] P. Stock and R. Gribonval. An Embedding of ReLU Networks and an Analysis of their Identifiability. *Constructive Approximation*, 2023.

[38] R. J. Tibshirani. Equivalences Between Sparse Models and Neural Networks, 2021. Online working notes.

[39] I. Tosic and P. Frossard. Dictionary Learning. *IEEE Signal Processing Magazine*, 28(2):27–38, Mar. 2011.

[40] S. A. Vavasis. On the Complexity of Nonnegative Matrix Factorization. *SIAM Journal on Optimization*, 20(3):1364–1377, Jan. 2010.

[41] L. Zheng, G. Puy, E. Riccietti, P. Pérez, and R. Gribonval. Butterfly factorization by algorithmic identification of rank-one blocks. In *XXIXème Colloque Francophone de Traitement du Signal et des Images*, Grenoble, France, Aug. 2023.