Phase Retrieval for Sparse Signals: Uniqueness Conditions

I Introduction

In many real-world scenarios, we naturally measure the Fourier transform (FT) of a signal of interest instead of the signal itself. During the measuring process, it may happen that the phase of the FT is lost or irremediably distorted. The recovery of the phase is fundamental to reconstructing the signal, and this recovery process is known as phase retrieval (PR). Phase loss problems occur in many scientific fields, particularly those involving optics and communications. For example, in X-ray crystallography, the measurements are the diffraction patterns of a crystallized molecule and we would like to recover the molecule itself.

Although the PR problem has a long history with a rich literature , there are still open questions regarding the uniqueness of the solution and the existence of reliable algorithms to recover the signal. We underline that in many applications the signal of interest is sparse: for example, the atoms of a molecule are distinct elements in the spatial domain. However, a review of the main results related to PR reveals that sparsity has only been rarely exploited to obtain uniqueness conditions or efficient reconstruction algorithms. Moreover, we note that the problem is usually defined on the discrete domain for simplicity, while most of the applications involve continuous signals.

In this paper, we present a uniqueness condition for 1–dimensional sparse signals exploiting a previous result related to the turnpike problem, a classic combinatorial problem. We show that the same uniqueness condition holds for multidimensional signals. Note that these results are valid both for discrete and continuous domains, the condition relies solely on the sparsity and on the characteristics of the support of the signal. The only difference between the discrete and the continuous problem is the probability of satisfying the uniqueness condition.

In what follows, we set the notation and we precisely state the PR problem for continuous signals. We then describe a number of applications, emphasizing the role of sparsity in these scenarios, and we show that this property can be further exploited to obtain a uniqueness condition.

II Problem statement and applications

In this section, we state the phase retrieval for signals defined on $D$–dimensional continuous domains. We underline the difficulties characterizing these problems and we define the non-trivial concept of unique solution. We introduce a sparse model for continuous signals and we present a number of applications that can exploit such model. Unless otherwise stated, we use the following notation:

bold lower case symbols, such as $\boldsymbol{x}$, for vectors,

bold capital symbols, such as $\boldsymbol{A}$, for matrices,

$x_{n}$ is the $n$-th element of the vector $\boldsymbol{x}$, $A_{m,n}$ is the element in the $m$-th row and the $n$-column of $\boldsymbol{A}$,

capital calligraphic letters, such as $\mathcal{X}$, for sets.

We define the FT of the signal $f(\boldsymbol{x})$ as

where $\phi(\boldsymbol{u})$ is the phase of the FT. Then, we define the PR problem as follows: given the magnitude $|\widehat{f}(\boldsymbol{u})|$ of the FT, recover the original signal $f(\boldsymbol{x})$. Since the FT is a bijective mapping, the problem is equivalent to recovering the phase term $\phi(\boldsymbol{u})$, hence the name phase retrieval. It is easy to show that the knowledge of $|\widehat{f}(\boldsymbol{u})|$ is equivalent to the knowledge of the autocorrelation function (ACF), defined as

More precisely, the ACF is the IFT of $|\widehat{f}(\boldsymbol{u})|^{2}$. Now, we have all the ingredients to state the PR problem for a continuous signal $f(\boldsymbol{x})$.

PR for continuous signals Given the magnitude $|\widehat{f}(\boldsymbol{u})|$ of the FT or the ACF $a(\boldsymbol{x})$ of a signal of interest $f(\boldsymbol{x})$, recover the signal itself.

We can show that in the general case the solution of Problem 1 is not unique. In fact, we can simply assign a random phase to the measured magnitudes. Moreover, some information on the signal $f(\boldsymbol{x})$ is entirely embedded into the phase of its FT and we cannot hope to recover this information from the magnitude only. Namely, the following transformations of $f(\boldsymbol{x})$,

do not influence the magnitude of the FT. Hence, time–reversal, sign change and absolute position cannot be recovered once the phase is lost. It is then appropriate to define what “unique” means with the following equivalence class

II-B Sparse signals

A natural question arises from Problem 1: “What are the conditions that $f(\boldsymbol{x})$ must satisfy to have a unique PR?”

In this paper, we constrain the PR problem using a sparse model for $f(\boldsymbol{x})$. More precisely, we define a $N$–sparse signal using the Dirac delta notation, that is

where the $n$-th delta has coefficient $c^{(n)}$ and is located at $\boldsymbol{x}^{(n)}$ and $N$ is finite. Then, the ACF is defined by the following linear combination of $N^{2}-N+1$ deltas,

where $d^{(n)}$ and $\boldsymbol{y}^{(n)}$ are the coefficients and the locations of the deltas in the ACF. Note that the ACF is centro-symmetric, meaning that for every delta located at $\boldsymbol{y}^{(n)}$, there is another one with the same coefficient located at $-\boldsymbol{y}^{(n)}$. Then, we can rewrite the ACF as

where we can consider only the second sum instead of the whole ACF, since it contains all the available information.

II-C Applications

While (2) may look too simple to model signals of interest for real-world applications, the following three scenarios are of interest and fit this model.

The is the primary technique to determine the structure of molecules. The experiment consists of the following steps: first, the molecule $e(\boldsymbol{x})$ of interest is crystallized. The obtained crystal $f(\boldsymbol{x})$ is simply a periodic repetition of the basic structure $e(\boldsymbol{x})$, called unit cell,

where $\boldsymbol{N}$ is a vector containing the sizes of the unit cell in each dimension. Second, the crystal is exposed to an X-ray beam, under different angles. For each angle, we have a diffraction pattern that, mathematically speaking, is a slice of the three–dimensional FT of the crystal. See Figure 1 for a graphical depiction of a unit cell, a crystal and a diffraction pattern. The diffraction patterns are recorded using traditional imaging techniques, such as CCDs, and only the magnitude is acquired. Hence, we aim at recovering the spatial distribution $e(\boldsymbol{x})$ of the molecule, called electron density, from the magnitude of the FT of $f(\boldsymbol{x})$.

Due to the periodicity of the crystal, we are in fact measuring the magnitude of the Fourier series coefficients of the crystal, $|\widehat{f}_{\boldsymbol{u}}|$. This set of coefficients is equal, up to a constant factor, to samples of the magnitude of the unit cell FT, $|e(\boldsymbol{x})|$.

The PR problem in X-ray crystallography is more complex than Problem 1. In fact, the set of measured samples is not sufficiently dense to reconstruct $|\widehat{e}(\boldsymbol{u})|^{2}$ using Shannon’s sampling theorem. More precisely, we have an undersampling factor of two for each dimension and we do not dispose of the entire ACF of $e(\boldsymbol{x})$. See for more details about this remark.

There is some a-priori information about the crystal that we can exploit. For example, we can realistically model the unit cell $e(\boldsymbol{x})$ as

where $\phi(\boldsymbol{x})$ is the electron density of a single atomIn reality, every atom has a different electron density $\phi(\boldsymbol{x})$. However, this assumption is reasonable and many reconstruction methods used in crystallography consider similar assumptions . that has a positive coefficient $c^{(n)}$ and is located at $\boldsymbol{x}^{(n)}$. We can now specify the PR problem for crystallography.

PR for Crystallography Consider a unit cell $e(\boldsymbol{x})$ with $N$ positive atoms on a bounded domain. Given a set of magnitudes of the Fourier series coefficients $|\widehat{f}_{\boldsymbol{u}}|$, estimate the locations and amplitudes of the atoms.

II-C2 Speckle imaging in astronomy

Another example of the PR problem can be found in astronomy, namely an imaging method known as speckle imaging. This technique attempts to mitigate the resolution downgrade introduced by atmospheric turbulences. Namely, the atmosphere blurs $M$ images $\{g^{(i)}(\boldsymbol{x})\}_{i=1}^{M}$ collected by a telescope and the blurring is modeled as a linear filter that may vary for each image, $\{s^{(i)}(\boldsymbol{x})\}_{i=1}^{M}$. The $i$-th measured image, also called speckle, is the convolution between the astronomic object $f(\boldsymbol{x})$ and the $i$-th linear filter

See Figure 2 for an example of the speckles $g^{\,(i)}(\boldsymbol{x})$ and the target of the astronomic observations $f(\boldsymbol{x})$.

We reduce the atmospheric distortion and some potential additive white noise by taking the average of the squared magnitudes of FT of the images as

where ${\widehat{g}}^{\,(i)}(\boldsymbol{u})$, $\widehat{f}(\boldsymbol{u})$ and $\widehat{s}^{\,(i)}(\boldsymbol{u})$ are the FT of the measured image, the object of interest, and the transfer function of the atmosphere, respectively. Note that the averaging strategy is effective since we assume that the atmospherical transfer functions generally affects only the phases of $\widehat{f}(\boldsymbol{u})$ . The averaged atmospherical transfer function $\sum_{k=1}^{M}|\widehat{s}^{\,(i)}(\boldsymbol{u})|^{2}$ is estimated using atmospheric models or images of a reference astronomical object.

PR is necessary to recover the high resolution image of the astronomic object of interest from $|\widehat{f}(\boldsymbol{u})|^{2}$. Note that we introduced the problem considering a continuous model for both the astronomical object and the measured images. However, since the images are generally measured and processed as sampled data, we may consider a set of discrete images $\{g^{(i)}_{\boldsymbol{k}}\}_{i=1}^{M}$ and the DFT is usually employed.

We assume that $f(\boldsymbol{x})$ is sparse,

since the astronomical object is composed of a set of stars that can be modeled as Dirac deltas. Therefore, we have the following PR problem for speckle imaging.

II-C3 Blind channel estimation

We conclude this section with an interesting example of 1–dimensional PR: blind estimation of a communication channel. The knowledge of the channel impulse response is fundamental for wireless communications systems, such as the ones based on Orthogonal Frequency Division Multiplexing . According to the theory of multi-path propagation , the channel $g(t)$ can be faithfully modeled as

where $t$ is the time variable and the $N$ deltas describe the multi-path phenomenon. More precisely, each delta represents a secondary communication path generated by a reflective body between the source and the receiver.

We would like to estimate the locations and the coefficients of the deltas, without having direct control of the channel input, hence the “blind estimation” terminology. We only measure samples of the output $y(t)$ of the channel, that is the convolution between the input $f(t)$ and the channel itself,

The input data $f(t)$, which is the result of the modulation of a discrete sequence $x_{n}$, is usually whitened to achieve the maximum capacity of the channel. If the input sequence is statistically white, then the magnitude of the FT of the output signal is in expectation equal to the magnitude of the FT of the channel. Once more, we cannot access the continuous-time output $y(t)$, but only a set of $M$ samples $\{y_{k}\}_{k=1}^{M}$. To recover the channel $g(t)$, we can take the DFT of the collected samples, keep the magnitudes and solve the following PR problem.

PR for Blind Channel Estimation Let $g(t)$ be a multi-path fading communication channel as defined in (5), where $N$ is finite and generally small. Assume the input of the channel to be properly whitened. Then, estimate the channel impulse response $g(t)$ from a set of samples $\{y_{k}\}_{k=1}^{M}$ of the output $y(t)$.

We conclude this list of applications emphasizing the leitmotif connecting all the different applications: we are interested in PR for $N$–sparse signal $f(\boldsymbol{x})$ defined on a $D$–dimensional continuous domain. We collect samples of the signal and the phase information is lost, as in blind channel estimation, or irreparably distorted, as in speckle imaging. We would like to recover the sparse components in the continuous domain, without discretizing the solution’s domain.

This approach already proved beneficial in other domains. For example, it has been shown that it is possible to recover a $N$-sparse signal $f(x)$ from only $2N+1$ samples of the filtered signal $(g*f)(x)$, see . Another example where the continuous-time model has been proven to be effective is in channel estimation. In , the authors demonstrated that the channel estimator based on the continuous-time model achieves better performance when compared to the state-of-the-art discrete approaches.

III Literature review

We present a literature review that covers theoretical and algorithmic results for both continuous and discrete PR. Note that most of the works focused on the latter, given the difficulties of treating the PR for continuous signals.

Most of the relevant works connected to the continuous sparse PR problem was developed in combinatorics for the turnpike problem . The turnpike problem deals with the recovery of the locations of a set of points from their unlabeled distances. Note that the recovery of the support of $f(\boldsymbol{x})$ from the support of $a(\boldsymbol{x})$ is an instance of such problem. A theorem presented by Piccard in 1939 gives a sufficient condition for the uniqueness of the turnpike problem. Unfortunately, a counterexample to the theorem was first found by Bloom et al. and its generalization was recently obtained by Bekir et al. . A similar but weaker condition for multidimensional signals has been recently obtained by Senechal in 2008 . Skiena et al. proposed a non-trivial algorithm for solving the problem. It is known as the backtracking algorithm and solves any instance of the turnpike problem providing the existence of (possibly multiple) valid solutions. The algorithm has a polynomial computational complexity when the set $\{\boldsymbol{x}^{(n)}\}$ is drawn at random. Zhang showed how to build sets of points achieving the worst case computational complexity, that is $\mathcal{O}(2^{n}n\log n)$.

An equivalent problem has been stated for restriction site mapping, an interesting task in computational molecular biology, where a particular enzyme is added to a DNA sample, so that the DNA is cut at particular locations $\{x^{(n)}\}_{n=1}^{N}$, known as restriction sites. One can find the distance between each pair of restriction sites, $\{x^{(n)}-x^{(m)}\}_{n,m=1}^{N}$, using gel electrophoresis. Given the distances, we would like to recover the locations of the sites. This technique is used for DNA mapping and it usually involves different enzymes. When a single enzyme is used, it is known as partial digest . Note that it has been shown by Cieliebak et al. that the partial digest problem with noisy measurements is $\mathsf{NP}$-hard. The partial digest problem is in fact a turnpike problem with integer locations—a bridge between continuous and discrete PR problems.

III-B Discrete PR

The uniqueness of the discrete PR problem has been studied for multidimensional discrete signals. One of the main results is given by Hayes : the set of positive finitely supported images $f_{\boldsymbol{k}}$ which are not uniquely recoverable has measure zero. The results are derived using the theory of multidimensional polynomials. A possible algorithm to recover signals from the magnitudes of the FT is also given, but it does not achieve satisfactory results according to the authors. Note that this uniqueness result cannot be directly applied to X-ray crystallography (see Remark 1).

On the algorithmic front, many reconstruction algorithms were developed for Problem 2 and a review is given in . Among them, ab-initio or direct methods were introduced in the late 50s and have the considerable advantage of not requiring any prior information regarding the crystals. The state of the art among these type of algorithms is charge flipping . It performs two operations iteratively, one in the spatial domain where it imposes the positiveness of the electron density and the bounded support, one in the Fourier domain where it imposes the measured magnitudes. This algorithm was first presented in 2004 , while some of the recent developments are described in . It can be seen as a version of the Gerchberg-Saxton algorithm , where the positivity constraint is enforced if the electron density is above a certain threshold or set to zero otherwise.

Recent work defined efficient convex relaxations for solving discrete PR problems. These approaches are potentially more stable w.r.t noise and are capable of avoiding local minima. This strategy has been introduced by Lu et al. , who described a necessary condition for the uniqueness of PR for discrete 1–dimensional sparse signals together with a reconstruction algorithm based on the lifting of the problem in a higher dimensional linear space that is solved by traditional convex optimization solvers. A similar approach, but introducing random masks to improve the redundancy of measurements, has been introduced by Candès et al., . Hassibi and his collaborators proposed an improved algorithm and sufficient probabilistic uniqueness conditions based on the sparsity of the signal of interest. Waldspurger et al. formulated another tractable convex relaxation similar to the classical MaxCut semidefinite program that achieves better reconstruction performance when compared to other convex relaxations.

PR has been generalized to any linear operator beyond the FT. More precisely, we choose a frame $\left\{\mathcal{F}_{i}\right\}_{i=1}^{K}$ and we collect measurements of an unknown signal $f$ using the elements of this frame $\left\{w_{i}\right\}_{i=1}^{K}=\left\{\left<\mathcal{F}_{i},f\right>\right\}_{i=1}^{K}$. Equivalently to the PR problem, we assume we can only rely on the magnitude of the measurements and obtain a generalized PR problem: from the magnitude of the expansion $\left\{|w_{i}|\right\}_{i=1}^{K}$, recover the original signal $f$. Note that if the chosen frame is the Fourier frame, then we have an instance of PR on a discrete domain. Balan et al. formulated the problem and studied different theoretical and algorithmic aspect of the recovery of signals from the magnitude of generic frame coefficients. Specifically, fast algorithms are given in , while the statement of equivalent problems and the construction of particular frames for which the reconstruction is unique are given in , respectively. Relevant uniqueness results are given by Chebira et al. , where a necessary and sufficient condition for uniqueness has been described, however it requires exponential time to be checked. Note that the aforementioned convex relaxations can be applied to this generalized PR.

Other generalization of the PR problem have been considered. Oppenheim et al. proved the uniqueness of PR problems up to the knowledge of the signs of the Fourier coefficients in . A general analysis of the phase loss and the magnitude loss is given in and an extension to the multidimensional case is given in . The main result concerns the uniqueness and the reconstruction of the magnitude loss problem.

IV Uniqueness of the sparse PR problem

We use a divide and conquer approach to derive the uniqueness condition. First, we notice that the locations of the deltas of the ACF contain more information than their coefficients. In fact, if all the deltas have the same coefficient, the coefficients of the ACF do not carry any information. Therefore, we consider the problem of recovering the support of $f(\boldsymbol{x})$ given the support of its ACF $a(\boldsymbol{x})$. We then use the coefficients to further restrict the possibility of having a non-unique solution.

We consider the set of locations $\left\{\boldsymbol{x}^{(n)}\right\}_{n=1}^{N}$ and derive the set of differences $\mathcal{D}=\left\{\boldsymbol{x}^{(n)}-\boldsymbol{x}^{(m)}\right\}_{n,m=1}^{N,N}$. Given the possibility of repeated elements, $\mathcal{D}$ is formally a multiset. Looking at (3), we also notice that all the elements of a sparse ACF are supported on the set of differences, that is $\boldsymbol{y}^{(n)}\in\cal{D}$.

If we attempt to recover the support of a sparse signal from the support of its ACF, we realize that the lack of labeling of the elements of $\cal{D}$ makes the problem combinatorial, that is all the possible labelings must be tested to find the optimal solution. Moreover, the solution may be even more complex if the ACF has “collisions”: two deltas of the ACF located at the same position due to two couples of equi-spaced deltas in the signal $f(\boldsymbol{x})$.

We say there is a collision in the ACF when $\exists\,n\neq m\text{ such that }\boldsymbol{y}^{(n)}=\boldsymbol{y}^{(m)}$.

In other words, let $\boldsymbol{x}^{(i)}$,$\boldsymbol{x}^{(j)}$,$\boldsymbol{x}^{(k)}$ and $\boldsymbol{x}^{(l)}$ be the locations of four distinct deltas of a sparse signal $f(\boldsymbol{x})$, then we say that we have a collision in the ACF if $\boldsymbol{x}^{(i)}-\boldsymbol{x}^{(j)}=\boldsymbol{x}^{(k)}-\boldsymbol{x}^{(l)}$.

The main reason why collisions are problematic is the impossibility of knowing a-priori how many $\boldsymbol{y}^{(n)}$ are colliding on the same element of the ACF that we observe. In what follows, we show that if the observed ACF $a(\boldsymbol{x})$ does not have collisions, we are able to recover uniquely the sparse signal $f(\boldsymbol{x})$ in most of the cases.

We assume that we first want to recover the support of $f(\boldsymbol{x})$ from the support of $a(\boldsymbol{x})$. Here we study the problem for a one–dimensional signal $f(x)$. We consider the set of differences $\mathcal{D}$ as the input and the locations of the deltas $\left\{x^{(n)}\right\}_{n=1}^{N}$ as the output of the following problem:

Support Recovery Given all the pairwise distances $\{x^{(n)}-x^{(m)}\}_{n,m=1}^{N,N}$ between a set of $N$ points lying on a 1–dimensional domain, recover their locations $\{x^{(n)}\}_{n=1}^{N}$.

Note that we have no information about the labeling of the pairwise differences in $\cal{D}$. Therefore, Problem 5 is combinatorial and is equivalent to an instance of the turnpike problem . If the naming was known, the problem could be easily solved by multidimensional scaling .

We introduce the definition of homometric sets to define the uniqueness of Problem 5.

Two sets $\mathcal{X}$ and $\mathcal{Y}$ are said to be homometric if and only if their difference sets are congruent, that is $\mathcal{D}_{\mathcal{X}}=\mathcal{D}_{\mathcal{Y}}$.

In this section, we assume that all the differences are different from each other, i.e. we do not have any collision in the ACF. This is equivalent to saying that the set $\mathcal{D}$ has no repeated elements. It turns out that the problem of the support recovery was first posed by Patterson and a possible solution was proposed by Piccard in 1939 . More precisely, Piccard suggested that if there are no collisions, the solution of the turnpike problem is always unique.

Unfortunately, a counterexample to this result was found in 1975 by Bloom . Consider $\mathcal{X}=\left\{0,1,4,10,12,17\right\}$ and $\mathcal{Y}=\left\{0,1,8,11,13,17\right\}$, then

Recently, this counterexample has been proved to belong to a unique parametric family of counterexamples by Bekir .

Note that the equations defining the elements of sets $\mathcal{X}$ and $\mathcal{Y}$ are linear combinations of $\boldsymbol{p}$. This suggests that we can geometrically characterize (with linear subspaces) the supports of signals that generate a turnpike problem without a unique solution.

The sets of points that generate a turnpike problem without a unique solution belong to the following 2–dimensional linear subspaces,

Even if these two linear subspaces define completely the sets $\cal{X}$ and $\cal{Y}$, we need to define a proper ordering of the elements. In fact, while the sets do not have by definition a defined order, the two linear subspaces induce a precise ordering. We can choose any permutation as soon as it is univocally defined: in what follows, we always consider the permutations $\boldsymbol{\Pi}_{{\mathcal{Y}},\boldsymbol{p}}$ and $\boldsymbol{\Pi}_{{\mathcal{X}},\boldsymbol{p}}$ that sort in an increasing order the vectors $\boldsymbol{\Pi}_{{\mathcal{Y}},\boldsymbol{p}}\boldsymbol{Q}_{\mathcal{Y}}\boldsymbol{p}$ and $\boldsymbol{\Pi}_{{\mathcal{X}},\boldsymbol{p}}\boldsymbol{Q}_{\mathcal{X}}\boldsymbol{p}$, respectively. Equivalently, we always consider an operator that takes the elements from the sets $\cal{X}$ and $\cal{Y}$ and sorts them in an increasing order.

Note that the permutation and the operator are unique for each $\boldsymbol{p}$ and depend only on the direction of $\boldsymbol{p}$: in fact, for different magnitudes of $\boldsymbol{p}$, the ordering does not change. Moreover, it is easy to show that we have a finite number of permutations for all the $\boldsymbol{p}$. In fact, the number of permutations is upper bounded by 6 factorial, being the number of possible permutation of 6 elements in a set.

We defined the geometry of these sets of points as a manifold generated by a linear model and a varying permutation. In what follows, we take the span of these permuted linear systems to describe the sets of supports without a unique solution to the turnpike problem.

The set of supports for a signal $f(x)$ that generate a turnpike problem without a unique solution is described as,

A graphical representation of this set of elements is given in Figure 3, where we observe the subset of linear subspaces mixed by the permutations.

This geometrical intuition is useful for two reasons:

if we have no collisions and unless $N=6$, we can always recover uniquely the support of the signal from the support of the ACF,

even if $N=6$, the supports without a unique recovery lie on the 2–dimensional manifold defined in Corollary 2, and this manifold has measure zero in the set of all the supports of 6 elements.

Note that we have not used the coefficient information so far. The following theorem merges the previous results in terms of phase retrieval for sparse signals and considers the coefficients of the ACF to obtain a sufficient condition for the uniqueness of 1–dimensional sparse PR problems.

Assume we measure the 1–dimensional ACF of a signal with $N$ deltas and the elements of the ACF have no collisions. Then,

If $N\neq 6$, the PR problem has a unique solution.

If $N=6$ and not all the $c^{(n)}$ have the same value, the PR problem has a unique solution.

If $N=6$ and all the $c^{(n)}$ have the same value, the PR problem has almost surely a unique solution.

The proof is given in Appendix -A. The three cases are due to the parametric family of turnpike problems without a unique solution previously described and the additional information that may be available from the coefficients $c^{(n)}$ of the deltas of the ACF.

IV-B Uniqueness condition: collision-free D𝐷D–dimensional ACFs

The previous analysis applies only to 1–dimensional signals, that is $D=1$. Senechal et al. proposed an analysis of the uniqueness of the turnpike problem in higher dimensions, $D\geq 2$. Unfortunately, their result is too conservative and cannot cover very simple examples such as 1D sets of points embedded in higher dimensional spaces.

In what follows, we describe the result of Senechal et al. and propose a sufficient condition for the uniqueness of the PR for multi–dimensional sparse signals.

A point belonging to a multi–dimensional difference set, or equivalently a delta of the ACF, is visible if the line through the origin and the point contains only the origin, the point itself and another point in the centro–symmetric position w.r.t. the origin.

Note that there is always a point in a centro–symmetric position, see (4). In Figure 4, we show an ACF with some visible deltas, such as $\boldsymbol{y}^{(j)}$, and some deltas that are not visible, such as $\boldsymbol{y}^{(n)}$.

A signal $f(\boldsymbol{x})$ composed of $N$ deltas has its elements in general position if every delta of the ACF is visible.

Senechal et al. showed that having the points in general position is a sufficient condition to recover uniquely the support of a sparse signal $f(\boldsymbol{x})$ from the support of its ACF.

For example, Theorem 3 does not guarantee a unique recovery of the signal generating the ACF given in Figure 4. Note that while the proposed uniqueness condition for the one–dimensional PR problem described in Theorem 2 and the one for multi–dimensional ones given in Theorem 3 are both sufficient, the latter is fairly constraining. For example, if we pick a 1–dimensional sets of points that generates a turnpike problem with a unique solution and embed it into a higher dimensional domain, then Theorem 3 cannot guarantee anymore the uniqueness of the solution.

This example inspired the idea of solving a $D$–dimensional PR as a set of multiple $1$–dimensional PR problems. In fact, if we can solve uniquely the sub–problems, then the original problem may have a unique solution. In this section, we show that this divide and conquer strategy leads to a tighter necessary and sufficient condition: given an ACF $a(\boldsymbol{x})$ of a $N$–sparse $D$–dimensional signal $f(\boldsymbol{x})$, the PR problem has a unique solution.

We show the result as follows: first we consider a set of projections of the ACF to many 1–dimensional subspacesIt is possible to consider projections onto subspaces of any dimensionality. Here, we consider only 1–dimensional projections for simplicity of notation.. See Figure 5 for an example of a projection over a subspace defined by a vector $\boldsymbol{P}$. Second, we show that the projection of the ACF is the ACF of the projected signal. Third, we show that a finite number of different projections is necessary and sufficient to recover the deltas of $f(\boldsymbol{x})$. We conclude showing that we can find these projections for every $N$-sparse $f(\boldsymbol{x})$ that is embedded in a $D$–dimensional space, with $D\geq 2$.

Given Proposition 1, we can project a $D$–dimensional ACF $a(\boldsymbol{x})$ onto $P$ different 1–dimensional subspaces. If the projected ACF satisfies the conditions given in Theorem 2, we recover the projected signal, namely the deltas with the proper coefficient and the projected locations. As shown in , we need $(N+1)^{(D-1)}$ projections to reconstruct exactly the location of $N$ deltas in a $D$–dimensional space. It is possible to reduce the number of required projection to $D+1$ accepting to take random projections and exactly recovering the support with probability one.

Finally, we show the existence of the projections with a unique PR for any $D$–dimensional ACF with $D\geq 2$. More precisely, a projection may unfortunately belong to the parametric family described in Theorem 1 or it can have collisions in the projected points. In what follows, we show that the projections without collisions and with a unique PR exist and are easy to find for all the $D$–dimensional signals with $D\geq 2$. Indeed, it is simple to show that for any ACF without collisions, the projected ACF on a random subspace has no collisions with probability one. In the following theorem, we show that the set of projections with a unique PR does not have measure zero in the set of all projections for any ACF with $D\geq 2$. In other words, the multi–dimensional sparse PR problem always has a unique solution.

Let $f(\boldsymbol{x})$ be a $D$–dimensional $N$–sparse signal and let $a(\boldsymbol{x})$ be its ACF. Assume that $a(\boldsymbol{x})$ has no collisions and $D\geq 2$. Then the set of projections onto 1–dimensional domains that generates a 1–dimensional PR problem with a unique solution has a measure larger than zero in the set of all the possible projections. Therefore, the solution to the PR of $f(x)$ is always unique.

We conclude this section underlying that the uniqueness of the solution of the sparse PR problem does not depend on the type of domain—that is continuous or discrete. The presence of collisions is a better characterizing feature. However, the amount of collisions depends on the domain: if we randomly distribute the deltas of $f(\boldsymbol{x})$ on a continuous domain, we have almost surely no collision, while the probability of collisions for discrete signals is always larger than zero for $N\geq 4$.

V Conclusions

As we described in Section II-C, many real-world problems require the recovery a sparse signal on a continuous domain from the magnitude of its Fourier transform. In this work, we provided an answer to the uniqueness of the solution of such problem, whether defined on a continuous or a discrete domain.

In particular, we showed that the uniqueness of the solution of such problem depends on the presence of collisions. More precisely, the solution is always unique for signals that are embedded on 2 or more dimensions and whose ACFs do not have collisions. For 1–dimensional signals, we showed the uniqueness for most of the signals with a collision-free ACF, the only counterexample being signals composed of $N=6$ elements with the same coefficient supported on a set of points belonging to the unique family of counterexamples .

Future work will be focused on developing algorithms to solve the 1–dimensional PR problem on the continuous domain and its extension to the multi–dimensional setup using multiple projections on different 1–dimensional domains.

We divide the proof in three parts, one for each case given in the theorem’s statement.

$\bullet$ For $N\neq 6$ we can always recover uniquely the locations of the deltas from the sets of differences, see Theorem 1. Once the locations are known, the coefficients are uniquely determined, see Appendix -D.

$\bullet$ If $N=6$ and all the $c^{(n)}$ have the same coefficient, and the signal does not have a unique PR, then its support lies on the manifold defined in Corollary 2 and we cannot use the coefficients to enforce uniqueness. Note that all possible signals just $N=6$ span a 6–dimensional space, representing the six locations. Given that the manifold containing the signals without a unique reconstruction is 2–dimensional, then the set of these signals has measure zero w.r.t. the set of all the signals with $N=6$.

$\bullet$ If $N=6$ and not all the $c^{(n)}$ have the same coefficient, there is always only one set between $\mathcal{X}$ and $\mathcal{Y}$ that is a possible support of the signal $f(x)$. To prove this statement, we assume without loss of generalityIf the deltas of the $f(x)$ are not positive, we can always take the absolute value given the absence of collisions by assumption. that all $d^{(n)}$ are positive and we would like to show that it is always possible to discern the support of $f(x)$ between the two possibile ones using the coefficients.

First, define two vectors $\boldsymbol{q}$ and $\boldsymbol{r}$ containing the logarithms of the coefficients of the ACF and of the signal, respectively. Then, given the absence of collision and the linearization introduced by the logarithm we can write two systems of equations. Each equation represents the coefficient of one element of the ACF given the coefficients of two elements of the original signal, see Figure 6. The two systems assume that $f(x)$ is supported on $\cal{X}$ and $\cal{Y}$, respectively. More precisely, if we assume that $f(x)$ is supported on $\mathcal{X}$ we have

while if we consider $\mathcal{Y}$ to be the support of $f(x)$ we get a different system of equations,

Successively, we note that if all the coefficients have the same value, that is $c^{(n)}=C\;\forall\,n=\{1,\ldots,N\}$, all the ACF elements have the same coefficient as well, $d^{(n)}=C^{2}$. In this case, (6) and (7) are equivalent and we cannot distinguish between the two possible supports, $\mathcal{X}$ and $\mathcal{Y}$. Note that the set of signals having $c^{(n)}=C\;\forall n$ form a 1–dimensional subspace in the coefficient domain. In what follows, we show that this is the only subspace where $\mathcal{X}$ and $\mathcal{Y}$ cannot be distinguished from the coefficients.

Consider the intersection between the two columns spaces of $\boldsymbol{C}_{\mathcal{X}}$ and $\boldsymbol{C}_{\mathcal{Y}}$. This intersection contains all the coefficients of the ACF that could be equivalently generated by two signals, one supported on $\mathcal{X}$ and one on $\mathcal{Y}$. The dimensionality of this intersection can be computed as

where $[\boldsymbol{C}_{\mathcal{X}},\boldsymbol{C}_{\mathcal{Y}}]$ is a symbol representing the horizontal concatenation of the two matrices $\boldsymbol{C}_{\mathcal{X}}$ and $\boldsymbol{C}_{\mathcal{Y}}$. It is possible to verify that (8) holds for the example of Bloom and by linearity to any other element of the set of counterexamples.

We deduce from (8) that the ACFs with coefficients $\left\{d^{(n)}\right\}_{n=0}^{N^{2}-N}$ that can be generated by the two different supports lie on a 1–dimensional subspace. Given that we have already characterized this 1–dimensional subspace as the one where all the $\left\{c^{(n)}\right\}_{n=0}^{5}$ have the same value, this concludes the proof. ∎

-B Proof of Proposition 1

First, we define the projection $\bar{f}(s)$ of the signal $f(\boldsymbol{x})$ over the subspace indicated by $\boldsymbol{P}$ as the following inner product

Then, we compute the ACF of the projected signal as

where the $\delta(s)$ are located on the 1–dimensional domain defined by the projection. Finally, we compute the projection of the ACF according to (3) as

that is equal to (9), which proves the proposition. ∎

-C Proof of Theorem 4

First, we introduce some notation. The projected deltas are located on $\{w^{(n)}\}_{n=0}^{N^{2}-N}$. These locations can be computed with the following matrix vector multiplications,

Similarly to the analysis proposed in Section IV-A, we can geometrically characterize the supports of the ACFs for which we cannot solve the turnpike problem uniquely. We build a linear model $\boldsymbol{Q}$ from the difference set induced by $\cal{X}$, or $\cal{Y}$ equivalently. Then, we define a permutation matrix $\Pi_{\boldsymbol{p}}$ to match the ordering given by the linear model with the ordering chosen for the set of points representing the ACF support. Assume that the projected points $\left\{w^{(n)}\right\}_{n=0}^{N^{2}-N}$ are supported on a set of points without a unique PR. Then, we have

where $\boldsymbol{Q}$ spans a 2-dimensional linear subspace and when combined with the $\boldsymbol{p}$-dependent permutation $\Pi_{\boldsymbol{p}}$, it defines a 2-dimensional manifold.

if the ACF is 2–dimensional and the intersection is 1–dimensional, there is only one projection that does not have a unique reconstruction and we can pick any other $\boldsymbol{P}$ to have a unique PR. On the other hand, if the intersection is 2–dimensional, then the signal subspace is the same as the one spanned by the matrix $\boldsymbol{Q}$. However, we can change $\boldsymbol{p}$ such that the permutation $\boldsymbol{\Pi}_{\boldsymbol{p}}$ changes and the projected points have a unique 1–dimensional PR. Given the structure of the manifold induced by the permutations, the good projections are relatively easy to find.

if the ACF is $D$–dimensional and $D>2$, the set of projections $\boldsymbol{P}$ that generates projected PR problems with a unique solutions is dense in the set of the projections. In fact, the set of projections without a unique solutions is a 2–dimensional manifold while all the projections are a $D$–dimensional set. We always find sufficient projections to recover the support of $f(\boldsymbol{x})$ and therefore all the $f(\boldsymbol{x})$ without collisions have a unique PR.

-D Recovering the coefficients from the support

In general, the problem of recovering the coefficients $\{c^{(n)}\}_{n=1}^{N}$ of a $N$ sparse signal $f(x)$ knowing the support $\{x^{(n)}\}_{n=1}^{N}$ and the ACF $a(x)$ is not trivial. It is equivalent to solve a system of quadratic forms and a possible convex relaxation is given in . Unfortunately, it is not possible to guarantee the success of reconstruction.

The problem has the same formulation whether the domain is discrete or continuous. In what follows we derive an algorithm that recovers the coefficients of the a sparse signal whose support and ACF are both known while requiring the absence of collision in its ACF.

Let $\mathbf{c}=[c^{(1)},c^{(2)},\ldots,c^{(N)}]^{\top}$ be a vector containing the coefficients of a $N$-sparse signal $f(\boldsymbol{x})$ and define a rank-one matrix $\mathbf{C}=\mathbf{cc}^{*}$

Given that we know the ACF and the support of $f(\boldsymbol{x})$, then we know all the off-diagonal elements of the matrix $\mathbf{C}$. Therefore, we can reformulate the task of recovering the coefficients from the support as follows.

Consider a rank-one matrix $\mathbf{C}$ whose elements are all known except the ones on the main diagonal. Can we uniquely reconstruct the elements on the main diagonal?

There exists many alternative approaches to solve Problem 6. In what follows, we describe a method that requires a single matrix inversion.

Let $\alpha_{k}\triangleq|c_{k}|^{2}$ $\forall k=1,\ldots,N$, and define a matrix $\mathbf{A}=\mathbf{C}-\operatorname{diag}(\alpha_{1},\ldots,\alpha_{N})$. Note that $\mathbf{A}$ is known and we aim to recover the matrix $\mathbf{C}$. The following result describes a method that requires a single matrix inversion to find the unique solution of Problem 6.

where $\mathbf{C}_{k,k}$ and ${\mathbf{A}^{-1}}_{k,k}$ are the $k$-th element on the main diagonal of $\mathbf{C}$ and $\mathbf{A}^{-1}$, respectively.

Denote by $\mathbf{D}=\operatorname{diag}(\alpha_{1},\ldots,\alpha_{N})$, then we have $\mathbf{A}=-\mathbf{D}+\mathbf{cc}^{*}$. Applying the matrix inversion Lemma, we obtain

We conclude the proof noticing that the $k$-th element on the main diagonal of $\mathbf{A}^{-1}$ is a function of $N$ and $\alpha_{k}$,

Once we have recovered $\mathbf{C}$, we obtain the coefficients $\mathbf{c}$ by taking the eigenvector of $\mathbf{C}$ corresponding to the largest eigenvalue.

Acknowledgment

The authors would like to thank Prof. Gervais Chapuis for the insightful discussions on X-ray crystallography. This research was supported by an ERC Advanced Grant – Support for Frontier Research – SPARSAM Nr: 247006.

References