Identifiability Conditions for Multi-channel Blind Deconvolution with Short Filters

TL;DR

This paper establishes conditions for the well-posedness of multi-channel blind deconvolution with short filters and demonstrates an effective nuclear norm minimization approach through numerical experiments.

Contribution

It derives necessary and sufficient conditions for identifiability and reformulates the problem as a low-rank matrix recovery task solved by nuclear norm minimization.

Findings

Effective algorithm under certain generative models

Successful in noiseless and noisy scenarios

Provides theoretical identifiability conditions

Abstract

This work considers the multi-channel blind deconvolution problem under the assumption that the channels are short. First, we investigate the ill-posedness issues inherent to blind deconvolution problems and sufficient and necessary conditions on the channels that guarantee well-posedness are derived. Following previous work on blind deconvolution, the problem is then reformulated as a low-rank matrix recovery problem and solved by nuclear norm minimization. Numerical experiments show the effectiveness of this algorithm under a certain generative model for the input signal and the channels, both in the noiseless and in the noisy case.

Equations

$${\mathbf{y}}_n=\mathbf{s}⊛{\mathbf{w}}_n,n\in [N],$$

$${\mathbf{w}}_n=\mathbf{C}{\mathbf{h}}_n=\left[\begin{array}{c}{\mathbf{h}}_n\\ 0_{L-K}\end{array}\right],{\mathbf{h}}_n\in {\mathbb{R}}^K,n\in [N],$$

$$\mathbf{F}{\mathbf{y}}_n=\sqrt{L}\mathbf{F}\mathbf{s}\odot \mathbf{F}\mathbf{C}{\mathbf{h}}_n,n\in [N],$$

$$\mathcal{B}(\hat{\mathbf{p}},\{{\mathbf{q}}_n{\}}_{n=0}^{N-1})=\{\sqrt{L}\hat{\mathbf{p}}[l]{\hat{\mathbf{q}}}_n[l]|(l,n)\in [L]\times [N]\},$$

$$f(\mathbf{x})=\frac{1}{2}{\Vert \mathcal{B}(\mathbf{x})-\hat{\mathbf{y}}\Vert }_2^2.$$

$${\nabla }^{2}g({\mathbf{x}}_0)={\mathbf{A}}^{\ast }\mathbf{A},$$

$$\mathbf{J}=\left[\begin{array}{ccccc}{\mathbf{D}}_{{\hat{\mathbf{w}}}_0}& {\mathbf{D}}_{\hat{\mathbf{s}}}{\mathbf{F}}_K& 0& \cdots & 0\\ {\mathbf{D}}_{{\hat{\mathbf{w}}}_1}& 0& {\mathbf{D}}_{\hat{\mathbf{s}}}{\mathbf{F}}_K& \cdots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ {\mathbf{D}}_{{\hat{\mathbf{w}}}_{N-1}}& 0& 0& \cdots & {\mathbf{D}}_{\hat{\mathbf{s}}}{\mathbf{F}}_K\end{array}\right]$$

$$\text{find}\mathbf{X}\text{s.t.}\hat{\mathbf{y}}=\mathcal{A}(\mathbf{X})\text{and}\mathrm{rank}(\mathbf{X})=1.$$

$$\underset{\mathbf{X}}{\mathrm{minimize}}$$

$$\underset{\mathbf{p},\mathbf{q}}{\mathrm{minimize}}$$

$$\mathcal{L}(\mathbf{p},\mathbf{q},\mathbf{\lambda },\sigma )=\frac{1}{2}\left({\Vert \mathbf{p}\Vert }_2^2+{\Vert \mathbf{q}\Vert }_2^2\right)-⟨\mathbf{\lambda },\mathcal{A}(\mathbf{p}{\mathbf{q}}^{\ast })-\hat{\mathbf{y}}⟩+\frac{\sigma }{2}{\Vert \mathcal{A}(\mathbf{p}{\mathbf{q}}^{\ast })-\hat{\mathbf{y}}\Vert }_2^2,$$

$$\left\lVert\boldsymbol{X}_{0}-\boldsymbol{p}\boldsymbol{q}^{*}\right\rVert_{F}\big{/}\left\lVert\boldsymbol{X}_{0}\right\rVert_{F},$$

$$\boldsymbol{v}_{n}=\sigma\cdot\left\lVert\boldsymbol{y}_{n}\right\rVert_{2}\cdot\boldsymbol{\nu}_{n}\big{/}\left\lVert\boldsymbol{\nu}_{n}\right\rVert_{2},$$

$$f({\mathbf{x}}_0+ϵ\mathbf{v})=f({\mathbf{x}}_0)+ϵ\nabla f({\mathbf{x}}_0{)}^T\mathbf{v}+\frac{{ϵ}^2}{2}{\mathbf{v}}^{\ast }{\nabla }^{2}f({\mathbf{x}}_0)\mathbf{v}+\cdots .$$

$$f({\mathbf{x}}_0+ϵ\mathbf{v})=0\text{or}{\mathbf{v}}^{\ast }{\nabla }^{2}f({\mathbf{x}}_0)\mathbf{v}=0$$

$${\sum }_{n=0}^{N-1}{\Vert {\mathbf{D}}_{{\hat{\mathbf{w}}}_n}\hat{\mathbf{p}}+{\mathbf{D}}_{\hat{\mathbf{s}}}{\mathbf{F}}_K{\mathbf{q}}_n\Vert }_2^2=0$$

$$\mathbf{v}=[-{\hat{\mathbf{s}}}^{\ast },{\mathbf{h}}_0^{\ast },\dots ,{\mathbf{h}}_{N-1}^{\ast }{]}^{\ast }.$$

$${\mathbf{x}}_0-{\mathbf{x}}^{\prime }=[{\hat{\mathbf{s}}}^{\ast }(1-{\alpha }^{-1}),{\mathbf{h}}_0^{\ast }(1-\alpha ),\dots ,{\mathbf{h}}_{N-1}^{\ast }(1-\alpha ){]}^{\ast }.$$

$${\mathbf{x}}_0-{\mathbf{x}}^{\prime }=[ϵ{\hat{\mathbf{s}}}^{\ast },-ϵ{\mathbf{h}}_0^{\ast },\dots ,-ϵ{\mathbf{h}}_{N-1}^{\ast }{]}^{\ast }=-ϵ\mathbf{v}$$

$${\mathbf{D}}_{{\hat{\mathbf{w}}}_n}^{-1}{\mathbf{D}}_{\hat{\mathbf{s}}}{\mathbf{F}}_K{\mathbf{q}}_n={\mathbf{D}}_{{\hat{\mathbf{w}}}_m}^{-1}{\mathbf{D}}_{\hat{\mathbf{s}}}{\mathbf{F}}_K{\mathbf{q}}_m,(n,m)\in [N]\times [N].$$

$${\mathbf{D}}_{{\hat{\mathbf{w}}}_m}{\mathbf{F}}_K{\mathbf{q}}_n={\mathbf{D}}_{{\hat{\mathbf{w}}}_n}{\mathbf{F}}_K{\mathbf{q}}_m,(n,m)\in [N]\times [N]$$

$${\mathbf{w}}_m⊛\mathbf{C}{\mathbf{q}}_n={\mathbf{w}}_n⊛\mathbf{C}{\mathbf{q}}_m,(n,m)\in [N]\times [N].$$

$$\mathcal{Z}({\mathbf{w}}_m)\mathcal{Z}(\mathbf{C}{\mathbf{q}}_n)=\mathcal{Z}({\mathbf{w}}_n)\mathcal{Z}(\mathbf{C}{\mathbf{q}}_m)$$

$$P_{{\mathbf{h}}_m}(z)P_{{\mathbf{q}}_n}(z)=P_{{\mathbf{h}}_n}(z)P_{{\mathbf{q}}_m}(z),(n,m)\in [N]\times [N].$$

$$P_{{\mathbf{q}}_n}(z)=\alpha (z-\beta )P_{{\mathbf{h}}_n}(z),n\in [N]$$

$$P_{{\mathbf{q}}_n}(z)=\alpha \frac{P_{{\mathbf{h}}_n}(z)}{z-\beta }(z-\gamma ),n\in [N]$$

$$\mathrm{Pr}({\mathbf{h}}_n[K-1]=0)\stackrel{a.s.}{=}0,n\in [N]$$

Full text

1 Introduction

The objective of multi-channel blind deconvolution is to recover unknown $L$-length vectors $\boldsymbol{s}$ and $\{\boldsymbol{w}_{n}\}_{n=0}^{N-1}$ from their circular111In blind deconvolution problems, circular convolution is generally preferred to linear convolution so that the convolution theorem for the Discrete Fourier Transform (DFT) applies. Even though convolutions arising in practical scenarios are not circular, a linear convolution can reasonably be approximated by a circular convolution if the support of the filter is sufficiently short with respect to the length of the input signal (i.e., if the filter decays quickly) [1]. The assumption that the convolution is circular is thus not limiting in most practical situations. convolutions

[TABLE]

where $[N]=\{0,\dots,N-1\}$. This problem has a lot of applications in many fields of science and engineering. In communications systems for example, blind deconvolution is referred to as blind equalization [2]. In this context, $\boldsymbol{s}$ is a message that needs to be transmitted to a receiver having $N$ antennas whereas $\boldsymbol{w}_{n}$ is the channel from the transmitter to the $n$-th receiving antenna. Traditionally, the transmitter would first send a known pilot message that allows the receiver to estimate the channels. The receiver would then use these estimates to invert the effect of the channels on subsequent messages. If the channels are varying over time (e.g., in mobile communications), this process needs to be periodically repeated, causing a non-negligible overhead. Performing blind channel equalization would allow to avoid this overhead.

Multi-channel blind deconvolution can also be used in noise imaging applications [3]. One such application in the field of geophysical imaging is “seismic-while-drilling” (SWD). In SWD, the objective is to perform subsurface imaging on a drilling site without disrupting drilling operations. One proposed way of achieving this is to use the vibrations generated by the drill-bit in the subsurface as the input signal $\boldsymbol{s}$. Multiple receivers recording the seismic traces at several locations would then allow to recover the “channel impulse responses of the Earth” $\{\boldsymbol{w}_{n}\}_{n=0}^{N-1}$ using multi-channel blind deconvolution [4].

Other applications of blind deconvolution include medical imaging [5], astronomical imaging [6] and computer vision [7].

1.1 Problem setup

Without any additional assumptions, the multi-channel blind deconvolution problem formulated above counts more unknowns than observations. The problem is thus under-determined and ill-posed. In this work, the channels (also often called filters) are assumed to live in the following $K$-dimensional subspace

[TABLE]

where $\boldsymbol{C}$ consists in the first $K$ columns of the $L\times L$ identity matrix. In other words, the channels are assumed to be short or supported on $[K]$. The input signal $\boldsymbol{s}$ is simply assumed to be an unknown $L$-length random vector. Recovering $\boldsymbol{s}$ and $\{\boldsymbol{w}_{n}\}_{n=0}^{N-1}$ is then equivalent to recovering $\boldsymbol{s}$ and $\{\boldsymbol{h}_{n}\}_{n=0}^{N-1}$ and a mandatory condition for successful recovery reads ${LN\geq L+KN}$, i.e., the information count must be favorable. The problem setup is illustrated in Fig. 1.

1.2 Related work and contributions

Blind deconvolution being an ubiquitous problem in science and engineering, there exists a vast amount of literature on the topic. This work was primarily inspired by [8, 9]. In [8], Ahmed, Recht and Romberg consider the single-input single-channel blind deconvolution problem. Both the input signal and the channel are assumed to live in known subspaces, one of which is random. In [9], Ahmed and Demanet leverage multiple diverse inputs – living in known random subspaces – to weaken the assumption on the channel; the subspace in which the channel lives is now also to be discovered. In this work, in a sense “dual” to the one studied in [9], we leverage multiple diverse and short channels to free the input signal from any assumption. It is interesting to note that no random subspaces are involved in this work. This more realistic setting is challenging as it prevents us to use the classical stability proofs based on, e.g., high dimensional statistics and measure concentration theory. The short channel assumption was already used in [10] for the single-input single-channel case and in [11, 12] for the multi-channel case. However, [12] makes extra-assumptions only valid for the specific case of underwater acoustic deconvolution.

In [11, 8, 9, 12], the problem is reformulated as a low-rank matrix recovery problem from linear measurements and solved via nuclear norm minimization using Burer and Monteiro’s algorithm [13, 14, 15]. However, [11] only shows “example” of successful recovery but lacks a detailed characterization of the performances of the algorithm, thus giving little insight on the feasibility of the problem. Robustness to noise is not discussed neither. In this work, we also use Burer and Monteiro’s algorithm. Yet, and at the opposite of prior works using this algorithm to solve blind deconvolution problems, we do not convexify the problem according to the work in [15] and use a heuristic strategy to deal with local minima instead. This also seems to be the approach followed in [12], although no strategy to deal with local minima is described. The first contribution of this work is to assess the performances of Burer and Monteiro’s algorithm to solve multi-channel blind deconvolution with short channels.

Making structural assumptions on the channels is not sufficient to guarantee that the problem is well-posed. In [16], necessary and sufficient identifiability conditions are derived for multi-channel blind deconvolution with short filters. The second contribution of this work is to give a new interpretation of those conditions by linking them to a localized ill-posedness through a local analysis – around the ground truth – of the $\ell_{2}$-norm misfit between the observations and a candidate solution.

1.3 Organization of the paper

Sec. 2 first investigates the ill-posedness issues arising in multi-channel blind deconvolution with short filters. Necessary and sufficient conditions on the channels that make the problem well-posed are derived. A class of channels satisfying the aforementioned conditions almost surely is also given. Next, Sec. 3 reformulates the multi-channel blind deconvolution problem as a low-rank matrix recovery problem and describes Burer and Monteiro’s algorithm and the heuristic strategy used to avoid getting trapped at local minima. Numerical experiments then assess the performances of this algorithm under a generative model for the input signal and the channels that ensures that the problem is well-posed. Both the noiseless and the noisy cases are investigated. Finally, Sec. 4 concludes this work by discussing the obtained results, the limits of the short channel assumption, and some perspectives for future research.

2 Local ill-posedness and identifiability conditions

As discussed in the introduction, assuming $K$-length channels and having $LN\geq L+KN$ is not sufficient to guarantee that the problem (1) is well-posed. Indeed, let $\boldsymbol{f}$ be a filter that admits an inverse $\boldsymbol{f}^{-1}$ such that $\boldsymbol{f}\circledast\boldsymbol{f}^{-1}=\boldsymbol{\delta}$. If $\boldsymbol{f}$ and $\boldsymbol{w}_{n}$ are such that $\boldsymbol{w}_{n}^{\prime}=\boldsymbol{w}_{n}\circledast\boldsymbol{f}$ is still supported on $[K]$ for all $n\in[N]$, then $\{\boldsymbol{w}^{\prime}_{n}\}_{n=0}^{N-1}$ and $\boldsymbol{s}^{\prime}=\boldsymbol{s}\circledast\boldsymbol{f}^{-1}$ also constitutes a valid solution of the multi-channel blind deconvolution problem. This ill-posedness will be referred to as convolutive ambiguity. As a particular case, $\{\alpha\boldsymbol{w}_{n}\}_{n=0}^{N-1}$ and $\alpha^{-1}\boldsymbol{s}$ also constitutes a valid solution for any scalar $\alpha\neq 0$. While this scalar ambiguity is acceptable in most applications, convolutive ambiguity is not. In this section, we derive sufficient and necessary conditions on the channels that guarantee that the scalar ambiguity is the only possible ambiguity, making the problem well-posed.

We start by looking at the observed circular convolutions in the Fourier domain. Noting $\boldsymbol{F}\in\mathbb{C}^{L\times L}$ the unitary DFT matrix and using the convolution theorem, (1) becomes

[TABLE]

where $\odot$ is the entry-wise product. We then define a bilinear map $\mathcal{B}:\mathbb{C}^{L}\times\mathbb{R}^{KN}\to\mathbb{C}^{LN}$ mapping any input signal in the Fourier domain $\hat{\boldsymbol{p}}\in\mathbb{C}^{L}$ and any set of $K$-length filters $\{\boldsymbol{q}_{n}\}_{n=0}^{N-1}$ to their corresponding observed circular convolutions in the Fourier domain, i.e.,

[TABLE]

where $\hat{\boldsymbol{q}}_{n}=\boldsymbol{F}\boldsymbol{C}\boldsymbol{q}_{n}$. Eq. (2) then reads $\hat{\boldsymbol{y}}=\mathcal{B}(\hat{\boldsymbol{s}},\{\boldsymbol{h}_{n}\}_{n=0}^{N-1})$ with $\hat{\boldsymbol{y}}$ the concatenation of all the observed convolutions in the Fourier domain and $\hat{\boldsymbol{s}}=\boldsymbol{F}\boldsymbol{s}$. To ease notations, let $\boldsymbol{q}$ be the concatenation of all the filters $\{\boldsymbol{q}_{n}\}_{n=0}^{N-1}$ and $\boldsymbol{x}=[\hat{\boldsymbol{p}}^{*},\boldsymbol{q}^{*}]^{*}$ be the concatenation of all the arguments of $\mathcal{B}$. The natural way to solve the multi-channel blind deconvolution problem is then to minimize the $\ell_{2}$-norm data misfit given by

[TABLE]

First, Lemma 1, proved in Appendix A, makes the link between “local ill-posedness” and the Hessian of the objective function $f$ around the ground truth.

Lemma 1**.**

Let $\boldsymbol{x}_{0}$ be the ground truth. For sufficiently small $\epsilon$, $\boldsymbol{x}_{0}+\epsilon\boldsymbol{v}$ also minimizes the objective function if $\boldsymbol{v}$ belongs to the null space of $\boldsymbol{\nabla}^{2}f(\boldsymbol{x}_{0})\in\mathbb{C}^{(L+KN)\times(L+KN)}$.

For ease of notation, the null space of the matrix $\boldsymbol{\nabla}^{2}f(\boldsymbol{x}_{0})$ is noted $\mathcal{N}$. Lemma 1 tells us that characterizing $\mathcal{N}$ is equivalent to characterizing the local ill-posedness of the problem around the ground truth. The following known lemma [17] is the first step toward the characterization of $\mathcal{N}$.

Lemma 2**.**

Let $\mathcal{A}$ be a differentiable operator such that $\mathcal{A}(\boldsymbol{x}_{0})=\boldsymbol{b}$. Let also $g(\boldsymbol{x})=\frac{1}{2}\left\lVert\mathcal{A}(\boldsymbol{x})-\boldsymbol{b}\right\rVert^{2}_{2}$. The Hessian of $g(\boldsymbol{x})$ at $\boldsymbol{x}_{0}$ is then given by

[TABLE]

where $\boldsymbol{A}=\frac{\delta\mathcal{A}(\boldsymbol{x})}{\delta\boldsymbol{x}}\rvert_{\boldsymbol{x}_{0}}$ is the linearization of $\mathcal{A}(\boldsymbol{x})$ around $\boldsymbol{x}_{0}$ (i.e., the Jacobian of $\mathcal{A}$ at $\boldsymbol{x}_{0}$).

Using Lemma 2, we then have $\boldsymbol{\nabla}^{2}f(\boldsymbol{x}_{0})=\boldsymbol{J}^{*}\boldsymbol{J}$ with ${\boldsymbol{J}\in\mathbb{C}^{LN\times(L+KN)}}$ the Jacobian of $\mathcal{B}$ at $\boldsymbol{x}_{0}$ given by

[TABLE]

where $\boldsymbol{D}_{\boldsymbol{x}}=\operatorname{diag}(\boldsymbol{x})$ and $\boldsymbol{F}_{K}\in\mathbb{C}^{L\times K}$ is the restriction of $\boldsymbol{F}$ to its $K$ first columns. The next lemma, proved in Appendix B, gives us a first insight on the structure of $\mathcal{N}$.

Lemma 3**.**

The kernel $\mathcal{N}$ is always at least one-dimensional and contains the scalar ambiguity inherent to any blind deconvolution problem.

The scalar ambiguity is related to Lemma 3 as it imposes the existence of a vector $\boldsymbol{v}\in\mathcal{N}$ such that, if ${\boldsymbol{x}_{0}=[\hat{\boldsymbol{s}}^{*},\boldsymbol{h}^{*}]^{*}}$ is the ground truth, then, for $\epsilon$ sufficiently small, ${\boldsymbol{x}^{\prime}=[(1+\epsilon)^{-1}\hat{\boldsymbol{s}}^{*},(1+\epsilon)\boldsymbol{h}^{*}]^{*}\approx\boldsymbol{x}_{0}+\epsilon\boldsymbol{v}}$ is also a valid solution, with $\boldsymbol{v}=[-\hat{\boldsymbol{s}}^{*},\boldsymbol{h}^{*}]^{*}$ by a first order approximation (see Appendix B). For the scalar ambiguity to be the only possible ambiguity, we thus would like $\mathcal{N}$ to be exactly one-dimensional. From the rank-nullity theorem, a first necessary condition is to have $LN\geq L+KN-1$. The next proposition, proved in Appendix C, gives two sufficient and necessary conditions on the filters for $\mathcal{N}$ to be exactly one-dimensional.

Proposition 1**.**

Assume that the signal $\boldsymbol{s}$ has no zero in the Fourier domain, i.e., $\hat{\boldsymbol{s}}[l]\neq 0$ for all $l\in[L]$, and $LN\geq L+KN-1$. The kernel $\mathcal{N}$ is one-dimensional if and only if

there exists an index $n^{\prime}\in[N]$ such that $\boldsymbol{h}_{n^{\prime}}[K-1]\neq 0$, 2. 2.

the polynomials $\{\sum_{k=0}^{K-1}\boldsymbol{h}_{n}[k]z^{k}\}_{n=0}^{N-1}$ are coprime, i.e., they do not share any common root.

In particular, Prop. 1 implies that the support of the channels must be complete on “both sides” to avoid ill-posedness issues. Indeed, the first condition prevents situations where $\boldsymbol{h}_{n}[K-1]=0$ for all $n\in[N]$, i.e., situations where the “end” of the support is not filled by any channel. Furthermore, the second condition prevents, in particular, situations where $\boldsymbol{h}_{n}[0]=0$ for all $n\in[N]$ (as 0 would then be a shared root), i.e., situations where the “beginning” of the support is not filled by any channel. As a consequence, the support size $K$ needs to be exactly known to avoid ill-posedness issues. This limits our short channels model to situations where the support size of the channels can be reliably estimated beforehand.

As mentioned in the introduction, similar sufficient and necessary conditions were already derived in [16]. To our knowledge, this is, however, the first time that those conditions are linked to $\mathcal{N}$ and to a notion of “local ill-posedness”. In [16], the authors also derive a condition involving the input signal $\boldsymbol{s}$. It is not obvious how this condition might appear from the above reasoning.

The conditions given in Prop. 1 are verified by specific filters. The next proposition, proved in Appendix D, determines a class of random vectors $\{\boldsymbol{h}_{n}\}_{n=0}^{N-1}$ that satisfy the two conditions of Prop. 1 almost surely.

Proposition 2**.**

$K$-length vectors $\{\boldsymbol{h}_{n}\}_{n=0}^{N-1}$ whose entries are independent continuous random variables whose distribution can be defined by density functions satisfy the conditions of Prop. 1 almost surely.

3 Algorithm and numerical experiments

Following the approach described in [9], this section reformulates the multi-channel blind deconvolution problem and describes the algorithm used to solve the obtained new problem.

The observed convolutions (1) can be written in the Fourier domain as $\hat{\boldsymbol{y}}=\mathcal{A}(\boldsymbol{s}\boldsymbol{h}^{*})$ where $\boldsymbol{h}$ is the concatenation of all the unknown filters and $\mathcal{A}$ is a linear operator. Multi-channel blind deconvolution can thus be recasted as a low-rank matrix recovery problem under linear equality constraints as follows

[TABLE]

This problem is however non-convex and NP-hard. Fortunately, it can be relaxed as a convex program using the nuclear norm heuristic [13]

[TABLE]

where $\left\lVert\boldsymbol{X}\right\rVert_{*}$ is the sum of the singular values of $\boldsymbol{X}$. Because $\boldsymbol{X}$ is known to be a rank-1 matrix, we then write $\boldsymbol{X}=\boldsymbol{p}\boldsymbol{q}^{*}$ with $\boldsymbol{p}$ and $\boldsymbol{q}$ the candidate solution for the input signal and the filters, respectively. Program (3) can be shown to be equivalent to

[TABLE]

except that the latter is non-convex and thus subject to local minima. In [14], Burer and Monteiro proposed an algorithm based on the method of multipliers to efficiently find a local minimum of (4). Starting from an initial guess, this algorithm iteratively minimizes the following augmented Lagrangian with respect to $\boldsymbol{p}$ and $\boldsymbol{q}$ [13, 14]

[TABLE]

where $\boldsymbol{\lambda}\in\mathbb{C}^{LN}$ is the vector of Lagrange multipliers and $\sigma$ is a positive penalty parameter. A way of updating $\boldsymbol{\lambda}$ and $\sigma$ at each iteration is also proposed in [14]. The algorithm stops when $\left\lVert\mathcal{A}(\boldsymbol{p}\boldsymbol{q}^{*})-\hat{\boldsymbol{y}}\right\rVert^{2}_{2}$ is below a given tolerance. In this work, the minimization is performed using L-BFGS with the MATLAB solver minFunc [18].

To avoid getting trapped at local minima and find a global minimum instead, the following heuristic strategy was implemented. When $\left\lVert\mathcal{A}(\boldsymbol{p}\boldsymbol{q}^{*})-\hat{\boldsymbol{y}}\right\rVert^{2}_{2}$ stopped decreasing for a given number of iterations and is largely above the tolerance, the algorithm decrees to be trapped at a local minimum. At this point, the algorithm simply starts again from another initial guess. Following the work of Burer and Monteiro [15], (4) can also be convexified by writing $\boldsymbol{X}=\boldsymbol{P}\boldsymbol{Q}^{*}$ with $\boldsymbol{P}\in\mathbb{R}^{L\times r}$, $\boldsymbol{Q}\in\mathbb{R}^{KN\times r}$ and $r>1$. In this case, the algorithm stops when either $\boldsymbol{P}$ or $\boldsymbol{Q}$ is rank deficient and the recovered vectors $\boldsymbol{p}$ and $\boldsymbol{q}$ are given by the leading singular vectors of $\boldsymbol{P}$ and $\boldsymbol{Q}$, respectively. This is for example the approach used in [8] for the single-input single-channel case. However, we empirically observed better performances with the first approach. Interestingly, this is also the approach followed in [12] to solve multi-channel blind deconvolution for underwater acoustic, although no strategy to handle the local minima is described.

In the absence of a better strategy, $\boldsymbol{p}$ and $\boldsymbol{q}$ are initialized as independent standard random Gaussian vectors. The deterministic nature of the linear operator $\mathcal{A}$ does not allow us to use a smarter initialization scheme such as, for example, the spectral method proposed in [10].

A sketch of the complete algorithm is given in Algorithm 1.

3.1 Numerical experiments

In the following numerical experiments, both the input signal $\boldsymbol{s}$ and the channels $\{\boldsymbol{h}_{n}\}_{n=0}^{N-1}$ are independent random Gaussian vectors width i.i.d. entries. In this way, the identifiability conditions of Prop. 1 are satisfied almost surely. To comply with the existing literature on blind deconvolution using low-rank matrix recovery, the metric used to assess the quality of the recovery is the relative outer product error, i.e,

[TABLE]

where $\boldsymbol{X}_{0}=\boldsymbol{s}\boldsymbol{h}^{*}$ is the ground truth, $\boldsymbol{p}$ is the recovered input signal and $\boldsymbol{q}$ is the concatenation of the recovered filters.

3.1.1 Noiseless case

We first assess the performances of the proposed algorithm in the noiseless case. Figures 2a and 2b depict the empirical probability of successful recovery and the average number of attempts before successful recovery (i.e., the number of times the algorithm has to restart because it got stuck at a local minimum) when $L=32$, $N$ is ranging from 2 to 10 and $K$ is ranging from 1 to 32. The recovery is considered to be successful if the relative outer product error is less than222This threshold has been chosen to comply with the existing literature. Most of the time, however, the actual relative error is orders of magnitude smaller. 2%. For each pair of parameters $N$ and $K$, the value displayed in Figs. 2a and 2b is the result of an average over 100 deconvolutions, each time involving new random instances of $\boldsymbol{s}$ and $\{\boldsymbol{h}_{n}\}_{n=0}^{N-1}$. On both figures, a red curve indicates the information-theoretic limit. Above this curve, the problem is under-determined and successful recovery is thus not possible. Below this curve, however, Fig. 2a shows that the recovery is successful with probability 1 almost everywhere. This result is particularly interesting when compared to “usual” results in blind deconvolution where the minimal number of observations is larger than the number of unknowns (usually within log factors). Next, Fig. 2b shows that the average number of attempts before successful recovery is usually small and decreases when the oversampling factor increases. The decision to not convexify (4) along with the heuristic strategy implemented to avoid getting trapped at local minima thus seem reasonable.

3.1.2 Robustness to noise

Each observed convolution $\boldsymbol{y}_{n}$ is now corrupted by an additive noise $\boldsymbol{v}_{n}$ given by

[TABLE]

where $\boldsymbol{\nu}_{n}\in\mathbb{R}^{L}$ is a standard Gaussian random vector. The SNR is then simply given by $\text{SNR}_{n}=-20\log_{10}\sigma$ and is assumed to be the same for each observed convolution. Fig. 3 depicts the average relative outer product error against the SNR. Each data point is the result of an average over 2800 deconvolutions. For each deconvolution, new random instances of $\boldsymbol{s}$, $\{\boldsymbol{h}_{n}\}_{n=0}^{N-1}$ and the noise were used. Except at low SNR (i.e., less than 10dB), the relative error scales linearly with the noise level in a log-log scale as desired. In addition, Fig. 3 shows that increasing the oversampling factor increases the robustness to noise.

4 Conclusion

Section 2 revisited known identifiability conditions for the channels and linked those conditions with a notion of “local ill-posedness” around the ground truth. It is not yet clear, however, how identifiability conditions involving the input signal could be derived with the same approach. This section also highlighted a limitation of the short channels model used in this work; the support size $K$ needs to be exactly known for the problem to be well-posed. In most applications, however, the support size is not known a priori but rather needs to be discovered. In such applications, other ways of encoding the short channels assumption should be employed (see for example [4]). What happens when the conditions of Prop. 1 are “almost not satisfied” (e.g., when the support is “almost not filled” or when the polynomials $\{\sum_{k=0}^{K-1}\boldsymbol{h}_{n}[k]z^{k}\}_{n=0}^{N-1}$ “almost share a root”) is another interesting question. Preliminary numerical experiments suggest that it degrades well-posedness. However, how to quantitatively measure this degradation remains an open question.

Next, Sec. 3.1 demonstrated the effectiveness of the proposed algorithm to solve multi-channel blind deconvolution under a generative model for the channels that ensures that the problem is well-posed. It would also be interesting to study the impact of channels similarities on the performances of the algorithm. To this effect, a measure of channels similarities and a way of controlling it must first be devised. The condition number of the $K\times N$ matrix whose columns are the channels might be a good candidate to measure channels similarities, as suggested by unreported experiments.

Acknowledgment

The authors would like to thank Laurent Demanet, Ali Ahmed and Pawan Bharadwaj from MIT and Luc Vandendorpe from UCLouvain for interesting discussions. Laurent Jacques is funded by the F.R.S.-FNRS.

Appendix A Proof of Lemma 1

We write the Taylor series of the objective function around the ground truth $\boldsymbol{x}_{0}$

[TABLE]

Assuming $\epsilon$ is small enough for the higher order terms to be negligible and because $f(\boldsymbol{x}_{0})=0$ and $\boldsymbol{\nabla}f(\boldsymbol{x}_{0})=0$, points $\boldsymbol{x}_{0}+\epsilon\boldsymbol{v}$ such that

[TABLE]

also minimize the objective function. The set of solutions of this last equation is given by the null space of $\boldsymbol{\nabla}^{2}f(\boldsymbol{x}_{0})$.

Appendix B Proof of Lemma 3

The null space of $\boldsymbol{\nabla}^{2}f(\boldsymbol{x}_{0})$ can be computed by solving (5). Letting $\boldsymbol{v}=[\hat{\boldsymbol{p}}^{*},\boldsymbol{q}^{*}_{0},\dots,\boldsymbol{q}_{N-1}^{*}]^{*}$, this can be developed as

[TABLE]

or, equivalently $\boldsymbol{D}_{\hat{\boldsymbol{w}}_{n}}\hat{\boldsymbol{p}}=-\boldsymbol{D}_{\hat{\boldsymbol{s}}}\boldsymbol{F}_{K}\boldsymbol{q}_{n}$ for all $n\in[N]$. The obvious solution to this system of equations is given by ${\hat{\boldsymbol{p}}=-\alpha\hat{\boldsymbol{s}}}$ and $\boldsymbol{q}_{n}=\alpha\boldsymbol{h}_{n}$ for all $n\in[N]$ and for any scalar $\alpha$. This gives us a basis vector of the null space of $\boldsymbol{\nabla}^{2}f(\boldsymbol{x}_{0})$

[TABLE]

We now want to prove that the space generated by this vector contains the scalar ambiguity inherent to any blind deconvolution problem. Let’s consider the ground truth $\boldsymbol{x}_{0}$ and another minimizer of the objective function resulting from the scalar ambiguity, $\boldsymbol{x}^{\prime}=[\alpha^{-1}\hat{\boldsymbol{s}}^{*},\alpha\boldsymbol{h}_{0}^{*},\dots,\alpha\boldsymbol{h}_{N-1}^{*}]^{*}$ for any scalar $\alpha\neq 0$. The vector joining $\boldsymbol{x}^{\prime}$ and $\boldsymbol{x}_{0}$ is given by

[TABLE]

For $\alpha=1+\epsilon$ and $\epsilon$ close to 0, this reduces to333Using a first order Taylor series around 0, $\frac{\epsilon}{1+\epsilon}\approx\epsilon$.

[TABLE]

and $\boldsymbol{x}^{\prime}=\boldsymbol{x}_{0}+\epsilon\boldsymbol{v}$.

Appendix C Proof of Prop. 1

The objective is to prove that the obvious solution to (6) is the only one. For other solutions to exist, we would need $\hat{\boldsymbol{p}}=-\boldsymbol{D}_{\hat{\boldsymbol{w}}_{n}}^{-1}\boldsymbol{D}_{\hat{\boldsymbol{s}}}\boldsymbol{F}_{K}{\boldsymbol{q}}_{n}$ for all $n\in[N]$ which then implies

[TABLE]

With the assumption that $\hat{\boldsymbol{s}}[l]\neq 0$ for $l\in[L]$, this gives

[TABLE]

which in the time domain becomes

[TABLE]

We now have to prove that $\boldsymbol{q}_{n}=\alpha\boldsymbol{h}_{n}$ for all $n\in[N]$ is the only solution to this equation if and only if conditions 1 and 2 are true. Taking the $\mathcal{Z}$-transform on both sides of (7) and using the convolution theorem leads to

[TABLE]

for all pairs $(n,m)\in[N]\times[N]$. After the inconsequential change of variable $z^{-1}\to z$, and noting $P_{\boldsymbol{x}}(z)=\sum_{k=0}^{K-1}\boldsymbol{x}[k]z^{k}$, this becomes

[TABLE]

The objective has now been transformed into proving that (8) only admits $P_{\boldsymbol{q}_{n}}(z)=\alpha P_{\boldsymbol{h}_{n}}(z)$ for all $n\in[N]$ as a solution if and only if conditions 1 and 2 are true.

Let us first prove the forward direction by contraposition. We want to prove that if the negation of either the first or the second condition is true, then there exists other solutions to (6) and $\mathcal{N}$ is thus more than one dimensional.

We first assume that the negation of the first condition is true: there does not exist any index $n\in[N]$ such that $\boldsymbol{h}_{n}[K-1]\neq 0$. Equivalently, there does not exist any index $n\in[N]$ such that $P_{\boldsymbol{h}_{n}}(z)$ is of degree $>K-2$. In this case,

[TABLE]

also constitutes a solution to (8) for any $\beta\in\mathbb{C}$ and the obvious solution is thus not the only one. 2. 2.

We now suppose that the negation of the second condition is true: there is a root $\beta\in\mathbb{C}$ shared by all the polynomials $P_{\boldsymbol{h}_{n}}(z)$. Then

[TABLE]

also constitutes a solution to (8) for any $\gamma\in\mathbb{C}$. For $\gamma\neq\beta$, this solution is different from the obvious one and the obvious solution is thus not the only one.

Next, concerning the backward direction, we want to prove that if the two conditions are true, then the only solution to (6) is the obvious one and $\mathcal{N}$ is one-dimensional. A necessary condition for (8) to be satisfied is that $P_{\boldsymbol{h}_{m}}(z)P_{\boldsymbol{q}_{n}}(z)$ and $P_{\boldsymbol{h}_{n}}(z)P_{\boldsymbol{q}_{m}}(z)$ have the same roots for all $(n,m)\in[N]\times[N]$.

We first assume that the second condition is satisfied: the polynomials $\{P_{\boldsymbol{h}_{n}}(z)\}_{n=0}^{N-1}$ do not share any common root. In this case, the necessary condition for (8) to be satisfied directly translates to $P_{\boldsymbol{q}_{n}}(z)$ having at least all the roots of $P_{\boldsymbol{h}_{n}}(z)$, $\forall n\in[N]$. This is however not sufficient to conclude that $P_{\boldsymbol{q}_{n}}(z)=\alpha P_{\boldsymbol{h}_{n}}(z)$ for all $n\in[N]$ as we could for example add a root $\beta\in\mathbb{C}$ in every $P_{\boldsymbol{q}_{n}}(z)$ if the first condition is not true (as was done in the preceding paragraph). 2. 2.

Next, we suppose that the first condition is also satisfied: there exists an index $n^{\prime}\in[N]$ such that $\boldsymbol{h}_{n^{\prime}}[K-1]\neq 0$. Equivalently, there exists an index $n^{\prime}\in[N]$ such that $P_{\boldsymbol{h}_{n^{\prime}}}(z)$ is of degree $K-1$. In this case, it is no longer possible to add a root $\beta\in\mathbb{C}$ in every $P_{\boldsymbol{q}_{n}}(z)$ as it would result in $P_{\boldsymbol{q}_{n^{\prime}}}(z)$ being of degree $K$, which is not possible.

We can conclude that, if the two conditions are true, $P_{\boldsymbol{q}_{n}}(z)$ and $P_{\boldsymbol{h}_{n}}(z)$ must have exactly the same roots for all $n\in[N]$. Because polynomials having exactly the same roots are equal up to a scalar factor, we can conclude that $P_{\boldsymbol{q}_{n}}(z)=\alpha P_{\boldsymbol{h}_{n}}(z)$ for all $n\in[N]$ is the only solution to (6) and that $\mathcal{N}$ is one dimensional and only contains the scalar ambiguity inherent to blind deconvolution.

Appendix D Proof of Prop. 2

Let us first prove that the first condition of Prop. 1 is satisfied with probability one. For all $n\in[N]$, $\boldsymbol{h}_{n}[K-1]$ is a continuous random variable whose distribution can be defined by a density function. As such, one can write

[TABLE]

and the first condition is satisfied almost surely. For the second condition, we work with the set of polynomials $\{P_{\boldsymbol{h}_{n}}(z)\}_{n=0}^{N-1}$. A proof that the probability for two such polynomials to have any common root is zero is given in [19, p. 442]. This proof can be extended to any number of polynomials and the second condition is thus also satisfied with probability one.