# Robust Deconvolution with Parseval Filterbanks

Rossen Nenov, Vincent Lostanlen, Peter Balazs

► To cite this version:

Rossen Nenov, Vincent Lostanlen, Peter Balazs. Robust Deconvolution with Parseval Filterbanks. IEEE International Conference on Sampling Theory and Applications (SampTA), Jul 2025, Vienne, Austria. hal-05006852

HAL Id: hal-05006852

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

Submitted on 26 Mar 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Robust Deconvolution with Parseval Filterbanks

Rossen Nenov*, Vincent Lostanlen†, and Peter Balazs*‡

*Acoustics Research Institute, Austrian Academy of Sciences, Vienna, Austria

†Nantes Université, École Centrale Nantes, LS2N, CNRS, UMR 6004, F-44000 Nantes, France

‡Interdisciplinary Transformation University Austria (IT:U), Linz, Austria

rossen.nenov@oeaw.ac.at

**Abstract**—This article introduces two contributions: Multi-band Robust Deconvolution (Multi-RDCP), a regularization approach for deconvolution in the presence of noise; and Subband-Normalized Adaptive Kernel Evaluation (SNAKE), a first-order iterative algorithm designed to efficiently solve the resulting optimization problem. Multi-RDCP resembles Group LASSO in that it promotes sparsity across the subband spectrum of the solution. We prove that SNAKE enjoys fast convergence rates and numerical simulations illustrate the efficiency of SNAKE for deconvolving noisy oscillatory signals.

**Index Terms**—Deconvolution, digital signal processing, gradient methods, signal reconstruction, sparse approximation.

## I. INTRODUCTION

Deconvolution is a well-known inverse problem with numerous applications in signal processing [1], [2]. The (circular) convolution between two finite sequences $x$ and $w$ in $\mathbb{R}^N$ is:

$$ (w * x)[n] = \sum_{m=0}^{N-1} w[m]x[(n-m) \mod N]. \quad (1) $$

Reciprocally, (non-blind) deconvolution operators assume the knowledge of two finite sequences $x$, $y$, and seek $w$ such that $w * x \approx y$; more precisely, such that a certain cost function $F$ is minimized. One such function is mean square error (MSE):

$$ f: w \in \mathbb{R}^N \mapsto \frac{1}{2} \|w * x - y\|^2. \quad (2) $$

If one simply chooses to minimize Equation (2), deconvolution can be solved directly via the fast Fourier transform (FFT) [3] by spectral division. Denoting the discrete Fourier transform of $x$ by $\mathcal{F}(x) = \hat{x}$, the convolution theorem [4] yields, for every frequency $k$: $\mathcal{F}(w * x)[k] = \hat{w}[k]\hat{x}[k]$. Let $w^* = \mathcal{F}^{-1}(\hat{w}^*)$ where, for every $k$, $\hat{w}^*[k] = \hat{y}[k]/\hat{x}[k]$ if $\hat{x}[k] \neq 0$ and $\hat{w}^*[k] = 0$ otherwise. One may verify that $w^*$ minimizes the MSE. Indeed, $f$ is convex and its gradient is zero at $w^*$ [4].

A different approach is needed for regularized deconvolution problems, whose cost function writes as $F = f + r$. Here, $f$ is a data fidelity term-like MSE-depending on $x$ and $y$, while $r$ penalizes a certain notion of “complexity” in the solution. The practical interest behind these problems lies in the need to reliably estimate the impulse response of the linear time-invariant (LTI) system which maps $x$ to $y$ in the presence of highly noisy observations. The complexity of the solution

Corresponding author: Rossen Nenov. The work of R. Nenov and P. Balazs is supported by the Austrian Science Fund (FWF) projects LoFT [10.55776/P34624], NoMASP [10.55776/P34922] and Voice Prints [10.55776/P36446] and the WWTF project EleCom [10.47379/LS23024]. The work of V. Lostanlen is supported by ANR project MuReNN (ANR-23-CE23-0007-01). To reproduce our numerical simulations, please visit: https://github.com/RosseN16/SNAKE

may arguably be quantified in terms of spectral spread, i.e., the notion that the entries $|\hat{w}[k]|^2$ are non-negligible for a wide range of frequency values $k$. Narrowband LTI systems are ubiquitous in applied sciences: formants in phonetics [5] and Helmholtz resonators in musical acoustics [6] are some well-known examples.

In this article, we introduce the multiband robust deconvolution problem (Multi-RDCP). The key idea is to employ a predefined filterbank $\Psi$ to decompose $x$ and $y$ into $J > 1$ subbands. This decomposition enables the formulation of individual subproblems and the regularization of the solution $w$ through a non-squared $\ell^2$-norm penalty applied to each subband component, promoting sparsity across subbands, akin to Group LASSO [7]. Since Multi-RDCP lacks closed-form solutions and full differentiability, an iterative method such as proximal gradient descent (PGD) [8] is required. Yet, we will see that general PGD is impractically slow if $\hat{x}$ exhibits large variations in magnitude across regions of the Fourier spectrum. In order to overcome this, we present a new algorithm: subband-normalized adaptive kernel evaluation (SNAKE). SNAKE applies PGD over each subband independently, with a step size that is adapted to the power spectral density of $x$ within that band.

Our main theoretical result (Theorem III.6) is that SNAKE enjoys a global linear convergence rate. We corroborate this result through numerical experiments, demonstrating significant convergence speedups compared to classical PGD, along with an improved robustness to noise.

**Note.** We write $[N]$ as a shorthand for $\{0, \dots, N-1\}$. We denote the time reversal of a finite sequence $x \in \mathbb{R}^N$ by $\tilde{x}$.

## II. PROXIMAL GRADIENT-BASED DECONVOLUTION

**Definition II.1.** Given $\lambda > 0$ and $x, y \in \mathbb{R}^N$, the robust deconvolution problem (RDCP) of $y$ from $x$ minimizes

$$ F: w \in \mathbb{R}^N \mapsto \frac{1}{2} \|w * x - y\|^2_2 + \lambda \|w\|_2, \quad (3) $$

i.e., $F = f + r$ with $f$ is the MSE and $r$ the $\ell^2$-norm times $\lambda$. **Note.** We square the $\ell^2$ norm in $f$ but not in $r$. This is an important design choice, which we will discuss in Section III.

A common choice for minimizing the sum of a smooth and a non-smooth convex function is proximal gradient descent.

**Definition II.2** ([9]). The proximal operator of a convex, lower semi-continuous function $g: \mathbb{R}^N \to \mathbb{R}$ is defined as

$$ \operatorname{Prox}_{\alpha g}: w \in \mathbb{R}^N \mapsto \arg \min_{v \in \mathbb{R}^N} \left\{ g(v) + \frac{1}{2\alpha} \|v - w\|^2_2 \right\}. \quad (4) $$
**Lemma II.3** ([10], Section 6.5.1). Given $\alpha > 0$ and $\lambda > 0$, the proximal operator of $r: \mathbf{w} \in \mathbb{R}^N \mapsto \lambda \|\mathbf{w}\|_2$ is given by

$$ \text{Prox}_{\alpha r}: \mathbf{w} \mapsto \max(0, 1 - \alpha \lambda / \|\mathbf{w}\|_2) \mathbf{w}. \quad (5) $$

**Definition II.4** ([8]). Given $\mathbf{x}, \mathbf{y} \in \mathbb{R}^N$ and RDCP (Def. II.1), an initial guess $\mathbf{w}^{(0)} \in \mathbb{R}^N$ and a step size $\alpha > 0$, proximal gradient descent (PGD) computes the sequence $\{\mathbf{w}^{(t)}\}_{t \in \mathbb{N}}$ recursively:

$$ \mathbf{w}^{(t+1)} = \text{Prox}_{\alpha r} \left( \mathbf{w}^{(t)} - \alpha \nabla f(\mathbf{w}^{(t)}) \right). \quad (6) $$

PGD theory has first been developed under the assumption of strong convexity. More recently, [11] have proved the linear convergence of PGD under a weaker assumption than strong convexity, called “proximal Polyak-Łojasiewicz” (Prox-PŁ) in reference to [12] and [13]. We proceed to establish the same result under an even weaker assumption than Prox-PŁ.

First we state some results on the MSE $f$ in RDCP.

**Definition II.5** ([14]). For a $\mathbf{x} \in \mathbb{R}^N$ the associated circulant matrix $\mathbf{C}_\mathbf{x} \in \mathbb{R}^{N \times N}$, such that $\mathbf{w} * \mathbf{x} = \mathbf{C}_\mathbf{x} \mathbf{w}$, is given by $\mathbf{C}_\mathbf{x}[n, n'] = \mathbf{x}[(n - n') \text{ mod } N]$ for $n, n' \in [N]$.

**Definition II.6.** Given $\mathbf{x} \in \mathbb{R}^N$, we denote its smallest nonzero and largest power spectrum entries by

$$ S_{\min}^{+}(\mathbf{x}) = \min_{k \in [N]} \left\{ |\hat{\mathbf{x}}|^2[k] \middle| |\hat{\mathbf{x}}|^2[k] > 0 \right\}, \quad (7) $$

$$ S_{\max}(\mathbf{x}) = \max_{k \in [N]} \left\{ |\hat{\mathbf{x}}|^2[k] \right\}, \quad (8) $$

where $\hat{\mathbf{x}}$ is the discrete Fourier transform of $\mathbf{x}$.

**Lemma II.7.** Given $\mathbf{w}, \mathbf{v} \in \mathbb{R}^N$, Equation (2) satisfies: $f(\mathbf{v}) = f(\mathbf{w}) + \langle \nabla f(\mathbf{w}), \mathbf{v} - \mathbf{w} \rangle + \frac{1}{2} \|\mathbf{w} * \mathbf{x} - \mathbf{w} * \mathbf{x}\|_2^2$.

**Proof.** We rewrite $f$ in terms of $\mathbf{C}_\mathbf{x}$ and expand the norm to obtain $f(\mathbf{v}) = \frac{1}{2} (\|\mathbf{C}_\mathbf{x} \mathbf{w} - \mathbf{y}\|_2^2 + 2 \langle \mathbf{C}_\mathbf{x} \mathbf{w} - \mathbf{y}, \mathbf{C}_\mathbf{x} \mathbf{v} - \mathbf{C}_\mathbf{x} \mathbf{w} \rangle + \|\mathbf{C}_\mathbf{x} \mathbf{v} - \mathbf{C}_\mathbf{x} \mathbf{w}\|_2^2)$. Since the gradient is given by $\nabla f(\mathbf{w}) = \mathbf{C}_\mathbf{x} (\mathbf{C}_\mathbf{x} \mathbf{w} - \mathbf{y})$, we conclude with shifting the transpose in $\langle \nabla f(\mathbf{w}), \mathbf{w} - \mathbf{v} \rangle = \langle \mathbf{C}_\mathbf{x} \mathbf{w} - \mathbf{y}, \mathbf{C}_\mathbf{x} \mathbf{w} - \mathbf{C}_\mathbf{x} \mathbf{v} \rangle$. $\square$

**Lemma II.8.** For $\mathbf{v} \in \mathbb{R}^N$, $\|\mathbf{v} * \mathbf{x}\|_2^2 \le S_{\max}(\mathbf{x}) \|\mathbf{v}\|_2^2$. For all $\mathbf{w} \in \ker(\mathbf{C}_\mathbf{x})^{\perp}$, $S_{\min}^{+}(\mathbf{x}) \|\mathbf{w}\|_2^2 \le \|\mathbf{w} * \mathbf{x}\|_2^2$.

**Proof.** The matrix $\mathbf{C}_\mathbf{x}$ is diagonalized by the discrete Fourier transform and its singular values are the entries of $\hat{\mathbf{x}}$ [15]. $\square$

We now introduce a measure of the distance of $\mathbf{w}$ to the set of optimal solutions in terms of function values and we establish its connection to the iterates of PGD for RDCP.

**Definition II.9.** The proximal gradient optimality gap of RDCP (Def. II.1), for $\mathbf{w} \in \mathbb{R}^N$ and $\gamma > 0$, is given by:

$$ D_r(\mathbf{w}, \gamma) := \begin{aligned} & -\gamma \min_{\mathbf{v} \in \mathcal{V}} \left\{ \langle \nabla f(\mathbf{w}), \mathbf{v} - \mathbf{w} \rangle + \frac{\gamma}{2} \|\mathbf{v} - \mathbf{w}\|_2^2 + r(\mathbf{v}) - r(\mathbf{w}) \right\}, \\ & \text{where } \mathcal{V} = \text{ran}(\mathbf{C}_\mathbf{x}^\top) \text{ is the row space of the matrix } \mathbf{C}_\mathbf{x}. \end{aligned} \quad (9) $$

*Note.* This definition is similar to [11, Equation 13] except that we restrict the minimum to $\mathcal{V}$, i.e., a linear subspace of $\mathbb{R}^N$. For $\mathbf{v} \in \mathcal{V}$, we have $\|\mathbf{v} * \mathbf{x} - \mathbf{w} * \mathbf{x}\|_2^2 > 0$ unless $\mathbf{v} = \mathbf{w}$. Indeed, $\mathcal{V}$ is orthogonal to $\ker(\mathbf{C}_\mathbf{x})$, along which $f$ is constant.

**Lemma II.10.** The proximal gradient optimality gap of RDCP (Def. II.1) is monotonically nondecreasing with respect to its second argument: $\forall \mathbf{w} \in \mathbb{R}^N, \forall \gamma, \gamma' > 0$:

$$ (\gamma < \gamma') \implies (D_r(\mathbf{w}, \gamma) \le D_r(\mathbf{w}, \gamma')). \quad (10) $$

**Proof.** This can be proven analogously to [11, Appendix E, Lemma 1] by ensuring that the same arguments apply if the minimum is taken over $\mathcal{V} = \text{ran}(\mathbf{C}_\mathbf{x}^\top)$ instead of over $\mathbb{R}^N$. $\square$

**Lemma II.11.** Given $\mathbf{x}, \mathbf{y} \in \mathbb{R}^N$ and RDCP (Def. II.1), an initial guess $\mathbf{w}^{(0)} \in \text{ran}(\mathbf{C}_\mathbf{x}^\top)$ and a step size $\alpha = \frac{1}{L}$. Then the iterates of PGD satisfy $\mathbf{w}^{(t)} \in \text{ran}(\mathbf{C}_\mathbf{x}^\top)$ for $t \in \mathbb{N}$. Furthermore, the following holds:

$$ \langle \nabla f(\mathbf{w}^{(t)}), \mathbf{w}^{(t+1)} - \mathbf{w}^{(t)} \rangle + \frac{L}{2} \|\mathbf{w}^{(t+1)} - \mathbf{w}^{(t)}\|_2^2 + r(\mathbf{w}^{(t+1)}) - r(\mathbf{w}^{(t)}) = -\frac{1}{L} D_r(\mathbf{w}^{(t)}, L). $$

**Proof.** Using Lemma II.3, we verify that, for any given $\mathbf{w} \in \text{ran}(\mathbf{C}_\mathbf{x}^\top)$, $\nabla f(\mathbf{w})$ and $\text{prox}_{\alpha r}(\mathbf{w})$ are both in $\text{ran}(\mathbf{C}_\mathbf{x}^\top)$. Thus, starting from $\mathbf{w}^{(0)} \in \text{ran}(\mathbf{C}_\mathbf{x}^\top)$, we prove by recurrence that $\mathbf{w}^{(t)} \in \text{ran}(\mathbf{C}_\mathbf{x}^\top)$ for every $t \in \mathbb{N}$. Furthermore, due to the definition of the proximal operator and $\mathbf{w}^{(t+1)}$, we have

$$ \begin{aligned} \mathbf{w}^{(t+1)} &= \arg \min_{\mathbf{v} \in \mathbb{R}^N} \{r(\mathbf{v}) + \frac{L}{2} \|\mathbf{v} - \mathbf{w}^{(t)} + \frac{1}{L} \nabla f(\mathbf{w}^{(t)})\|_2^2\} \\ &= \arg \min_{\mathbf{v} \in \mathbb{R}^N} \{r(\mathbf{v}) + \langle \nabla f(\mathbf{w}^{(t)}), \mathbf{v} - \mathbf{w}^{(t)} \rangle + \frac{L}{2} \|\mathbf{v} - \mathbf{w}^{(t)}\|_2^2\} \end{aligned} $$

with the omission of $\frac{1}{L} \|\nabla f(\mathbf{w}^{(t)})\|_2^2$, a term which does not affect the argmin. From Definition II.9, we see that:

$$ -\frac{1}{L} D_r(\mathbf{w}^{(t)}, L) + r(\mathbf{w}^{(t)}) = \min_{\mathbf{v} \in \mathcal{V}} \left[ r(\mathbf{v}) + \langle \nabla f(\mathbf{w}^{(t)}), \mathbf{v} - \mathbf{w}^{(t)} \rangle + \frac{L}{2} \|\mathbf{v} - \mathbf{w}^{(t)}\|_2^2 \right]. $$

The minimum is reached for $\mathbf{v} = \mathbf{w}^{(t+1)}$, which we have proven to be in $\mathcal{V}$; thus concluding the second statement. $\square$

**Lemma II.12.** Given $\mathbf{x}, \mathbf{y} \in \mathbb{R}^N$ and RDCP (Def. II.1). Then $\arg \min_{\mathbf{v} \in \mathbb{R}^N} F(\mathbf{v}) \subset \text{ran}(\mathbf{C}_\mathbf{x}^\top)$.

**Proof.** Since $F(\mathbf{v}) \to +\infty$ if $\|\mathbf{v}\|_2 \to +\infty$, $\arg \min_{\mathbf{v} \in \mathbb{R}^N} F$ is nonempty [16]. Let $\mathbf{v}^* \in \arg \min_{\mathbf{v} \in \mathbb{R}^N} F(\mathbf{v})$. If $\mathbf{v}^* = 0$, the result is trivial. Otherwise, we write the optimality condition: $0 = \nabla f(\mathbf{v}^*) + \partial r(\mathbf{v}^*) = \mathbf{C}_\mathbf{x}^\top (\mathbf{C}_\mathbf{x} \mathbf{v}^* - \mathbf{y}) + (\lambda / \|\mathbf{v}^*\|_2) \mathbf{v}^*$. $\square$

**Proposition II.13.** Given $\mathbf{x}, \mathbf{y} \in \mathbb{R}^N$ and RDCP (Def. II.1), for every $\mathbf{w} \in \text{ran}(\mathbf{C}_\mathbf{x}^\top)$:

$$ D_r(\mathbf{w}, L) \ge \mu(F(\mathbf{w}) - F^*), \quad (11) $$

where $\mu = S_{\min}^{+}(\mathbf{x})$, $L = S_{\max}(\mathbf{x})$, and $F^* = \min_{\mathbb{R}^N} F$.

**Proof.** For $\mathbf{v}, \mathbf{w} \in \text{ran}(\mathbf{C}_\mathbf{x}^\top) = \ker(\mathbf{C}_\mathbf{x})^{\perp}$, Lemma II.8 yields $\frac{1}{2} \|\mathbf{v} * \mathbf{x} - \mathbf{w} * \mathbf{x}\|_2^2 \ge \frac{\mu}{2} \|\mathbf{v} - \mathbf{w}\|_2^2$. Plugging this into Lemma II.7 yields: $f(\mathbf{v}) \ge f(\mathbf{w}) + \langle \nabla f(\mathbf{w}), \mathbf{v} - \mathbf{w} \rangle + \frac{\mu}{2} \|\mathbf{v} - \mathbf{w}\|_2^2$. We add $r(\mathbf{v})$, take the minimum over all elements $\mathbf{v} \in \text{ran}(\mathbf{C}_\mathbf{x}^\top)$, and recognize Definition II.9:

$$ \min_{\mathbf{v} \in \text{ran}(\mathbf{C}_\mathbf{x}^\top)} F(\mathbf{v}) \ge F(\mathbf{w}) - \frac{1}{\mu} D_r(\mathbf{w}, \mu). $$

On the left, we recognize $F^*$ (Lemma II.12). On the right, we use Lemma II.10 to replace $D_r(\mathbf{w}, \mu)$ by $D_r(\mathbf{w}, L)$. $\square$
**Theorem II.14.** Given $\mathbf{x}, \mathbf{y} \in \mathbb{R}^N$ and RDCP (Def. II.1), an initial guess $\mathbf{w}^{(0)} \in \text{ran}(\mathbf{C}_\mathbf{x}^\top)$ and a step size $\alpha = 1/S_{\text{max}}(\mathbf{x})$, the iterates $\{\mathbf{w}^{(t)}\}_{t \in \mathbb{N}}$ of PGD satisfy:

$$ \frac{F(\mathbf{w}^{(t)}) - F^*}{F(\mathbf{w}^{(0)}) - F^*} \le \left(1 - \frac{S_{\text{min}}^+(\mathbf{x})}{S_{\text{max}}(\mathbf{x})}\right)^t. \quad (12) $$

*Proof.* We adapt the proof of [11, Theorem 5]. Applying Lemma II.7 for $\mathbf{v} = \mathbf{w}^{(t+1)}$ and $\mathbf{w} = \mathbf{w}^{(t)}$, and Lemma II.7 for $\frac{1}{2}\|\mathbf{C}_\mathbf{x}(\mathbf{w}^{(t+1)} - \mathbf{w}^{(t)})\|_2^2 \le \frac{L}{2}\|\mathbf{w}^{(t+1)} - \mathbf{w}^{(t)}\|_2^2$ yield:

$$ f(\mathbf{w}^{(t+1)}) \le f(\mathbf{w}^{(t)}) + \langle \nabla f(\mathbf{w}^{(t)}), \mathbf{w}^{(t+1)} - \mathbf{w}^{(t)} \rangle + \frac{L}{2} \|\mathbf{w}^{(t+1)} - \mathbf{w}^{(t)}\|_2^2 + r(\mathbf{w}^{(t+1)}) - r(\mathbf{w}^{(t)}) \\ + \frac{L}{2} \|\mathbf{w}^{(t+1)} - \mathbf{w}^{(t)}\|_2^2. $$

Combining this with Lemma II.11 and Proposition II.13 yields:

$$ \begin{aligned} F(\mathbf{w}^{(t+1)}) &= f(\mathbf{w}^{(t+1)}) + r(\mathbf{w}^{(t)}) + r(\mathbf{w}^{(t+1)}) - r(\mathbf{w}^{(t)}) \\ &\le F(\mathbf{w}^{(t)}) + \langle \nabla f(\mathbf{w}^{(t)}), \mathbf{w}^{(t+1)} - \mathbf{w}^{(t)} \rangle \\ &\quad + \frac{L}{2} \|\mathbf{w}^{(t+1)} - \mathbf{w}^{(t)}\|_2^2 + r(\mathbf{w}^{(t+1)}) - r(\mathbf{w}^{(t)}) \\ &= F(\mathbf{w}^{(t)}) - \frac{1}{L} D_r(\mathbf{w}^{(t)}, L) \\ &\le F(\mathbf{w}^{(t)}) - \frac{\mu}{L} [F(\mathbf{w}^{(t)}) - F^*]. \end{aligned} $$

This gives us $F(\mathbf{w}^{(t+1)}) - F^* \le (1 - \frac{\mu}{L})[F(\mathbf{w}^{(t)}) - F^*]$, from which we can conclude the statement by recurrence. □

Theorem II.14 guarantees that PGD applied to robust deconvolution of $\mathbf{y}$ from $\mathbf{x}$ converges in terms of optimal function value with a rate of $(1 - S_{\text{min}}^+(\mathbf{x})/S_{\text{max}}(\mathbf{x}))^t$ or better.

### III. FILTERBANK ANALYSIS AND NORMALIZATION

In many applications, the ratio between the smallest nonzero and large power spectrum coefficients of $\mathbf{x}$ is close to zero, hence a poor convergence rate guarantee in Equation (12). Against this issue, we use a filterbank to partition RDCP into subproblems, each focusing on a subband, i.e., a subset of frequency components. By assuming that the power spectral density of $\mathbf{x}$ has small relative variations within each subband, we obtain subband-wise power spectrum ratios which are closer to one, yielding improved convergence rate guarantees.

**Definition III.1.** A Parseval filterbank is a collection of finite sequences $\psi_0, \dots, \psi_{J-1} \in \mathbb{R}^N$ such that, for every $\mathbf{x} \in \mathbb{R}^N$:

$$ \sum_{j=0}^{J-1} \|\psi_j * \mathbf{x}\|_2^2 = \|\mathbf{x}\|_2^2. \quad (13) $$

**Theorem III.2 ([17]).** Given $\mathbf{w} \in \mathbb{R}^N$ and a Parseval filterbank $(\psi_0, \dots, \psi_{J-1})$: $\mathbf{w} = \sum_{j=0}^{J-1} \mathbf{C}_{\psi_j}^\top \mathbf{C}_{\psi_j} \mathbf{w}$.

**Definition III.3.** Let $\Psi = (\psi_0, \dots, \psi_{J-1})$ be a Parseval filterbank. Choose $(\lambda_j)_{j=0}^{J-1}$ and $(L_j)_{j=0}^{J-1}$ such that $\lambda_j \ge 0$ and $L_j > 0$ for every $j \in [J]$. Given $\mathbf{x}, \mathbf{y} \in \mathbb{R}^N$, for each $j \in [J]$, define the data fidelity terms

$$ f_j(\mathbf{w}_j) = \frac{1}{2L_j} \left\| (\mathbf{w}_j * \psi_j * \mathbf{x}) - \sqrt{L_j}(\psi_j * \mathbf{y}) \right\|_2^2 \quad (14) $$

and the regularization terms $r_j(\mathbf{w}_j) = \lambda_j \|\mathbf{w}_j\|_2$. The multi-band robust deconvolution (Multi-RDCP) cost function at $\mathbf{W} = (\mathbf{w}_j)_{j=0}^{J-1}$ is given by:

$$ F_{\Psi}(\mathbf{W}) = \sum_{j=0}^{J-1} f_j(\mathbf{w}_j) + r_j(\mathbf{w}_j). \quad (15) $$

*Note.* The double convolutions above, $(\mathbf{w}_j * \psi_j * \mathbf{x})$, form a hybrid filterbank [18]; with dyadic subsampling, it would resemble a multiresolution neural network (MuReNN) [19].

*Note.* The term $\sum_{j=0}^{J-1} \lambda_j \|\mathbf{w}_j\|_2$ corresponds to weighted Group LASSO [7] where the groups are the subbands of $\Psi$.

In order to synthesize the global $\mathbf{w}^*$ from $\mathbf{W}$, we use the formula motivated by the following proposition.

**Proposition III.4.** Given $\mathbf{x}, \mathbf{y} \in \mathbb{R}^N$, let $\mathbf{W}^* = (\mathbf{w}_j^*)_{j=0}^{J-1}$ be a global optimum of unregularized Multi-RDCP (Def. III.3), i.e. $\lambda_j = 0$ for $j \in [J]$. Then

$$ \mathbf{w}^* = \sum_{j=0}^{J-1} \frac{1}{\sqrt{L_j}} \mathbf{C}_{\psi_j}^\top \mathbf{C}_{\psi_j} \mathbf{w}_j^* \quad (16) $$

is a global optimum for the deconvolution of $\mathbf{y}$ from $\mathbf{x}$.

*Proof.* At the global optimum, one has $\nabla F_{\Psi}(\mathbf{W}^*) = 0$, thus:

$$ L_j \nabla f_j(\mathbf{w}_j^*) = \mathbf{C}_{(\psi_j * \mathbf{x})}^\top (\mathbf{C}_{(\psi_j * \mathbf{x})} \mathbf{w}_j^* - \sqrt{L_j} \mathbf{C}_{\psi_j}^\top \mathbf{y}) = 0 $$

for all $j \in [J]$. The convolution is associative and commutative:

$$ \mathbf{C}_{(\psi_j * \mathbf{x})}^\top \mathbf{C}_{(\psi_j * \mathbf{x})} = \mathbf{C}_{\psi_j}^\top \mathbf{C}_\mathbf{x}^\top \mathbf{C}_{\psi_j} \mathbf{C}_\mathbf{x} = \mathbf{C}_{\mathbf{x}}^\top \mathbf{C}_\mathbf{x} \mathbf{C}_{\psi_j}^\top \mathbf{C}_{\psi_j}. $$

Thus, by dividing each identity by $\sqrt{L_j}$ and summing over $j$:

$$ \sum_{j=0}^{J-1} \mathbf{C}_{\mathbf{x}}^\top \mathbf{C}_{\mathbf{x}} \left( \frac{1}{\sqrt{L_j}} \mathbf{C}_{\psi_j}^\top \mathbf{C}_{\psi_j} \mathbf{w}_j^* \right) = \sum_{j=0}^{J-1} \mathbf{C}_{\mathbf{x}}^\top \mathbf{C}_{\psi_j}^\top \mathbf{C}_{\psi_j} \mathbf{y}. $$

We recognize $\mathbf{C}_{\mathbf{x}}^\top \mathbf{C}_{\mathbf{x}} \mathbf{w}^*$ on the left and $\mathbf{C}_{\mathbf{x}}^\top \mathbf{y}$ on the right (see Theorem III.2). Hence: $\mathbf{C}_{\mathbf{x}}^\top (\mathbf{C}_{\mathbf{x}} \mathbf{w} - \mathbf{y}) = 0$, i.e., $\nabla f(\mathbf{w}^*) = 0$ (see Equation 2). We conclude by convexity of the MSE. □

As the objective of Multi-RDCP is separable in $(\mathbf{w}_j)_{j=0}^{J-1}$, it can be solved with PGD for each $\mathbf{w}_j$ for $j \in [J]$, since it corresponds to the cost function of a rescaled RDCP.

**Proposition III.5.** For $j \in [J]$, $F_j = f_j + r_j$ is the cost function of RDCP (Def. II.1) for some given $\mathbf{x}_j \in \mathbb{R}^N$ and $\mathbf{y}_j \in \mathbb{R}^N$. For an initial guess $\mathbf{w}_j^{(0)} \in \text{ran}(\mathbf{C}_\mathbf{x}^\top)$ and the step size $\alpha = L_j/S_{\text{max}}(\psi_j * \mathbf{x})$, the iterates $\{\mathbf{w}_j^{(t)}\}_{t \in \mathbb{N}}$ of PGD applied to $F_j$ satisfy for every $j \in [J]$:

$$ \frac{F_j(\mathbf{w}^{(t)}) - F_j^*}{F_j(\mathbf{w}^{(0)}) - F_j^*} \le \left(1 - \frac{S_{\text{min}}^+(\psi_j * \mathbf{x})}{S_{\text{max}}(\psi_j * \mathbf{x})}\right)^t. \quad (17) $$

*Proof.* For $\mathbf{y}_j = L_j^{1/2}(\psi_j * \mathbf{y})$ and $\mathbf{x}_j = L_j^{-1/2}(\psi_j * \mathbf{x})$ we see that $f_j(\mathbf{w}_j) = \frac{1}{2L_j} \|(\mathbf{w}_j * \psi_j * \mathbf{x}) - \sqrt{L_j}(\psi_j * \mathbf{y})\|_2^2 = \frac{1}{2} \|\mathbf{w}_j * \mathbf{x}_j - \mathbf{y}_j\|_2^2$. Theorem II.14 concludes this result. □

*Note.* Given the result above, we set $L_j = S_{\text{max}}(\psi_j * \mathbf{x})$ in order to have a unit step size of $\alpha = 1$.

Algorithm 1 summarizes our method: across the different subbands $\psi_j$, we normalize adaptively with respect to the given signal $\mathbf{x}$ and derive kernel estimates $\mathbf{w}_j$ (SNAKE).

**Theorem III.6.** Given $\mathbf{x}, \mathbf{y} \in \mathbb{R}^N$, a Parseval filterbank $\Psi = (\psi_0, \dots, \psi_{J-1})$ and Multi-RDCP (Def. III.3), the iterates $\{\mathbf{W}^{(t+1)}\}_{t \in \mathbb{N}}$ of SNAKE (Algorithm 1) satisfy:

$$ \frac{F_{\Psi}(\mathbf{W}^{(t)}) - F_{\Psi}^*}{F_{\Psi}(\mathbf{W}^{(0)}) - F_{\Psi}^*} \le \left(1 - \min_{j \in [J]} \frac{S_{\text{min}}^+(\psi_j * \mathbf{x})}{S_{\text{max}}(\psi_j * \mathbf{x})}\right)^t. \quad (18) $$
**Algorithm 1** Subband-Normalized Adaptive Kernel Estimation
Given $\mathbf{x}$ and $\mathbf{y} = \mathbf{h} * \mathbf{x} + \mathbf{\eta}$ where $\mathbf{\eta}$ is noise, finds $\mathbf{w} \approx \mathbf{h}$.

Require: $(\psi_j)_{j=0}^{J-1}$ (Parseval), $(\lambda_j)_{j=0}^{J-1}$ positive, $T \in \mathbb{N}$

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

2: $\mathbf{x}_j = \psi_j * \mathbf{x}$

3: $L_j = S_{\max}(\mathbf{x}_j)$ (see Equation (8))

4: $\nabla f_j : \mathbf{w}_j \mapsto \frac{1}{L_j}(\mathbf{w}_j * \mathbf{x}_j - \sqrt{L_j}(\psi_j * \mathbf{y})) * \bar{\mathbf{x}}_j$

5: $\text{Prox}_{r_j} : \mathbf{v} \mapsto \max(0, 1 - \lambda_j / \|\mathbf{v}\|_2)\mathbf{v}$

6: $\mathbf{w}_j^{(0)} = \mathbf{0}$

7: **for** $t = 0 : T - 1$ **do**

8: $\mathbf{w}_j^{(t+1)} = \text{Prox}_{r_j}(\mathbf{w}_j^{(t)} - \nabla f_j(\mathbf{w}_j^{(t)})$,

9: **return** $\mathbf{w} = \sum_{j=0}^{J-1} \frac{1}{\sqrt{L_j}} \mathbf{w}_j^{(T)} * \psi_j * \bar{\psi}_j$

*Proof.* We write $F_{\Psi}(\mathbf{W}) = \sum_{j=0}^{J-1} F_j(\mathbf{w}_j)$. Proposition III.5 gives us, for each $j \in [J]$:

$$F_j(\mathbf{w}_j^{(t)}) - F_j^* \le \left(1 - \frac{S_{\min}^+(\psi_j * \mathbf{x})}{S_{\max}(\psi_j * \mathbf{x})}\right)^t (F_j(\mathbf{w}^{(0)}) - F_j^*) \\ \le \left(1 - \min_{j \in [J]} \frac{S_{\min}^+(\psi_j * \mathbf{x})}{S_{\max}(\psi_j * \mathbf{x})}\right)^t (F_j(\mathbf{w}^{(0)}) - F_j^*)$$

Summing over $j \in [J]$ concludes the result. $\square$

<figure><img src="image_2.png" /><figcaption>Figure 1. Plot of relative error to optimal cost function value of $F$ (as in LHS of (12)) for PGD (green) and of $F_{\Psi}$ (as in LHS of (18)) for SNAKE (blue) as a function of CPU time. We report the median and quartiles across 100 independent samples of $\mathbf{x}$ and $\mathbf{\eta}$.</figcaption></figure>

This convergence rate does not depend on the ratio of the smallest to largest positive power spectrum entry $S_{\min}^{+}(\mathbf{x})/S_{\max}^{+}(\mathbf{x})$ like PGD for RDCP, but on the ratios within subbands. This represents an improvement, even with a simple filterbank, as we will demonstrate in the next section.

## IV. NUMERICAL SIMULATION

We consider a kernel $\mathbf{h} \in \mathbb{R}^N$, $N = 256$ (see Figure 1, left). The input $\mathbf{x}$ is a sample of Brownian noise [20]. We observe $\mathbf{x}$ and $\mathbf{y} = \mathbf{h} * \mathbf{x} + \mathbf{\eta}$, where $\mathbf{\eta}$ is an sample of white Gaussian noise with $\|\mathbf{\eta}\|_2 = \frac{1}{20} \|\mathbf{x}\|_2$. Note that $\mathbf{h}$ is narrowband while $\mathbf{\eta}$ is broadband. In the absence of regularization, deconvolution of $\mathbf{y}$ from $\mathbf{x}$ would yield an estimate $\tilde{\mathbf{w}}$ such that $\tilde{\mathbf{w}} * \mathbf{x} \approx \mathbf{y}$. Instead, by promoting sparsity across subbands of the filterbank $\Psi$, multi-RDCP aims at $\mathbf{w}^*$ such that $\mathbf{w}^* \approx \mathbf{h}$.

We build a Parseval filterbank $\tilde{\Psi} = (\psi_j)_{j=0}^{J-1}$ of $J = 7$ “real Shannon wavelets” [21]. For every $j \in [J]$ and $k \le N/2$, $\hat{\psi}_j[k] = 1$ if $2^j < k \le 2^{j+1}$ and 0 elsewhere. For $k > N/2$, we complete the construction by Hermitian symmetry.

After a preliminary stage of hyperparameter calibration, we select $\lambda = 40 \|\mathbf{x}\|_2$ for RDCP and $\lambda_j = 150/L_j \|\mathbf{x}\|_2$ for Multi-RDCP. We leave the important question of sensitivity to hyperparameters to future work.

The CPU time per iteration is roughly $J = 7$ times greater for SNAKE than for PGD. Yet, this discrepancy is compensated by the faster convergence guarantee of SNAKE: see Theorem III.6. Figure 1 illustrates this phenomenon.

To assess the quality of our reconstructions and justify our choice of regularization, we have employed an oracle benchmark based on Ridge Regression, where the optimal solution $\mathbf{w}^*$ is computed with full knowledge of the ground truth $\mathbf{h}$. Ridge Regression is defined by the regularizer $r(\mathbf{w}) = \lambda \|\mathbf{w}\|_2^2$ and since the optimality condition is $(\mathbf{w} * \mathbf{x} - \mathbf{y}) * \bar{\mathbf{x}} + \frac{\lambda}{2} \mathbf{w} = 0$, this admits a closed-form solution via the FFT. For each

<figure><img src="image_3.png" /><figcaption>Figure 2. Left: a narrowband signal $\mathbf{h}$ (green) and its unregularized deconvolution $\tilde{\mathbf{w}}$ (gray) given Brownian noise $\mathbf{x}$ and $\mathbf{y} = \mathbf{h} * \mathbf{w} + \mathbf{\eta}$ where $\mathbf{\eta}$ is white noise. Right: typical solutions $\mathbf{w}^*$ or $\mathbf{w}^{(T)}$ from regularized deconvolutions. For clarity, we have zoomed onto the red rectangle.</figcaption></figure>

instance of $\mathbf{x}$ and $\mathbf{\eta}$, we test across different $\lambda > 0$, to see what the best possible $\mathbf{w}^*$ achievable by Ridge Regression is, while measuring the reconstruction error as $\|\mathbf{w}^* - \mathbf{h}\|_2 / \|\mathbf{h}\|_2$. Figure 2 illustrates our results.

Over 100 instances, the best possible Ridge Regression reconstruction has a median error of 6.7%. In comparison, the reconstructions of RDCP achieves a median error of 6.7% and Multi-RDCP attains 2.3%.

## V. CONCLUSION

We have presented two contributions: a new regularization technique named multiband robust deconvolution (Multi-RDCP) and a dedicated algorithm named Subband-Normalized Adaptive Kernel Evaluation (SNAKE). SNAKE benefits from filterbank analysis in Multi-RDCP to adapt the step size of subband-wise proximal gradient descents to local spectral characteristics, hence an improved convergence rate. Furthermore, Multi-RDCP promotes sparsity across the subband spectrum of the solution, surpassing the best possible reconstruction performance of Ridge Regression.

Looking ahead, SNAKE’s adaptability opens possibilities for other regularization schemes, where classical methods struggle due to high spectral variability. Additionally, exploring adaptive strategies for selecting the subband regularization parameters $(\lambda_j)_{j=0}^{J-1}$ based on prior knowledge of the target signal holds promise to further enhance performance.
REFERENCES

[1] J. Kovacevic, V. K. Goyal, and M. Vetterli, *Fourier and Wavelet Signal Processing*. Cambridge University Press, 2014.

[2] S. M. Riad, "The deconvolution problem: An overview," *Proceedings of the IEEE*, vol. 74, no. 1, pp. 82–85, 1986.

[3] T. Ceccherini-Silberstein, F. Scarabotti, and F. Tolli, *Discrete Harmonic Analysis: Representations, Number Theory, Expanders, and the Fourier Transform*. Cambridge University Press, 2018.

[4] A. V. Oppenheim and R. W. Schafer, *Discrete-Time Signal Processing*, 2nd ed. Upper Saddle River, NJ: Prentice Hall, 1999.

[5] G. Fant, *Speech Acoustics and Phonetics*. Kluwer Academic Publishers, 2004.

[6] N. H. Fletcher and T. D. Rossing, *The physics of musical instruments*. Springer Science & Business Media, 2012.

[7] M. Yuan and Y. Lin, "Model selection and estimation in regression with grouped variables," *Journal of the Royal Statistical Society: Series B (Statistical Methodology)*, vol. 68, no. 1, pp. 49–67, 2006.

[8] A. Beck and M. Teboulle, "A fast iterative shrinkage-thresholding algorithm for linear inverse problems," *SIAM Journal on Imaging Sciences*, vol. 2, no. 1, pp. 183–202, 2009.

[9] J.-J. Moreau, "Fonctions convexes duales et points proximaux dans un espace hilbertien," *Comptes Rendus de l'Académie des Sciences de Paris*, vol. A255, no. 22, pp. 2897–2899, Nov. 1962.

[10] N. Parikh and S. Boyd, "Proximal algorithms," *Foundations and Trends in Optimization*, vol. 1, no. 3, pp. 127–239, Jan. 2014.

[11] H. Karimi, J. Nutini, and M. Schmidt, "Linear convergence of gradient and proximal-gradient methods under the Polyak–Łojasiewicz condition," in *Proceedings of the European Conference on Machine Learning and Principles and Practice of Knowledge Discovery in Databases (ECML PKDD)*. Springer, 2016, pp. 795–811.

[12] B. T. Polyak, "Gradient methods for the minimisation of functionals," *USSR Computational Mathematics and Mathematical Physics*, vol. 3, no. 4, pp. 864–878, 1963.

[13] S. Łojasiewicz, "Une propriété topologique des sous-ensembles analytiques réels," in *Les équations aux dérivées partielles*, vol. 117. Colloques internationaux du CNRS, 1963, pp. 87–89.

[14] R. A. Horn and C. R. Johnson, *Matrix Analysis*, 2nd ed. Cambridge University Press, 2012.

[15] R. M. Gray, *Toeplitz and Circulant Matrices: A Review*. Foundations and Trends in Communications and Information Theory, 2006, vol. 2, no. 3.

[16] R. T. Rockafellar and R. J.-B. Wets, *Variational Analysis*, ser. Grundlehren der mathematischen Wissenschaften. Berlin: Springer-Verlag, 1998, vol. 317.

[17] P. G. Casazza and G. Kutyniok, *Finite Frames: Theory and Applications*. Springer Science & Business Media, 2012.

[18] D. Haider, F. Perfler, V. Lostanlen, M. Ehler, and P. Balazs, "Hold me tight: Stable encoder-decoder design for speech enhancement," in *Proceedings of the International Speech Communication Association Conference (INTERSPEECH)*. ISCA, 2024.

[19] V. Lostanlen, D. Haider, H. Han, M. Lagrange, P. Balazs, and M. Ehler, "Fitting auditory filterbanks with multiresolution neural networks," in *Proceedings of the IEEE Workshop on Applications of Signal Processing to Audio and Acoustics*, 2023.

[20] C. Gardiner, *Stochastic methods*. Springer Berlin Heidelberg, 2009, vol. 4.

[21] M. W. Frazier, *An Introduction to Wavelets through Linear Algebra*. Springer Science & Business Media, 2006.