# Signal Recovery From Product of Two Vandermonde Matrices

Dževdan Kapetanović
Huawei Technologies Sweden AB
Lund Research Center, Sweden

###### Abstract

In this work, we present some new results for compressed sensing and phase
retrieval. For compressed sensing, it is shown that if the unknown
$n$-dimensional vector can be expressed as a linear combination of $s$ unknown
Vandermonde vectors (with Fourier vectors as a special case) and the
measurement matrix is a Vandermonde matrix, exact recovery of the vector with
$2s$ measurements and $O(\mathrm{poly}(s))$ complexity is possible when $n\geq
2s$. Hence, in this case, there is no need for random measurements (with the
number of measurements scaled by $\log(n)$) as required by conventional
compressed sensing theory. From these results, a measurement matrix is
constructed from which it is possible to recover $s$-sparse $n$-dimensional
vectors for $n\geq 2s$ with as few as $2s$ measurements and with a recovery
algorithm of $O(\mathrm{poly}(s))$ complexity. This seems to be the first
measurement matrix construction achieving exact recovery of $s$-sparse vectors
with $O(s)$ measurements (achieving the lower bound of $2s$) and
$O(\mathrm{poly}(s))$ complexity (typically, the recovery complexity is
$O(\mathrm{poly}(n))$), significantly improving upon a classical construction
by Bourgain et al. In the second part of the work, these results are extended
to the challenging problem of phase retrieval (compressed sensing with
magnitude-only observations). The most significant discovery in this direction
is that if the unknown $n$-dimensional vector is composed of $s$ frequencies
with at least one being non-harmonic, $n\geq 4s-1$ and we take at least $8s-3$
Fourier measurements, there are, remarkably, only two possible vectors
producing the observed measurement values and they are easily obtainable from
each other. The two vectors can be found by an algorithm with only
$O(\mathrm{poly}(s))$ complexity. This result is interesting when put in
relation to recent work by Xu which shows that solving Fourier phase retrieval
with $O(\mathrm{poly}(n))$ complexity proves that P = NP holds for the well
known P vs NP problem. An immediate application of the new result is
construction of a measurement matrix from which it is possible to recover all
$s$-sparse $n$-dimensional signals (up to a global phase) from $O(s)$
magnitude-only measurements and $O(\mathrm{poly}(s))$ recovery complexity when
$n\geq 4s-1$. As a consequence, the proposed matrix construction solves (and
significantly improves upon) an open problem by Xu on finding a measurement
matrix construction from which one can (provably) recover any $s$-sparse
vector with $O(s\log(n/s))$ phase-less measurements and 1-norm minimization.

###### Index Terms:

Compressed sensing, sparse signal estimation, phase retrieval, sparse phase
retrieval.

##  I Introduction

We start by considering the conventional compressed sensing problem, which
amounts to solving for $\bm{x}$ the following linear system

$\displaystyle\bm{y}=\bm{A}\bm{x},\quad\bm{x}\in\mathbb{C}_{\leq s}^{n},$ (1)

where $\mathbb{C}^{n}_{\leq s}$ is the space of $n$-dimensional complex
vectors with at most $s$ elements being non-zero and $\bm{A}$ is a known
$m\times n$ measurement matrix ($m$ is the number of measurements). This
problem has a rich history. Seminal works [1, 2, 3] and references therein
show that, remarkably, there exist matrices $\bm{A}$ such that $\bm{x}$ can be
recovered with far less measurements than $n$. We will refer to exact recovery
of $\bm{x}$ as recovery with probability 1 (abbreviated as w.p.1) while
probabilistic recovery refers to recovery at some probability strictly smaller
than 1. The obvious 0-norm minimization that searches across all subsets of
$s$ columns from $\bm{A}$ and inverts them to obtain a candidate solution is
of exponential complexity. Fortunately, the work in [4] shows that there
exists a recovery algorithm for (1) with $O(\mathrm{poly}(n))$ complexity,
where $O(.)$ stands for big O notation. More specifically, [4] shows that, for
a given $\bm{y}=\bm{A}\bm{x}$, the unique solution to the convex 1-norm
optimization problem

$\displaystyle\min_{\bm{z}}\|\bm{z}\|_{1}$ subject to
$\displaystyle\bm{y}=\bm{A}\bm{z}$ (2)

is $\bm{z}=\bm{x}$ if $\bm{A}$ satisfies the restricted isometry property
(RIP). Some random constructions of $\bm{A}$ with $m=O(s\log(n/s))$ number of
rows satisfy the RIP property with high probability. These existence and low
recovery complexity results is what made compressed sensing a breakthrough
technique with vast practical applications [5]. Since random constructions of
$\bm{A}$ only result in probabilistic recovery of $\bm{x}$, constructions that
guarantee exact recovery are also of high interest. In this direction, the
work in [7] shows that for large $n$ there are deterministic constructions of
$\bm{A}$ with $m=o(s^{2})$ rows ($o(.)$ stands for Little o notation) and
$n\geq m$ satisfying the RIP property and thus solving (2). To the best of my
knowledge, the quadratic guarantee of $m=o(s^{2})$ is the lowest to date.

The problem in (1) assumes full knowledge of the complex valued vector
$\bm{y}$. A related problem is the recovery of $\bm{x}\in\mathbb{C}^{n}_{s}$
from phase-less measurements

$\displaystyle\bm{y}=|\bm{A}\bm{x}|^{2},\quad\bm{x}\in\mathbb{C}_{s}^{n},$ (3)

where $|\bm{x}|^{2}$ is the vector of squared magnitudes of the elements in
$\bm{x}$. The problem in (3) is known as sparse phase retrieval in the
literature [8]-[10]. The recovery is allowed to be correct up to a common
phase, i.e., obtaining $e^{i\alpha}\bm{x}$ for any $\alpha$ is considered as
exact recovery. Clearly, due to the absence of phase information, this problem
is significantly more challenging than the compressed sensing problem in (1)
(which measurements we refer to as phase-aware measurements); as such, more
measurements are expected for exact recovery. Injectivity studies have shown
that if $\bm{A}$ is a generic matrix [15] then $m=4s-4$ measurements suffice
for guaranteeing a unique solution. It is conjectured that the necessary lower
bound is of the form $4s-O(s)$. The injectivity studies do not propose an
efficient algorithm to check whether a certain $\bm{A}$ is generic nor do they
provide an efficient recovery algorithm. A reformulation of (2) to (3) results
in the recovery algorithm

$\displaystyle\min_{\bm{z}}\|\bm{z}\|_{1}$ subject to $\displaystyle\bm{y}$
$\displaystyle=|\bm{A}\bm{x}|^{2}=|\bm{A}\bm{z}|^{2}.$ (4)

Unfortunately, since (4) is a highly non-convex problem, it is no longer
possible to prove that the unique solution to (4) is
$\bm{z}=e^{i\alpha}\bm{x}$ for some $\alpha$. The work in [16] shows that if
$\bm{A}$ is a circularly symmetric complex Gaussian matrix and
$m=O(s\log(n/s))$, then solving (4) results in $\bm{z}=e^{i\alpha}\bm{x}$ with
probability at least $1-2e^{-c_{0}m}$ where $c_{0}$ is a constant; hence, only
probabilistic recovery is guaranteed. When $\bm{A}$ is a discrete Fourier
matrix, an iterative projection scheme is presented in [11] that solves (4)
with high probability (i.e., exact recovery cannot be guaranteed) when
$m=O(s^{2}\log(n/s^{2}))$. However, the scheme is highly complex since it uses
many iterations (200 000) for guaranteeing satisfactory performance and it
assumes knowledge of $\|\bm{x}\|_{1}$ which is an impractical assumption. The
complexity of solving (4) for a general $\bm{A}$ is at least
$O(\mathrm{poly}(n))$ (probably significantly higher and possibly
exponential), since the problem is at least as difficult as (2).

For general $\bm{A}$, PhaseLiftOff [12] is a state-of-art algorithm that
provides high probabilistic recovery with $O(s\log(n/s))$ measurements when
elements in $\bm{A}$ are chosen from a complex Gaussian distribution. In order
to avoid random constructions, PhaseCode [13] utilizes the theory of sparse-
graph codes to provide a deterministic construction of $\bm{A}$ that is able
to recover with high probability a large fraction of elements in $\bm{x}$ with
$O(s)$ measurements. As a particular operating point, PhaseCode is able to
recover a fraction of $1-10^{-7}$ non-zero components of $\bm{x}$ with $14s$
measurements and recovery probability $1-O(1/m)$. Requiring a fraction of
recovered elements close to 1 and a recovery probability close to 1 (i.e.,
almost exact recovery) drives the number of measurements to infinity; hence,
as for PhaseLiftOff, exact recovery eludes PhaseCode as well.

In many practical instances of compressed sensing and phase retrieval,
$\bm{A}$ is constrained to be a Fourier matrix [1, 17, 11] which in practice
corresponds to an apparatus that takes frequency measurements of an underlying
signal. Compressed sensing and (sparse) phase retrieval with Fourier
measurement matrices, known as Fourier compressed sensing and (sparse) Fourier
phase retrieval, respectively, are notoriously difficult problems, as
described in [1, 6, 17]. The main challenge is in obtaining an affordable
recovery complexity (e.g. $O(\mathrm{poly}(n))$) with few measurements (e.g.
$O(s\log(n/s))$). The breakthrough work in compressed sensing [1] considered
random Fourier measurements in order to escape the deterministic nature of the
problem and instead provide probabilistic recovery with complexity
$O(\mathrm{poly}(n))$ and $O(s\log(n/s))$ measurements. The contributions in
this work on compressed sensing show that if $\bm{x}$ is instead $s$-sparse in
a Vandermonde basis, with $n\geq 2s$, then it can be recovered exactly from
$O(s)$ deterministic Fourier measurements and $O(\mathrm{poly}(s))$ recovery
complexity; hence, there is no need for random measurements in this case. In
case of (non-sparse) Fourier phase retrieval, as shown in [6], the reason for
why no exact $O(\mathrm{poly}(n))$ complexity algorithms have been proposed so
far is because in doing so one solves the well known P vs NP problem, proving
that P = NP holds. Interestingly, the contributions in this work on Fourier
phase retrieval show that if $\bm{x}$ is $s$-sparse in a non-DFT basis, with
$n\geq 4s-1$, then it can be recovered from $O(s)$ measurements and
$O(\mathrm{poly}(s))$ recovery complexity.

###  I-A Contributions

Let $\bm{V}(\bm{z})$ denote a Vandermonde matrix of $n$ rows parameterized by
the elements from vector $\bm{z}=[z_{0},\ldots,z_{m-1}]$

$\displaystyle\bm{V}(\bm{z})=\left[\begin{array}[]{ccc}1&\ldots&1\\\
z_{0}&\ldots&z_{m-1}\\\ \vdots&\ddots&\vdots\\\
z_{0}^{n-1}&\ldots&z_{m-1}^{n-1}\end{array}\right].$ (9)

Let $\mathbb{T}=\\{z;z\in\mathbb{C},|z|=1\\}$ be the set of all points on the
complex unit circle and $(.)^{T}$ the matrix transpose operator. When
$z_{k}\in\mathbb{T}$, $0\leq k\leq m-1$, the Vandermonde matrix becomes a
Fourier matrix, which we for clarity henceforth denote as $\bm{F}(\bm{z})$
(with DFT matrix as a special case when $z_{k}$ are $n$th roots of unity and
$m=n$). For Vandermonde and Fourier measurement matrices, we have the
following main contributions

  1. 1.

When $\bm{A}=\bm{V}(\bm{z})^{T}$, $\bm{x}$ is of the form
$\bm{x}=\bm{V}(\bm{\theta})\bm{g}$ for unknown
$\bm{\theta}\in\mathbb{C}^{s},\bm{g}\in\mathbb{C}^{s}$, and $n\geq 2s$, we
have exact recovery of $\bm{\theta}$ and $\bm{g}$ from $\bm{y}=\bm{A}\bm{x}$
with $O(\mathrm{poly}(s))$ recovery complexity and

    1. a.

When the elements in $\bm{z}$ are $n$th roots of unity, $m=2s$ measurements
suffice.

    2. b.

For arbitrary $\bm{z}\in\mathbb{C}^{m}$, $m\geq 3s$ suffices.

  2. 2.

Assume $\bm{A}=\bm{F}(\bm{z})^{T}$, $\bm{x}$ is of the form
$\bm{x}=\bm{F}(\bm{\theta})\bm{g}$, with $\bm{\theta}\in\mathbb{T}^{s}$ and
$\bm{g}\in\mathbb{C}^{s}$ being unknown, and $n\geq 4s-1$. Let
$\bm{y}=|\bm{A}\bm{x}|^{2}$.

    1. a.

If the elements in $\bm{\theta}$ are $n$th roots of unity
($\bm{F}(\bm{\theta})$ equals a subset of $s$ columns from a $n\times n$ DFT
matrix), we have exact recovery of $\bm{\theta}$ and $|\bm{g}|$ (up to a
global real scalar) from $\bm{y}$ and $\bm{A}$ with $O(\mathrm{poly}(s))$
recovery complexity. There are $2^{s-1}$ different $\bm{g}$ (up to a global
phase), all having the same value of $|\bm{g}|$, producing the same $\bm{y}$.

    2. b.

If not all elements in $\bm{\theta}$ are $n$th roots of unity then

      1. i.

When the elements in $\bm{z}$ are $n$th roots of unity and $m=4s-1$, we have
exact recovery of $\bm{\theta}$ and $|\bm{g}|$ (up to a global real scalar)
from $\bm{y}$ and $\bm{A}$ with $O(\mathrm{poly}(s))$ recovery complexity.
There are $2^{s-1}$ different $\bm{g}$ (up to a global phase), all having the
same value of $|\bm{g}|$, producing the same $\bm{y}$.

      2. ii.

For arbitrary $\bm{z}\in\mathbb{T}^{m}$ with $m\geq 8s-3$, we have exact
recovery of $\bm{\theta}$ and $|\bm{g}|$ (up to a global real scalar) from
$\bm{y}$ and $\bm{A}$ with $O(\mathrm{poly}(s))$ recovery complexity. There
are only two different $\bm{g}$ (up to a global phase), both having the same
value of $|\bm{g}|$, that produce the same $\bm{y}$; moreover, the two
solutions are easily obtainable from each other.

By an arbitrary vector $\bm{a}\in\mathbb{C}^{k}$, we mean a vector chosen at
random with some probability distribution that is continuous on
$\mathbb{C}^{k}$.

The result in 1. states that we can obtain exact recovery of $n$-dimensional
$\bm{x}$ with $O(s)$ phase-aware Vandermonde measurements and
$O(\mathrm{poly}(s))$ complexity if $\bm{x}$ is a linear combination of $s$
Vandermonde vectors and $n\geq 2s$. Thus, although efficient and exact
recovery of sparse $\bm{x}$ from Fourier measurements is challenging, as
illustrated in [1], 1. shows that it is surprisingly simple when $\bm{x}$ is
sparse in a Vandermonde basis (note that the sparsity level is all the way up
to $n/2$).

For phase-less measurements, the result in 2. states that $\bm{\theta}$ and
$c|\bm{g}|^{2}$, for unknown real scalar $c$, can be found with $O(s)$
measurements and $O(\mathrm{poly}(s))$ recovery complexity. Thus, in scenarios
where only the frequency components $\bm{\theta}$ and their relative
magnitudes are sought for, the needed number of measurements and the recovery
complexity are essentially lowest possible. Moreover, the different vectors
$\bm{g}$ that produce the same measurements $\bm{y}$ all have same magnitudes
of elements on same positions, meaning that only the phases of the elements
differ.

On the other hand, when it comes to recovering $\bm{g}$ from phase-less
measurements, 2.2.a. is in line with the work in [6], which showed that
solving Fourier phase retrieval with $O(\mathrm{poly}(n))$ recovery complexity
for general $\bm{x}$ proves that P = NP holds for the well known P vs NP
problem. As seen, 2.2.a. states that recovering an oversampled signal composed
of $s$ harmonic frequencies with an oversampling factor of at least 4 ($n\geq
4s-1$) still results in an exponential recovery complexity in $s$ since there
are $2^{s-1}$ different $\bm{g}$ giving rise to the measurements in $\bm{y}$.
Therefore, this work does not prove that P = NP, since it doesn’t provide
polynomial recovery complexity for all $\bm{x}$. Interestingly, as noted
above, the $2^{s-1}$ different vectors $\bm{g}$ only differ in the phases of
their elements, as if there are only two choices for the phase of each element
(once the phase of one element is fixed).

Conversely, 2.2.b. states that when $\bm{x}$ is composed of $s$ frequencies
with at least one being non-harmonic, where $n\geq 4s-1$ and the number of
measurements $m\geq 8s-3$, there are only two possible solutions (one of them
is of course the true vector $\bm{x}$) and they can be recovered with
$O(\mathrm{poly}(s))$ complexity. Thus, Fourier phase retrieval is in
principle solvable with $O(s)$ measurements and $O(\mathrm{poly}(s))$ recovery
complexity for signals containing non-harmonic frequencies and oversampled
with at least an oversampling factor of 4. Obviously, to find which of the two
possible solutions is $\bm{x}$, one can take an additional measurement with an
arbitrary measurement vector for resolving the ambiguity; this is used next to
construct $\bm{A}$ for recovering $s$-sparse $\bm{x}$ from phase-less
measurements.

Removing the assumption of $\bm{x}$ being spanned by a Vandermonde or Fourier
basis, we look at a measurement matrix construction $\bm{A}$ for recovering
$\bm{x}\in\mathbb{C}_{\leq s}^{n}$. Let ; denote the vertical concatenation of
matrices. In this regard, the main contributions of this work can be
summarized as follows

  1. 3.

When $\bm{A}=\bm{V}(\bm{z})^{T}\bm{V}(\bm{\theta})$,
$\bm{z}\in\mathbb{C}^{m}$, $\bm{\theta}\in\mathbb{C}^{n}$, and $n\geq 2s$, we
get exact recovery of $\bm{x}\in\mathbb{C}_{\leq s}^{n}$ from
$\bm{y}=\bm{A}\bm{x}$ with $O(\mathrm{poly}(s))$ recovery complexity and

    1. a.

When elements in $\bm{z}$ are chosen as $n$th roots of unity and for almost
all $\bm{\theta}$, $m=2s$ measurements suffice.

    2. b.

For almost all $\bm{z}$ and $\bm{\theta}$, $m\geq 3s$ suffices.

  1. 4.

Let $\bm{y}=|\bm{A}\bm{x}|^{2}$, where $\bm{x}\in\mathbb{C}_{\leq s}^{n}$,
$\bm{z}\in\mathbb{T}^{m}$, $\bm{\theta}\in\mathbb{T}^{n}$ and $n\geq 4s-1$.

    1. a.

Assume $\bm{A}=\bm{F}(\bm{z})^{T}\bm{F}(\bm{\theta})$.

      1. i.

When $m=4s-1$, the elements in $\bm{z}$ are $n$th roots of unity and for
almost all $\bm{\theta}$, we have exact recovery of the support of
$\bm{x}\in\mathbb{C}_{\leq s}^{n}$ and $|\bm{x}|$ (up to a global real scalar)
from $\bm{y}$ and $\bm{A}$. There are $2^{s-1}$ different $\bm{x}$ (all with
the same value of $|\bm{x}|$) producing the same $\bm{y}$.

      2. ii.

For almost all $\bm{z}$ and $\bm{\theta}$ with $m\geq 8s-3$, we have exact
recovery of the support of $\bm{x}\in\mathbb{C}_{\leq s}^{n}$ and $|\bm{x}|$
(up to a global real scalar) from $\bm{y}$ and $\bm{A}$. There are only two
different $\bm{x}$ (with the same value of $|\bm{x}|$) producing the same
$\bm{y}$ and they are easily obtainable from each other.

    2. b.

Let $\bm{A}=[\bm{F}(\bm{z})^{T}\bm{F}(\bm{\theta});\bm{a}^{T}]$ where $\bm{a}$
is an arbitrary vector. Then,

      1. i.

When the elements in $\bm{z}$ are $n$th roots of unity and for almost all
$\bm{\theta}$, $m=4s-1$ measurements suffice to recover $\bm{x}$ from $\bm{y}$
and $\bm{A}$ with a recovery complexity of $O(\mathrm{poly}(s)2^{s})$.

      2. ii.

For almost all $\bm{z}$ and $\bm{\theta}$, $m\geq 8s-3$ suffices to recover
$\bm{x}$ from $\bm{y}$ and $\bm{A}$ with $O(\mathrm{poly}(s))$ recovery
complexity.

Almost all vectors $\bm{a}$ from a set $\mathbb{A}$ equals the set
$\mathbb{A}$ up to a null set.

The result in 3. significantly improves upon the result in [7]: we provide an
explicit construction of $\bm{A}$ with $m=O(s)$ rows (even achieving the lower
bound $2s$) that solves (2). Moreover, for this specific construction of
$\bm{A}$, we present an algorithm that recovers $\bm{x}$ with
$O(\mathrm{poly}(s))$ recovery complexity instead of $O(\mathrm{poly}(n))$.
The reason for our significant improvement in recovery complexity is that we
do not use compressed sensing theory for recovery, but instead develop a
novel, tailor-suited algorithm for the proposed matrix construction.

The results in 4. are obtained rather easily from 2. which can be seen from
the similarity of the statements. The result in 4.4.b. directly solves open
problem 3.5 in [18] and significantly improves upon it. Namely, that problem
asks for an explicit construction of $\bm{A}$ with $m=O(s\log(n/s))$ such that
solving (4) results in exact recovery of $\bm{x}$. The combination of the
deterministic construction $\bm{F}(\bm{z})^{T}\bm{F}(\bm{\theta})$ and a
single measurement as described in 4.4.b. shows that after only $4s$
measurements there is a unique solution to $\bm{x}$ (up to a global phase) in
(3); thus even (4) will find it. Noteworthy is that the number of measurements
$4s$ in the proposed construction almost achieves the injectivity lower bound
of $4s-4$. Furthermore, taking $8s-2$ measurements instead of $4s$, we are
able to reduce the recovery complexity to $O(\mathrm{poly}(s))$ with a novel
recovery algorithm, drastically reducing the recovery complexity for large
scale problems (large $n$) compared to (4).

###  I-B Organization

The rest of the article is structured as follows. In Section II we present two
results described in Section I for phase-aware measurements and their
derivations. Section III presents three results described in Section I for
phase-less measurements and their derivations.

##  II Phase-aware Measurements

We start by formulating the Vandermonde measurement result outlined in Section
I.

R1: Assume an $n$-dimensional vector $\bm{x}=\bm{V}(\bm{\theta})\bm{g}$, where
$\bm{\theta}\in\mathbb{C}^{s}$, $\bm{g}\in\mathbb{C}^{s}$ and $n\geq 2s$.
There are samples $\bm{z}\in\mathbb{C}^{m}$, with $m\geq 2s$, such that
$\bm{\theta}$ and $\bm{g}$ (and thus $\bm{x}$) can be exactly recovered from

$\displaystyle\bm{y}=\bm{V}(\bm{z})^{T}\bm{x}=\bm{V}(\bm{z})^{T}\bm{V}(\bm{\theta})\bm{g}$
(10)

with $O(\mathrm{poly}(s))$ recovery complexity.

Since $\theta_{j}$, $0\leq j\leq s-1$, are assumed to be continuous, R1
achieves grid-less compressed sensing [19] with the assumptions that the
measurement matrix is $\bm{V}(\bm{z})^{T}$ and that the $n$-dimensional
complex vector $\bm{x}$ is inside an $s$-dimensional complex subspace (with
$n\geq 2s$) spanned by a $n\times s$ Vandermonde basis matrix. A version of R1
in the specific case of Fourier matrices is presented in [20]. Hence, R1
extends the results in [20] to general Vandermonde matrices.

R1 is derived by first solving for the unknown vector $\bm{\theta}$. Once
$\bm{\theta}$ is known, it is proved that there is a unique $\bm{g}$
satisfying (10) and it can, e.g., be obtained as
$\bm{g}=\left(\bm{V}(\bm{z})^{T}\bm{V}(\bm{\theta})\right)^{\dagger}\bm{y}$
where $(.)^{\dagger}$ denotes the matrix pseudo-inverse. Deriving R1 in this
particular order gives us

R2: Assume $n\geq 2s$ and $m\geq 2s$. There is an $m\times n$ measurement
matrix construction $\bm{A}$ such that we can exactly recover any
$\bm{x}\in\mathbb{C}_{\leq s}^{n}$ from $\bm{y}=\bm{A}\bm{x}$ with
$O(\mathrm{poly}(s))$ recovery complexity.

R2 is a significant improvement to the classical result in [7] where an
$\bm{A}$ with $m=o(s^{2})$ rows was constructed (for large $n$) that satisfies
RIP and thereby solves (1) with $O(\mathrm{poly}(n))$ complexity.

###  II-A Derivation of R1

As observed in [20], element $y_{j}$ in R1 can be seen as a sample
$y_{j}=f(z_{j})$ of the rational function

$\displaystyle f(z)$
$\displaystyle=\sum_{l=0}^{s-1}g_{l}\sum_{j=0}^{n-1}(z\theta_{l})^{j}=\sum_{l=0}^{s-1}g_{l}\frac{(z\theta_{l})^{n}-1}{z\theta_{l}-1}=\frac{u(z)}{v(z)}$
(11) $\displaystyle u(z)$
$\displaystyle\stackrel{{\scriptstyle\triangle}}{{=}}\sum_{l=0}^{s-1}g_{l}((z\theta_{l})^{n}-1)\prod_{\begin{subarray}{c}i=0\\\
i\neq l\\\ g_{i}\neq 0\end{subarray}}^{s-1}(z\theta_{i}-1)$ (12)
$\displaystyle v(z)$
$\displaystyle\stackrel{{\scriptstyle\triangle}}{{=}}\prod_{\begin{subarray}{c}l=0\\\
g_{l}\neq 0\end{subarray}}^{s-1}(z\theta_{l}-1).$ (13)

The polynomial $u(z)$ can be further expressed as

$\displaystyle u(z)$ $\displaystyle=z^{n}\hat{u}(z)+\tilde{u}(z)$
$\displaystyle\hat{u}(z)$
$\displaystyle\stackrel{{\scriptstyle\triangle}}{{=}}\sum_{l=0}^{s-1}g_{l}\theta_{l}^{n}t_{l}(z)$
$\displaystyle\tilde{u}(z)$
$\displaystyle\stackrel{{\scriptstyle\triangle}}{{=}}-\sum_{l=0}^{s-1}g_{l}t_{l}(z)$
$\displaystyle t_{l}(z)$
$\displaystyle\stackrel{{\scriptstyle\triangle}}{{=}}\prod_{\begin{subarray}{c}i=0\\\
i\neq l\\\ g_{i}\neq 0\end{subarray}}^{s-1}(z\theta_{i}-1).$ (14)

Obviously, if some elements in $\bm{g}$ are zero then $s$ is just reduced to a
smaller value. Thus, without loss of generality (WLOG), we will at first
assume that $\bm{g}$ contains no zeros and later, when necessary, comment on
the impact of an unknown reduced $s$ on some developments in this work.
Furthermore, throughout this work, we also assume that $\bm{\theta}$ contains
no zeros.

Next, some results on the introduced polynomials are given that will be
important for this work.

######  Lemma 1.

Assume that $\bm{g}\in\mathbb{C}^{s}$ has no zeros. The polynomials
$t_{0}(z),\ldots,t_{s-1}(z)$ are linearly independent.

###### Proof.

Since no element in $\bm{g}$ is zero, $t_{j}(z),0\leq j\leq s-1$, are all
distinct and of degree $s-1$. Assume that these polynomials are linearly
dependent, so that we can write $\sum_{j=0}^{s-1}c_{j}t_{j}(z)=0$ for some
coefficients $c_{0},\ldots,c_{s-1}$ not all being 0. Since none of
$t_{j}(z),0\leq j\leq s-1$, are zero, at least two of $c_{j},0\leq j\leq s-1$,
are non-zero. Assume then that $c_{j}\neq 0$ and $c_{k}\neq 0$ for some $j\neq
k$. Hence, we can write
$t_{k}(z)=-\frac{1}{c_{k}}\sum_{\begin{subarray}{c}j=0\\\ j\neq
k\end{subarray}}^{s-1}c_{j}t_{j}(z)$. Since $t_{k}(\theta_{k}^{-1})\neq 0$ and
$t_{j}(\theta_{k}^{-1})=0$ for $j\neq k$, it follows that $0\neq
t_{k}(\theta_{k}^{-1})=-\frac{1}{c_{k}}\sum_{\begin{subarray}{c}j=0\\\ j\neq
k\end{subarray}}^{s-1}c_{j}t_{j}(\theta_{k}^{-1})=0$, a contradiction,
implying that $t_{0}(z),\ldots,t_{s-1}(z)$ are linearly independent. ∎

######  Lemma 2.

For an arbitrary $\bm{g}\in\mathbb{C}^{s}$, the roots of $\tilde{u}(z)$ have
multiplicity one w.p.1.

###### Proof.

An arbitrary $\bm{g}\in\mathbb{C}^{s}$ contains no zeros, hence the
polynomials $t_{0}(z),\ldots,t_{s-1}(z)$ are linearly independent from Lemma
1. Varying $\bm{g}$ across $\mathbb{C}^{s}$ makes $\tilde{u}(z)$ span all
polynomials of degree $s-1$. An arbitrary $\bm{g}$ produces therefore an
arbitrary realization of $\tilde{u}(z)$. Each polynomial of degree $s-1$ with
roots $r_{0},\ldots,r_{s-2}$ (where some roots are equal in case of
multiplicity larger than one) can be identified (up to a complex scalar) with
the vector of its roots $\bm{r}=[r_{0},\ldots,r_{s-2}]\in\mathbb{C}^{s-1}$.
Polynomials with roots that have multiplicity larger than one can therefore be
represented as a union of subspaces in $\mathbb{C}^{s-1}$ of dimensions lower
than $s-1$ and thus, for an arbitrary $\bm{g}$, these polynomials occur w.p.0.
Hence, the roots of $\tilde{u}(z)$ have multiplicity one w.p.1. ∎

######  Lemma 3.

Given a $\bm{\theta}$ such that $\theta_{j}^{n},0\leq j\leq s-1$, are not all
equal. Then there exists polynomials $\tilde{u}(z)$, $\hat{u}(z)$ that do not
have any roots in common.

###### Proof.

From (II-A) it follows that

$\displaystyle
g_{j}=\frac{\tilde{u}(\theta_{j}^{-1})}{t_{j}(\theta_{j}^{-1})},\quad 0\leq
j\leq s-1.$ (15)

Hence, the $g_{j},0\leq j\leq s-1$, (and thus $\tilde{u}(z)$, $\hat{u}(z)$ as
well) are uniquely determined from the $s$ samples
$\tilde{u}(\theta_{j}^{-1}),0\leq j\leq s-1$. The identity in (15) gives

$\displaystyle\tilde{u}(z)=\sum_{j=0}^{s-1}\frac{\tilde{u}(\theta_{j}^{-1})}{t_{j}(\theta_{j}^{-1})}t_{j}(z)$
$\displaystyle\hat{u}(z)=\sum_{j=0}^{s-1}\frac{\tilde{u}(\theta_{j}^{-1})}{t_{j}(\theta_{j}^{-1})}\theta_{j}^{n}t_{j}(z).$
(16)

Choose $\tilde{u}(\theta_{j}^{-1})=1,0\leq j\leq s-1$, which gives
$\tilde{u}(z)=1$. Clearly, $\hat{u}(z)\neq 0$ for this choice since the
assumption on $\bm{\theta}$ is that it has no zeros. Assume then instead that
$\hat{u}(z)$ is a constant polynomial, $\hat{u}(z)=c$, for some complex scalar
$c$. Then, from (II-A)

$\displaystyle\hat{u}(z)-c\tilde{u}(z)=\sum_{j=0}^{s-1}\frac{\theta_{j}^{n}-c}{t_{j}(\theta_{j}^{-1})}t_{j}(z)=0$
(17)

which holds if and only if $\theta_{j}^{n}=c,0\leq j\leq s-1$, contradicting
the assumption on $\bm{\theta}$. Hence, $\hat{u}(z)$ has at least degree one.
∎

Next, we show that for arbitrary $\bm{g}$, the resulting
$\tilde{u}(z),\hat{u}(z)$ share no common roots w.p.1 if and only if
$\theta_{j}^{n},0\leq j\leq s-1$, are not all equal. To do this, we utilize
Lemma 3 and the well known result that two polynomials $p(x),q(x)$ of degrees
$m,n$, respectively, share a common root if and only if their resultant
$\mathrm{res}(p,q)$ is zero, where

$\displaystyle\mathrm{res}(p,q)=\left|\begin{array}[]{cccccccc}p_{m}&p_{m-1}&\ldots&p_{0}&0&\ldots&\ldots&0\\\
0&p_{m}&p_{m-1}&\ldots&p_{0}&0&\ldots&0\\\ \vdots&&&\ddots&&\ddots&&\vdots\\\
0&\ldots&\ldots&0&p_{m}&p_{m-1}&\ldots&p_{0}\\\
q_{n}&q_{n-1}&\ldots&q_{0}&0&\ldots&\ldots&0\\\
0&q_{n}&q_{n-1}&\ldots&q_{0}&0&\ldots&0\\\ \vdots&&&\ddots&&\ddots&&\vdots\\\
0&\ldots&\ldots&0&q_{n}&q_{n-1}&\ldots&q_{0}\end{array}\right|$ (26)

and $|.|$ denotes the determinant. The matrix in (26) has dimension $m+n$,
where coefficients of $p(z)$ appear in $n$ rows and coefficients of $q(z)$
appear in $m$ rows.

######  Theorem 1.

For an arbitrary $\bm{g}\in\mathbb{C}^{s}$, $\tilde{u}(z)$ and $\hat{u}(z)$
share no common roots w.p.1. if and only if $\bm{\theta}$ is such that
$\theta_{j}^{n},0\leq j\leq s-1$, are not all equal.

###### Proof.

Note from (II-A) that if $\theta_{j}^{n},0\leq j\leq s-1$, are all equal then
$\tilde{u}(z)=\hat{u}(z)$ and thus they have the same roots. To finish the
proof, we now show that when not all $\theta_{j}^{n},0\leq j\leq s-1$, are
equal, then $\tilde{u}(z)$ and $\hat{u}(z)$ share no common roots w.p.1. for
an arbitrary $\bm{g}$.

From (26), it follows that the elements of the matrix in
$\mathrm{res}(\tilde{u},\hat{u})$ are a linear combination of the elements in
$\bm{g}$. Thus, $r(\bm{g})=\mathrm{res}(\tilde{u},\hat{u})$ is a multivariate
polynomial in the elements of $\bm{g}$. Denote by $\bm{g}^{1}$ the $\bm{g}$
that produces the polynomials of Lemma 3, i.e., $\tilde{u}(z)=1$ and
$\hat{u}(z)$ is of degree at least one. For these polynomials, the matrix in
(26) is an identity matrix and thus $r(\bm{g}^{1})=1$, showing that
$r(\bm{g})$ is not the zero polynomial. The set
$\mathcal{S}=\\{\bm{g}|\bm{g}\in\mathbb{C}^{s},r(\bm{g})=0\\}$ is either empty
or is a hypersurface of dimension at most $s-1$ in the ambient space
$\mathbb{C}^{s}$. An arbitrary $\bm{g}\in\mathbb{C}^{s}$ belongs to
$\mathcal{S}$ w.p.0, showing that $r(\bm{g})\neq 0$ w.p.1, which proves the
theorem. ∎

The proof for the bounds on $n$ and $m$ in [20] holds here as well, and thus
we have the necessary conditions

######  Theorem 2.

If R1 has a (unique) solution then $n\geq 2s$ and $m\geq 2s$.

Denote by $y_{j}^{+}=f^{+}(z_{j})$, $0\leq j\leq m-1$, the observed
measurements which stem from polynomials
$v^{+}(z),\hat{u}^{+}(z),\tilde{u}^{+}(z)$ that correspond to realizations
$\bm{\theta}=\bm{\theta}^{+},\bm{g}=\bm{g}^{+}$. Furthermore, denote by
$\bm{v},\hat{\bm{u}},\tilde{\bm{u}}$ the coefficient vectors of polynomials
$v(z),\hat{u}(z),\tilde{u}(z)$ with degrees $s,s-1,s-1$, respectively.
Consider the set of linear equations

$\displaystyle y_{j}^{+}v(z_{j})-z_{j}^{n}\hat{u}(z_{j})-\tilde{u}(z_{j})$
$\displaystyle=0,\quad 0\leq j\leq m-1$ (27)

in $\bm{v},\hat{\bm{u}},\tilde{\bm{u}}$. Let $\bm{0}_{n}$ denote the vector of
$n$ zeros (in case its dimension is implicit from the context, we will just
denote it as $\bm{0}$). The equation system in (27) corresponds to solving the
linear equation system

$\displaystyle\bm{A}(\bm{z})\bm{w}$ $\displaystyle=\bm{0}$ where
$\displaystyle\bm{A}(\bm{z})$
$\displaystyle=\left[\begin{array}[]{ccccccccc}y_{0}^{+}z_{0}^{s}&\ldots&y_{0}^{+}&-z_{0}^{n+s-1}&\ldots&-z_{0}^{n}&-z_{0}^{s-1}&\ldots&-1\\\
&\vdots&&\vdots&&\vdots&&\vdots&\\\
y_{m-1}^{+}z_{m-1}^{s}&\ldots&y_{m-1}^{+}&-z_{m-1}^{n+s-1}&\ldots&-z_{m-1}^{n}&-z_{m-1}^{s-1}&\ldots&-1\end{array}\right]$
(31)

with $\bm{A}(\bm{z})$ being of dimension $m\times(3s+1)$ and
$\bm{w}=[\bm{v};\hat{\bm{u}};\tilde{\bm{u}}]$. Assume that $S$ values from
$g_{0},\ldots,g_{s-1}$ are non-zero; hence,
$v^{+}(z),\hat{u}^{+}(z),\tilde{u}^{+}(z)$ are of degrees $S,S-1,S-1$,
respectively. Since
$v(z)=v^{+}(z),\hat{u}(z)=\hat{u}^{+}(z),\tilde{u}(z)=\tilde{u}^{+}(z)$ solve
(27), $\bm{w}=\bm{w}^{+}$ with

$\displaystyle\bm{w}^{+}=[\bm{0}_{s-S};\bm{v}^{+};\bm{0}_{s-S};\hat{\bm{u}}^{+};\bm{0}_{s-S};\tilde{\bm{u}}^{+}]$
(32)

solves (II-A) and vice versa. This implies that $\bm{A}(\bm{z})$ cannot be
full rank and thus it has a null space of at least one dimension. If
$\bm{A}(\bm{z})$ has a one-dimensional null space, meaning that all solutions
to (II-A) are of the form $\bm{w}=c\bm{w}^{+}$ for any complex scalar $c$,
then (27) has the polynomial solution $v(z)=cv^{+}(z)$,
$\hat{u}(z)=c\hat{u}^{+}(z)$, $\tilde{u}(z)=c\tilde{u}^{+}(z)$ and vice versa.

In the special case of $z_{j}=e^{i2\pi j/n}e^{i\gamma/n}$ and $m\leq n$, i.e.,
rotated $n$th roots of unity which we for brevity refer to as shifted
harmonics (or $\gamma$-shifted harmonics if needing to explicitly state
$\gamma$) since the same $n$th root is used throughout this work, we have
$z_{j}^{n}=e^{i\gamma}$ so that (27) becomes

$\displaystyle
y^{+}(z_{j})v(z_{j})-e^{i\gamma}\hat{u}(z_{j})-\tilde{u}(z_{j})$
$\displaystyle=0,\quad 0\leq j\leq m-1$ (33)

which corresponds to solving the linear system

$\displaystyle\bm{B}(\bm{z})\bm{w}$ $\displaystyle=\bm{0}$ where
$\displaystyle\bm{B}(\bm{z})$
$\displaystyle=\left[\begin{array}[]{cccccc}y_{0}^{+}z_{0}^{s}&\ldots&y_{0}^{+}&-z_{0}^{s-1}&\ldots&-1\\\
&\vdots&&&\vdots&\\\
y_{m-1}^{+}z_{m-1}^{s}&\ldots&y_{m-1}^{+}&-z_{m-1}^{s-1}&\ldots&-1\end{array}\right]$
(37)

with $\bm{B}(\bm{z})$ being of dimension $m\times(2s+1)$ and
$\bm{w}=[\bm{v};e^{i\gamma}\hat{\bm{u}}+\tilde{\bm{u}}]$. As with (27) and
(II-A), a one-dimensional null space of $\bm{B}(\bm{z})$ gives all solutions
to (II-A) as $\bm{w}=c\bm{w}^{+}$ for $c\in\mathbb{C}$ with

$\displaystyle\bm{w}^{+}=[\bm{0}_{s-S};\bm{v}^{+};\bm{0}_{s-S};e^{i\gamma}\hat{\bm{u}}^{+}+\tilde{\bm{u}}^{+}],$
(38)

which is equivalent to the polynomial solution $v(z)=cv^{+}(z)$,
$e^{i\gamma}\hat{u}(z)+\tilde{u}(z)=c(e^{i\gamma}\hat{u}^{+}(z)+\tilde{u}^{+}(z))$
of (33) and vice versa.

Next, we generalize Corollary 2 in [20].

######  Theorem 3.

Assume that $\bm{g}^{+}$ has no zeros. Furthermore, assume that
$z_{j}^{n}=e^{i\gamma}$ for $0\leq j\leq 2s-1$ and $(\theta_{k}^{+})^{n}\neq
e^{-i\gamma}$ for $0\leq k\leq s-1$. Solving (33) (or equivalently (II-A))
results in $v(z)=cv^{+}(z)$ and
$e^{i\gamma}\hat{u}(z)+\tilde{u}(z)=c(e^{i\gamma}\hat{u}^{+}(z)+\tilde{u}^{+}(z))$
for some complex scalar $c$.

###### Proof.

Since $z_{j}$, $0\leq j\leq 2s-1$, are shifted harmonics and
$(\theta_{k}^{+})^{n}\neq e^{-i\gamma}$, it follows that $v^{+}(z_{j})\neq 0$
and $t_{k}((\theta_{k}^{+})^{-1})\neq 0$ for $0\leq j\leq 2s-1$, $0\leq k\leq
s-1$. Inserting $y^{+}_{j}=\frac{u^{+}(z_{j})}{v^{+}(z_{j})}$, $0\leq j\leq
2s-1$, into (33) and using $v^{+}(z_{j})\neq 0$ we get

$\displaystyle\frac{u^{+}(z_{j})}{v^{+}(z_{j})}v(z_{j})-e^{i\gamma}\hat{u}(z_{j})-\tilde{u}(z_{j})$
$\displaystyle=0,\quad 0\leq j\leq 2s-1$
$\displaystyle(e^{i\gamma}\hat{u}(z_{j})+\tilde{u}(z_{j}))v^{+}(z_{j})$
$\displaystyle=u^{+}(z_{j})v(z_{j}),\quad 0\leq j\leq 2s-1.$ (39)

Note from (II-A) that
$u^{+}(z_{j})=e^{i\gamma}\hat{u}^{+}(z_{j})+\tilde{u}^{+}(z_{j})$. Since the
polynomials on both sides of the identity in (II-A) are of degree $2s-1$, the
identity theorem for polynomials implies that

$\displaystyle(e^{i\gamma}\hat{u}(z)+\tilde{u}(z))v^{+}(z)=(e^{i\gamma}\hat{u}^{+}(z)+\tilde{u}^{+}(z))v(z),\forall
z$ (40)

The roots of $v(z)$ are
$(\theta_{0}^{+})^{-1},\ldots,(\theta_{s-1}^{+})^{-1}$. Inserting
$z=(\theta_{k}^{+})^{-1}$, $0\leq k\leq s-1$, into
$e^{i\gamma}\hat{u}(z)+\tilde{u}(z)$ gives
$e^{i\gamma}\hat{u}((\theta_{k}^{+})^{-1})+\tilde{u}((\theta_{k}^{+})^{-1})=g_{k}^{+}t_{k}^{+}((\theta_{k}^{+})^{-1})(e^{i\gamma}(\theta_{k}^{+})^{n}-1)$
which is not 0 since $g_{k}^{+}\neq 0$, $(\theta_{k}^{+})^{n}\neq
e^{-i\gamma}$ and $t_{k}^{+}((\theta_{k}^{+})^{-1})\neq 0$. Hence, from (40)
it follows that the roots of $v(z)$ are roots of $v^{+}(z)$ and thus
$v(z)=cv^{+}(z)$ for some complex scalar $c$. Inserting this relation back
into (40) and dividing away $v^{+}(z)$ we obtain
$e^{i\gamma}\hat{u}(z)+\tilde{u}(z)=c(e^{i\gamma}\hat{u}^{+}(z)+\tilde{u}^{+}(z))$.
∎

Hence, Theorem 3 shows that we can recover $\bm{\theta}^{+}$ after $2s$
measurements with shifted harmonics (if there are more than $2s$ measurements,
it is enough if a subset of $2s$ measurements satisfies the assumptions of the
theorem). Note that there is no constraint on $\theta_{k}^{+},0\leq k\leq
s-1,$ to lie inside $\mathbb{T}$ as in [20]. For a $\bm{z}$ satisfying the
assumptions of Theorem 3, arbitrary $\bm{\theta}^{+}$ and $\bm{g}^{+}$ satisfy
the assumptions of Theorem 3 w.p.1. Moreover, for any
$\bm{\theta}^{+},\bm{g}^{+}$ (with $s$ non-zero elements), it is evidently
easy to find a $\bm{z}$ such that the assumptions of Theorem 3 are satisfied.

From Theorem 3 it is possible to obtain $\bm{g}^{+}$.

######  Proposition 1.

Given the assumptions in Theorem 3, the unique solution to (10) is
$\bm{g}=\bm{g}^{+}$ and it can be found with $O(\mathrm{poly}(s))$ complexity.

###### Proof.

From Theorem 3 one obtains $\bm{\theta}^{+}$ and
$q(z)=e^{i\gamma}\hat{u}(z)+\tilde{u}(z)=c(e^{i\gamma}\hat{u}^{+}(z)+\tilde{u}^{+}(z))$.
It follows that
$q((\theta_{k}^{+})^{-1})=cg_{k}^{+}t_{k}^{+}((\theta_{k}^{+})^{-1})(e^{i\gamma}(\theta_{k}^{+})^{n}-1)$,
$0\leq k\leq s-1$, from which we get $\hat{\bm{g}}=c\bm{g}^{+}$ as

$\displaystyle\hat{g}_{k}=cg_{k}^{+}=\frac{q((\theta_{k}^{+})^{-1})}{t_{k}^{+}((\theta_{k}^{+})^{-1})(e^{i\gamma}(\theta_{k}^{+})^{n}-1)},\quad
0\leq k\leq s-1.$

It follows trivially from (10) that

$\hat{\bm{y}}=\bm{V}(\bm{z})^{T}\bm{V}(\bm{\theta}^{+})\hat{\bm{g}}=c\bm{y}^{+}$

from which we easily obtain $c$ as $c=\hat{y}_{k}/y_{k}^{+}$ (for any $k$,
$0\leq k\leq m-1$, where $y_{k}^{+}\neq 0$) and thereby $\bm{g}^{+}$ as
$\bm{g}^{+}=\hat{\bm{g}}/c$. It is clear that these operations for finding
$\bm{g}^{+}$ are of $O(\mathrm{poly(s)})$ complexity. ∎

Clearly, another way to obtain $\bm{g}^{+}$ is by applying the pseudo-inverse:
$\bm{g}^{+}=\left(\bm{V}(\bm{z})^{T}\bm{V}(\bm{\theta}^{+})\right)^{\dagger}\bm{y}^{+}$.

To remove the constraint of shifted harmonic samples, we consider the more
general equation system in (27). First, we show the following.

######  Theorem 4.

Given the assumptions in Theorem 3 and $z_{j}^{n}\neq e^{i\gamma}$,
$v^{+}(z_{j})\neq 0$ for $2s\leq j\leq 3s-1$ with $\gamma$ from Theorem 3.
Solving (27) (or equivalently (II-A)) gives that $v(z)=cv^{+}(z)$,
$\hat{u}(z)=c\hat{u}^{+}(z)$ and $\tilde{u}^{+}(z)=c\tilde{u}(z)$ for some
complex scalar $c$.

###### Proof.

Inserting $y_{j}^{+}=\frac{u^{+}(z_{j})}{v^{+}(z_{j})}$, $0\leq j\leq 3s-1$,
into (27) and using $v^{+}(z_{j})\neq 0$, we get

$\displaystyle(z_{j}^{n}\hat{u}^{+}(z_{j})+\tilde{u}^{+}(z_{j}))v(z_{j})=(z_{j}^{n}\hat{u}(z_{j})+\tilde{u}(z_{j}))v^{+}(z_{j}),\quad
0\leq j\leq 3s-1.$ (41)

From Theorem 3, it follows that after the $2s$ measurements $z_{j}=e^{i2\pi
j/n}e^{i\gamma/n}$, $0\leq j\leq 2s-1$, we obtain that $v(z)=cv^{+}(z)$ and
$e^{i\gamma}\hat{u}(z)+\tilde{u}(z)=c(e^{i\gamma}\hat{u}^{+}(z)+\tilde{u}^{+}(z))$
hold for some complex scalar $c$. Inserting $v(z_{j})=cv^{+}(z_{j})$ into (41)
for $j\geq 2s$ and dividing both sides of the equality with $v^{+}(z_{j})$, we
get

$\displaystyle c(z_{j}^{n}\hat{u}^{+}(z_{j})+\tilde{u}^{+}(z_{j}))$
$\displaystyle=z_{j}^{n}\hat{u}(z_{j})+\tilde{u}(z_{j}),\quad 2s\leq j\leq
3s-1$ $\displaystyle\Rightarrow z_{j}^{n}(c\hat{u}^{+}(z_{j})-\hat{u}(z_{j}))$
$\displaystyle=-(c\tilde{u}^{+}(z_{j})-\tilde{u}(z_{j})),\quad 2s\leq j\leq
3s-1.$ (42)

Rearranging the equality
$e^{i\gamma}\hat{u}(z)+\tilde{u}(z)=c(e^{i\gamma}\hat{u}^{+}(z)+\tilde{u}^{+}(z))$
as $e^{i\gamma}(c\hat{u}^{+}(z)-\hat{u}(z))=-(c\tilde{u}^{+}(z)-\tilde{u}(z))$
and inserting it into (II-A), it follows that

$\displaystyle(c\hat{u}^{+}(z_{j})-\hat{u}(z_{j}))(z_{j}^{n}-e^{i\gamma})=0,\quad
2s\leq j\leq 3s-1.$ (43)

Utilizing the $s$ samples $z_{j}^{n}\neq e^{i\gamma}$, $2s\leq j\leq 3s-1$,
and the fact that $\hat{u}^{+}(z)$, $\hat{u}(z)$ are polynomials of degree
$s-1$, it follows that $c\hat{u}^{+}(z)=\hat{u}(z)$ holds for all $z$.
Inserting this equality into
$e^{i\gamma}(c\hat{u}^{+}(z)-\hat{u}(z))=-(c\tilde{u}^{+}(z)-\tilde{u}(z))$ it
also follows that $c\tilde{u}^{+}(z)=\tilde{u}(z)$. ∎

Clearly, if there are more than $3s$ samples, it is enough if $3s$ of them
satisfy the assumptions of Theorem 4 (while the rest can be arbitrary) to
obtain $v(z)=cv^{+}(z)$, $\hat{u}(z)=c\hat{u}^{+}(z)$ and
$\tilde{u}^{+}(z)=c\tilde{u}(z)$.

For a $\bm{z}$ satisfying the assumptions of Theorem 4, arbitrary
$\bm{\theta}^{+}$ and $\bm{g}^{+}$ will satisfy the assumptions of Theorem 4
w.p.1. Theorem 4 can also be made to hold for any $\bm{\theta}^{+}$ and
$\bm{g}^{+}$ with $s$ non-zero elements. This can be realized by noting that
for such $\bm{\theta}^{+}$ and $\bm{g}^{+}$, it is easy to find a $\bm{z}$
satisfying the assumptions of Theorem 4. The following result shows that we
can relax $\bm{z}$ to be arbitrary.

######  Theorem 5.

For any $\bm{\theta}^{+}$ and $\bm{g}^{+}$ with $s$ non-zero elements, solving
(27) (or equivalently (II-A)) for arbitrary $\bm{z}$, where $m=3s$ and $n\geq
2s$, gives that $v(z)=cv^{+}(z)$, $\hat{u}(z)=c\hat{u}^{+}(z)$ and
$\tilde{u}^{+}(z)=c\tilde{u}(z)$ (for some complex scalar $c$) w.p.1.

###### Proof.

For any $\bm{\theta}^{+}$, $\bm{g}^{+}$ with $s$ non-zero elements, let
$\hat{\bm{z}}=[\hat{z}_{0},\ldots,\hat{z}_{3s-1}]$ be a vector of samples that
satisfies the assumptions of Theorem 4. Thus, $\bm{A}(\hat{\bm{z}})$ has a
one-dimensional null space w.p.1. Since $\bm{A}(\hat{\bm{z}})$ is a
$3s\times(3s+1)$ matrix, a one-dimensional null space of
$\bm{A}(\hat{\bm{z}})$ implies that $\bm{A}(\hat{\bm{z}})$ has rank $3s$ so
that there exists a $3s\times 3s$ sub-matrix $\bm{M}(\hat{\bm{z}})$ in
$\bm{A}(\hat{\bm{z}})$ with full rank $3s$. $\bm{M}(\hat{\bm{z}})$ is obtained
from $\bm{A}(\hat{\bm{z}})$ by removing one column from
$\bm{A}(\hat{\bm{z}})$. We will show that we can replace $\hat{\bm{z}}$ with
the arbitrary $\bm{z}$ while still maintaining a full rank of $\bm{M}(\bm{z})$
w.p.1. From $y^{+}(z_{j})=u^{+}(z_{j})/v^{+}(z_{j})$, it follows that the
determinant $\left|\bm{M}(\bm{z})\right|$ is a multivariate rational function
$p(\bm{z})/q(\bm{z})$ in the variables $z_{0},\ldots,z_{3s-1}$. Since
$\bm{M}(\hat{\bm{z}})$ is invertible, $p(\hat{\bm{z}})\neq 0$, so $p(\bm{z})$
is not zero. Hence, the zero set
$\mathcal{S}=\\{\bm{z}|\bm{z}\in\mathbb{C}^{3s},p(\bm{z})=0\\}$ is either
empty or is a hypersurface in the ambient space $\mathbb{C}^{3s}$ of dimension
at most $3s-1$. Clearly, an arbitrarily chosen $\bm{z}$ is inside
$\mathcal{S}$ w.p.0, proving the theorem. ∎

Obviously, Theorem 5 holds also when $m\geq 8s-3$, since it can be applied to
any subset of $8s-3$ arbitrary samples. Theorem 5 also gives the following
version of Proposition 1

######  Proposition 2.

Given the assumptions in Theorem 5, the unique solution to (10) is
$\bm{g}=\bm{g}^{+}$ and it can be found with $O(\mathrm{poly}(s))$ complexity.

###### Proof.

From the assumptions in Theorem 5, one obtains $\bm{\theta}^{+}$ and
$\hat{u}(z)=c\hat{u}^{+}(z)$. It follows that
$\hat{u}((\theta_{k}^{+})^{-1})=cg_{k}^{+}t_{k}^{+}((\theta_{k}^{+})^{-1})$,
$0\leq k\leq s-1$, from which we get $\hat{\bm{g}}=c\bm{g}^{+}$ as

$\displaystyle\hat{g}_{k}=cg_{k}^{+}=\frac{u((\theta_{k}^{+})^{-1})}{t_{k}^{+}((\theta_{k}^{+})^{-1})},\quad
0\leq k\leq s-1.$

It follows trivially from (10) that

$\hat{\bm{y}}=\bm{V}(\bm{z})^{T}\bm{V}(\bm{\theta}^{+})\hat{\bm{g}}=c\bm{y}^{+}$

from which we easily obtain $c$ as $c=\hat{y}_{k}/y_{k}^{+}$ (for any $k$,
$0\leq k\leq m-1$, where $y_{k}^{+}\neq 0$) and thereby $\bm{g}^{+}$ as
$\bm{g}^{+}=\hat{\bm{g}}/c$. It is clear that these operations for finding
$\bm{g}^{+}$ are of $O(\mathrm{poly(s)})$ complexity. ∎

In the above results, we assumed that $g_{k}\neq 0$ for $0\leq k\leq s-1$;
i.e., we assumed exact knowledge of the number of non-zero elements in
$\bm{g}$. If $g_{k}=0$ for some $k$, $\bm{A}(\bm{z})$ and $\bm{B}(\bm{z})$ can
have null spaces of more than one dimension. Assume that $S$, $S\leq s$, of
the $g_{0},\ldots,g_{s-1}$ are non-zero. One can now iterate through different
dimensions of $\bm{A}(\bm{z})$ and $\bm{B}(\bm{z})$, reducing the value of $s$
in each iteration by 1 until a one-dimensional null space is encountered.
Since for each such dimension of $\bm{A}(\bm{z})$ and $\bm{B}(\bm{z})$,
$\bm{w}^{+}$ in (32) and (38) belongs to their null space, respectively, the
solution will not be missed and the iteration terminates after at most $s-S+1$
steps (when $s$ is reduced to $S$, giving the unique solution $c\bm{w}^{+}$).
Hence, we can formulate the following recovery algorithm

R1-Alg:

  1. Input: $s$, $n\geq 2s$, $\bm{y}^{+}$, samples $\bm{z}\in\mathbb{C}^{m}$ with $m\geq 2s$.

  2. Assumptions: If $2s\leq m<3s$, then $m\leq n$ and all samples in $\bm{z}$ are shifted harmonics. If $m\geq 3s$, $\bm{z}$ can be arbitrary.

  3. Output: $\bm{\theta}^{+}$, $\bm{g}^{+}$.

  4. 1.

If $2s\leq m<3s$, let $\bm{C}=\bm{B}(\bm{z})$ from (II-A); otherwise, let
$\bm{C}=\bm{A}(\bm{z})$ from (II-A).

  5. 2.

Keep reducing the value of $s$ (and thus the dimension of $\bm{C}$) by one
until $\bm{C}$ has a one-dimensional null space (assume $s=S$ when this occurs
for the first time) of the form $c[\bm{v}^{+};\bm{q}]$ for some complex scalar
$c$ and vector $\bm{q}$.

  6. 3.

With vector $c\bm{v}^{+}$ construct the polynomial $cv^{+}(z)$ which has roots
$(\theta_{j}^{+})^{-1}$, $0\leq j\leq S-1$, giving us $\bm{\theta}^{+}$.

  7. 4.

Compute
$\bm{g}^{+}=\left(\bm{V}(\bm{z})^{T}\bm{V}(\bm{\theta}^{+})\right)^{\dagger}\bm{y}^{+}$
or use the procedures outlined in the proofs of Proposition 1 and Proposition
2.

  8. 5.

Output $\bm{\theta}^{+},\bm{g}^{+}$.

For $m=O(s)$, it is obvious that the complexity of all these operations is
$O(\mathrm{poly}(s))$.

###  II-B Derivation of R2

Choose $\bm{z}\in\mathbb{C}^{m}$ and $\bm{\theta}\in\mathbb{C}^{n}$ such that
$\bm{z}$ and any subset of $s$ elements from $\bm{\theta}$ satisfy assumptions
in Theorem 3 or Theorem 4. As seen from these theorems, it is easy to choose
such $\bm{z}$ and $\bm{\theta}$. Construct $\bm{A}$ in R2 as
$\bm{A}=\bm{V}(\bm{z})^{T}\bm{V}(\bm{\theta})$. In R2, denote the (at most)
$s$ non-zero elements in $\bm{x}$ as vector $\bm{g}$ and assume they are on
positions $p_{0},\ldots,p_{s-1}$. Let $\bm{p}=[p_{0},\ldots,p_{s-1}]$ and
$\bm{\theta}_{\bm{p}}=[\theta_{p_{0}},\ldots,\theta_{p_{s-1}}]$. Hence, it
follows that
$\bm{A}\bm{x}=\bm{V}(\bm{z})^{T}\bm{V}(\bm{\theta})\bm{x}=\bm{V}(\bm{z})^{T}\bm{V}(\bm{\theta}_{\bm{p}})\bm{g}$.
Thus, using R1 we obtain $\bm{\theta}_{\bm{p}}$ (from which we infer the
positions $\bm{p}$) and $\bm{g}$ from which we easily reconstruct $\bm{x}$.
Essentially, what is done in R2 is to first create a grid
$\bm{\theta}=[\theta_{0},\ldots,\theta_{n-1}]$ through which $\bm{x}$ is
measured. Thereafter, by using R1, we are able to find the $s$ activated grid
elements which gives us the support of $\bm{x}$ from which we easily recover
its elements through a simple matrix inversion. Basically, we have encoded the
support of $\bm{x}$ onto the grid elements in $\bm{\theta}$ so that recovering
the support becomes equivalent to recovering the active grid elements. Thus,
an algorithm with $O(\mathrm{poly}(s))$ complexity that recovers $\bm{x}$ in
R2 can be summarized as

R2-Alg:

  1. Input: $s$, $n\geq 2s$, $\bm{y}$, $\bm{z}\in\mathbb{C}^{m}$ with $m\geq 2s$, $\bm{\theta}\in\mathbb{C}^{n}$, $\bm{A}=\bm{V}(\bm{z})^{T}\bm{V}(\bm{\theta})$.

  2. Assumptions: If $2s\leq m<3s$, then $m\leq n$ and all samples in $\bm{z}$ are shifted harmonics; furthermore, $\bm{z}$ and $\bm{\theta}$ are chosen such that $\bm{z}$ and each subset of $s$ samples from $\bm{\theta}$ satisfy Theorem 3. If $m\geq 3s$, $\bm{z}$ and $\bm{\theta}$ are chosen such that some $3s$ samples from $\bm{z}$ and each subset of $s$ samples from $\bm{\theta}$ satisfy Theorem 4.

  3. Output: $\bm{x}$.

  4. 1.

If $2s\leq m<3s$, let $\bm{C}=\bm{B}(\bm{z})$ from (II-A); otherwise, let
$\bm{C}=\bm{A}(\bm{z})$ from (II-A).

  5. 2.

Keep reducing the value of $s$ (and thus the dimension of $\bm{C}$) by one
until $\bm{C}$ has a one-dimensional null space (assume $s=S$ when this occurs
for the first time) of the form $c[\bm{v}^{+};\bm{q}]$ for some complex scalar
$c$ and vector $\bm{q}$.

  6. 3.

From vector $c\bm{v}^{+}$ construct the polynomial $cv^{+}(z)$ which has roots
$r_{0},\ldots,r_{S-1}$.

  7. 4.

Compute

$p_{j}=\arg\min_{k}|r_{j}-\theta_{k}^{-1}|,\quad 0\leq j\leq S-1$

and let $\bm{\theta}_{\bm{p}}=[\theta_{p_{0}},\ldots,\theta_{p_{S-1}}]$.

  8. 5.

Compute
$\bm{g}=\left(\bm{V}(\bm{z})^{T}\bm{V}(\bm{\theta}_{\bm{p}})\right)^{\dagger}\bm{y}$
or use the procedures outlined in the proofs of Proposition 1 and Proposition
2.

  9. 6.

Let $\bm{x}=\bm{0}_{n}$ and set $x_{p_{k}}=g_{k}$, $0\leq k\leq S-1$. Output
$\bm{x}$.

Since R2-Alg is essentially the same as R1-Alg, the complexity of R2-Alg is
also $O(\mathrm{poly}(s))$.

##  III Phase-less Measurements

In this section, we look at phase-less measurements. Related to result R2 is

R3: There exists a construction of an $m\times n$ measurement matrix $\bm{A}$,
with $m\geq 4s$ and $n\geq 4s-1$, such that there is a unique
$\bm{x}\in\mathbb{C}_{\leq s}^{n}$ satisfying

$\displaystyle\bm{y}=\left|\bm{A}\bm{x}\right|^{2}.$ (44)

Furthermore, there exists a corresponding recovery algorithm such that

  * •

If $4s\leq m<\min\\{n+1,8s-2\\}$, the recovery complexity is
$O(\mathrm{poly}(s)2^{s})$.

  * •

If $m\geq 8s-2$, the recovery complexity is $O(\mathrm{poly}(s))$.

This result provides a direct solution to open problem 3.5 in [18], as
explained in Section I. To obtain R3, we begin by considering the phase-less
equivalent of R1 with a Fourier matrix

R4: Assume an $n$-dimensional vector $\bm{x}=\bm{F}(\bm{\theta})\bm{g}$, where
$\bm{\theta}\in\mathbb{T}^{s}$, $\bm{g}\in\mathbb{C}^{s}$ and $n\geq 4s-1$.
Vector $\bm{x}$ is being measured with $m$ phase-less linear measurements of
the form

$\displaystyle\bm{y}=\left|\bm{F}(\bm{z})^{T}\bm{x}\right|^{2}=\left|\bm{F}(\bm{z})^{T}\bm{F}(\bm{\theta})\bm{g}\right|^{2}$
(45)

for some $\bm{z}\in\mathbb{T}^{m}$. We have the following results

  * •

If $4s-1\leq m<\min\\{n+1,8s-3\\}$ and samples in $\bm{z}$ are shifted
harmonics, $\bm{\theta}$ and $|\bm{g}|$ (up to a global real scalar) can be
exactly recovered with $O(\mathrm{poly}(s))$ complexity and there are
$2^{s-1}$ different $\bm{g}$ satisfying (45).

  * •

If $m\geq 8s-3$ and $\bm{z}$ is arbitrary, $\bm{\theta}$ and $|\bm{g}|$ (up to
a global real scalar) can be exactly recovered with $O(\mathrm{poly}(s))$
complexity and there are only two different $\bm{g}$ satisfying (45).

To derive R4 we will start (as with R1) by first recovering $\bm{\theta}$ from
(45). However, in contrast to R1, we can no longer recover $\bm{g}$ with a
simple matrix inversion due to the non-linear magnitude operator. After
support recovery $\bm{\theta}$, we are (at first glance) essentially back to a
classical phase retrieval problem. Below we will present an approach, similar
in nature to the approach in Section II, from which certain Laurent
polynomials are solved for. It turns out that there is enough information
about $\bm{g}$ embedded in these Laurent polynomials.

R4 shows that there is an ambiguity when measuring as in (45) which cannot be
resolved by simply taking more measurements (increasing $m$). In order to
resolve the ambiguity of which $\bm{g}$ is the correct one, we take one
additional measurement with some arbitrary measurement vector $\bm{a}$.
Therefore, we have the following reformulation of R4

R5: Assume an $n$-dimensional vector $\bm{x}=\bm{F}(\bm{\theta})\bm{g}$, where
$\bm{\theta}\in\mathbb{T}^{s}$, $\bm{g}\in\mathbb{C}^{s}$ and $n\geq 4s-1$.
Vector $\bm{x}$ is being observed through $m+1$ phase-less linear measurements
of the form

$\displaystyle\bm{y}=\left|\bm{F}(\bm{z})^{T}\bm{x}\right|^{2}=\left|\bm{F}(\bm{z})^{T}\bm{F}(\bm{\theta})\bm{g}\right|^{2}$
$\displaystyle
y_{m}=|\bm{a}^{T}\bm{x}|^{2}=\left|\bm{a}^{T}\bm{F}(\bm{\theta})\bm{g}\right|^{2}$
(46)

for some $\bm{z}\in\mathbb{T}^{m}$ where $\bm{a}$ is an arbitrary vector. We
have the following results

  * •

If $4s-1\leq m<\min\\{n+1,8s-3\\}$ and samples in $\bm{z}$ are shifted
harmonics, $\bm{\theta}$ and $\bm{g}$ can be recovered exactly with
$O(\mathrm{poly}(s)2^{s})$ complexity.

  * •

If $m\geq 8s-3$ and $\bm{z}$ is arbitrary, $\bm{\theta}$ and $\bm{g}$ can be
recovered exactly with $O(\mathrm{poly}(s))$ complexity.

###  III-A Derivation of R4 and R5

Sample $y_{j}$ in R4 is now seen as a phase-less sample of the rational
function $f(z)$ in (11)

$\displaystyle
y_{j}=|f(z_{j})|^{2}=\left|\frac{z_{j}^{n}\hat{u}(z_{j})+\tilde{u}(z_{j})}{v(z_{j})}\right|^{2}.$
(47)

Let $\overline{(.)}$ denote the complex conjugate operator. In order to deal
with the non-linear magnitude operation, we express it as

$\displaystyle
y(z)=\left(\frac{z^{n}\hat{u}(z)+\tilde{u}(z)}{v(z)}\right)\left(\frac{\overline{z^{n}\hat{u}(z)+\tilde{u}(z)}}{{\bar{v}(z)}}\right)=\left(\frac{z^{n}\hat{u}(z)+\tilde{u}(z)}{v(z)}\right)\left(\frac{z^{-n}\overline{\hat{u}(z)}+\overline{\tilde{u}(z)}}{{\overline{v(z)}}}\right).$
(48)

Henceforth, it is assumed that $z\in\mathbb{T}$ and
$\bm{\theta}\in\mathbb{T}^{s}$. Thus, $\bar{z}=1/z$ and (48) becomes

$\displaystyle
y(z)\stackrel{{\scriptstyle\triangle}}{{=}}\frac{\hat{u}(z)\overline{\hat{u}(z)}+\tilde{u}(z)\overline{\tilde{u}(z)}+z^{-n}\overline{\hat{u}(z)}\tilde{u}(z)+z^{n}\hat{u}(z)\overline{\tilde{u}(z)}}{v(z)\overline{v(z)}}.$
(49)

Due to $\overline{z}=1/z$, the number of unknowns is reduced to just one ($z$)
and the numerator and denominator in (49) become Laurent polynomials in
$z\in\mathbb{T}$. Defining the Laurent polynomials

$\displaystyle L(z)$
$\displaystyle\stackrel{{\scriptstyle\triangle}}{{=}}\hat{u}(z)\overline{\hat{u}(z)}+\tilde{u}(z)\overline{\tilde{u}(z)}=|\hat{u}(z)|^{2}+|\tilde{u}(z)|^{2}$
$\displaystyle\tilde{L}(z)$
$\displaystyle\stackrel{{\scriptstyle\triangle}}{{=}}\hat{u}(z)\overline{\tilde{u}(z)}$
$\displaystyle\hat{L}(z)$
$\displaystyle\stackrel{{\scriptstyle\triangle}}{{=}}v(z)\overline{v(z)}$ (50)

where $L(z)$ and $\tilde{L}(z)$ are of at most degree $s-1$ in $z$ and $1/z$,
respectively, while $\hat{L}_{s}(z)$ is of at most degree $s$ in $z$ and
$1/z$, we can express $y(z)$ as a ratio of Laurent polynomials

$\displaystyle
y(z)=\frac{L(z)+z^{n}\tilde{L}(z)+z^{-n}\overline{\tilde{L}(z)}}{\hat{L}(z)}.$
(51)

Two properties of the introduced Laurent polynomials, important for this work,
are given next.

######  Proposition 3.

A Laurent polynomial $P(z)$ with smallest degree $-s$ and largest degree $t$
can be expressed, for $z\neq 0$, as $\frac{1}{z^{s}}Q(z)$ where $Q(z)$ is a
complex polynomial of degree $t+s$; thus, $P(z)$ and $Q(z)$ have the same non-
zero roots.

###### Proof.

Write
$P(z)=\sum_{k=-s}^{t}p_{k}z^{k}=\frac{1}{z^{s}}\sum_{k=0}^{s+t}p_{k-s}z^{k}=\frac{1}{z^{s}}Q(z)$
where $Q(z)=\sum_{k=0}^{s+t}p_{k-s}z^{k}$. ∎

######  Proposition 4.

Given two complex polynomials $p(z)$ and $r(z)$ with non-zero roots
$p_{0},\ldots,p_{k-1}$ and $r_{0},\ldots,r_{l-1}$, respectively. Compute the
Laurent polynomial $Q(z)=p(z)\overline{r(z)}$ with the assumption that
$z\in\mathbb{T}$. It follows that the roots of $Q(z)$ are
$p_{0},\ldots,p_{k-1}$ and $1/\overline{r_{0}},\ldots,1/\overline{r_{l-1}}$.

###### Proof.

The assumption $z\in\mathbb{T}$ gives that $\overline{r(z)}$ can equivalently
be expressed as
$\overline{r(z)}=\prod_{j=0}^{l-1}\overline{(z-r_{j})}=\prod_{j=0}^{l-1}\left(\frac{1}{z}-\overline{r_{j}}\right)=\frac{1}{z^{l}}\prod_{j=0}^{l-1}\left(1-z\overline{r_{j}}\right)=\frac{1}{z^{l}}\hat{r}(z)$
where $\hat{r}(z)$ has $1/\overline{r_{0}},\ldots,1/\overline{r_{l-1}}$ as
roots. Hence, we can write $Q(z)$ as $Q(z)=\frac{1}{z^{l}}p(z)\hat{r}(z)$
which completes the proof. ∎

Proposition 4 shows that the roots of $\hat{L}(z)$,
$1/\theta_{k}=\overline{\theta_{k}}$, $0\leq k\leq s-1$, appear with
multiplicity two.

When $n>2s-2$, the Laurent polynomial in the numerator of (51) contains at
most $3(2s-1)=6s-3$ parameters which is independent of $n$. Furthermore,
$\hat{L}(z)$ contains at most $2s+1$ parameters. Hence, in total, $y(z)$ is
described by $6s-3+2s+1=8s-2$ parameters. Thus, as in [20], we can hope that
$8s-3$ values of $z$ from $\mathbb{T}$ are enough for recovering the $s$
values in $\bm{\theta}$ that identify the support of $\bm{x}$ (up to a scalar,
due to one less measurement than number of unknowns, which is irrelevant for
support identification).

Let $\bm{l},\tilde{\bm{l}},\hat{\bm{l}}$ denote the coefficient vectors of the
Laurent polynomials $L(z),\tilde{L}(z),\hat{L}(z)$ with degrees $s,s-1,s-1$ in
$z$ and $1/z$, respectively. Assume that $y^{+}_{j}=f^{+}(z_{j})$, $0\leq
j\leq m-1$, are the observed measurements which are produced by Laurent
polynomials $L^{+}(z),\tilde{L}^{+}(z),\hat{L}^{+}(z)$ corresponding to
realizations $\bm{g}=\bm{g}^{+}$ and $\bm{\theta}=\bm{\theta}^{+}$. Consider
the equation system

$\displaystyle
y^{+}_{j}\hat{L}(z_{j})-L(z_{j})-z_{j}^{n}\tilde{L}(z_{j})-z_{j}^{-n}\overline{\tilde{L}(z_{j})}=0,\quad
0\leq j\leq m-1$ (52)

in the vectors $\bm{l},\tilde{\bm{l}},\hat{\bm{l}}$. Equation system (52)
corresponds to the following linear equation system

$\displaystyle\bm{G}(\bm{z})\bm{w}$ $\displaystyle=\bm{0}$ where
$\displaystyle\bm{G}(\bm{z})$
$\displaystyle=\left[\begin{array}[]{cccccc}\bm{B}(\bm{z})&\bm{y}^{+}&\mathrm{fliplr}\\{\overline{\bm{B}(\bm{z})}\\}&-\bm{C}(\bm{z})&-\bm{1}&-\mathrm{fliplr}\\{\overline{\bm{C}(\bm{z})}\\}\end{array}\right]$
(54) $\displaystyle\bm{B}(\bm{z})$
$\displaystyle=\left[\begin{array}[]{cccc}y_{0}^{+}z_{0}^{s}&y_{0}^{+}z_{0}^{s-1}&\ldots&y_{0}^{+}z_{0}\\\
&\vdots&\vdots&\\\
y^{+}_{m-1}z_{m-1}^{s}&y^{+}_{m-1}z_{m-1}^{s-1}&\ldots&y^{+}_{m-1}z_{m-1}\end{array}\right]$
(58) $\displaystyle\bm{C}(\bm{z})$
$\displaystyle=\left[\begin{array}[]{cccccccc}z_{0}^{n+s-1}&z_{0}^{n+s-2}&\ldots&z_{0}^{n-s+1}&z_{0}^{s-1}&z_{0}^{s-2}&\ldots&z_{0}\\\
&\vdots&&\vdots&&\vdots&\\\
z_{m-1}^{n+s-1}&z_{m-1}^{n+s-2}&\ldots&z_{m-1}^{n-s+1}&z_{m-1}^{s-1}&z_{m-1}^{s-2}&\ldots&z_{m-1}\end{array}\right].$
(62)

where the operator $\mathrm{fliplr}\\{.\\}$ flips a matrix vertically (i.e.,
flips its columns), $\bm{G}(\bm{z})$ has dimensions $m\times(8s-2)$ and
$\bm{w}=[\hat{\bm{l}};\tilde{\bm{l}};\bm{l};\overline{\tilde{\bm{l}}}]$.
Assume that $S$ values from $g_{0},\ldots,g_{s-1}$ are non-zero. Since
$\hat{L}(z)=\hat{L}^{+}(z),L(z)=L^{+}(z),\tilde{L}(z)=\tilde{L}^{+}(z)$ solves
(52),

$\displaystyle\bm{w}^{+}=[\bm{0}_{s-S};\hat{\bm{l}}^{+};\bm{0}_{s-S};\bm{0}_{s-S};\tilde{\bm{l}}^{+};\bm{0}_{s-S};\bm{0}_{s-S};\bm{l}^{+};\bm{0}_{s-S};\bm{0}_{s-S};\overline{\tilde{\bm{l}}^{+}};\bm{0}_{s-S}]$
(63)

solves (III-A) and vice versa. This implies that $\bm{G}(\bm{z})$ cannot be
full rank and therefore has a null space of at least one dimension. If
$\bm{G}(\bm{z})$ has a one-dimensional null space, meaning that (III-A) has
the solution $\bm{w}=c\bm{w}^{+}$ for any real scalar $c$, then (52) has the
unique solution
$\hat{L}(z)=c\hat{L}^{+}(z),L(z)=cL^{+}(z),\tilde{L}(z)=c\tilde{L}^{+}(z)$ and
vice versa.

In the special case of $z_{j}=e^{i2\pi j/n}e^{i\gamma/n}$ and $m\leq n$, (52)
reduces to

$\displaystyle
y^{+}_{j}\hat{L}(z_{j})-L(z_{j})-e^{i\gamma/n}\tilde{L}(z_{j})-e^{-i\gamma/n}\overline{\tilde{L}(z_{j})}=0,\quad
z_{j}=e^{i2\pi j/n}e^{i\gamma/n},0\leq j\leq m-1$ (64)

which corresponds to the following linear equation system

$\displaystyle\tilde{\bm{G}}(\bm{z})\bm{w}$ $\displaystyle=\bm{0}$ where
$\displaystyle\tilde{\bm{G}}(\bm{z})$
$\displaystyle=\left[\begin{array}[]{cccccc}\bm{B}(\bm{z})&\bm{y}^{+}&\mathrm{fliplr}\\{\overline{\bm{B}(\bm{z})}\\}&-\tilde{\bm{C}}(\bm{z})&-\bm{1}&-\mathrm{fliplr}\\{\overline{\tilde{\bm{C}}(\bm{z})}\\}\end{array}\right]$
(66) $\displaystyle\bm{B}(\bm{z})$
$\displaystyle=\left[\begin{array}[]{cccc}y_{0}^{+}z_{0}^{s}&y_{0}^{+}z_{0}^{s-1}&\ldots&y_{0}^{+}z_{0}\\\
&\vdots&\vdots&\\\
y^{+}_{m-1}z_{m-1}^{s}&y^{+}_{m-1}z_{m-1}^{s-1}&\ldots&y^{+}_{m-1}z_{m-1}\end{array}\right]$
(70) $\displaystyle\tilde{\bm{C}}(\bm{z})$
$\displaystyle=\left[\begin{array}[]{cccc}z_{0}^{s-1}&z_{0}^{s-2}&\ldots&z_{0}\\\
&&\vdots&\\\ z_{m-1}^{s-1}&z_{m-1}^{s-2}&\ldots&z_{m-1}\end{array}\right].$
(74)

where $\tilde{\bm{G}}(\bm{z})$ has dimensions $m\times 4s$ and
$\bm{w}=[\hat{\bm{l}};\bm{l}+e^{i\gamma/n}\tilde{\bm{l}}+e^{-i\gamma/n}\overline{\tilde{\bm{l}}}]$.
As with (52) and (III-A), a one-dimensional null space of
$\tilde{\bm{G}}(\bm{z})$ gives all solutions to (III-A) as
$\bm{w}=c\bm{w}^{+}$, where

$\displaystyle\bm{w}^{+}=[\bm{0}_{s-S};\hat{\bm{l}};\bm{0}_{s-S};\bm{0}_{s-S};\bm{l}+e^{i\gamma/n}\tilde{\bm{l}}+e^{-i\gamma/n}\overline{\tilde{\bm{l}}};\bm{0}_{s-S}],$
(75)

for any real scalar $c$, which is equivalent to the solution
$\hat{L}(z)=c\hat{L}^{+}(z)$,
$L(z)+e^{i\gamma/n}\tilde{L}(z)+e^{-i\gamma/n}\overline{\tilde{L}(z)}=c\left(L^{+}(z)+e^{i\gamma/n}\tilde{L}^{+}(z)+e^{-i\gamma/n}\overline{\tilde{L}(z)^{+}}\right)$
of (64) and vice versa.

Before proceeding to solve (52) and (64), we show that from each solution to
these equations we can construct another solution to the equations.

######  Theorem 6.

Assume that $\theta_{i}\in\mathbb{T}$, $0\leq i\leq s-1$, are given. If
$\bm{g}^{a}$ is a solution to (52) and (64), then $\bm{g}^{b}$ with
$g_{l}^{b}=\overline{g_{l}^{a}}\theta_{l}^{-n}\left(\prod_{\begin{subarray}{c}i=0\\\
i\neq l\end{subarray}}^{s-1}\overline{\theta_{i}}\right)$, $0\leq l\leq s-1$,
is also a solution to (52) and (64).

###### Proof.

With $\bm{g}^{a}$ as a solution to (52) and (64), the coefficients of
$\hat{u}(z),\tilde{u}(z)$ are $g_{k}=g_{k}^{a}$. From the expressions of
$L(z),\tilde{L}(z)$ in (III-A), it follows that if we let
$\hat{u}^{c}(z)=\overline{\tilde{u}(z)}$ and
$\tilde{u}^{c}(z)=\overline{\hat{u}(z)}$, then these polynomials also satisfy
$L(z)=|\hat{u}^{c}(z)|^{2}+|\tilde{u}^{c}(z)|^{2}$ and
$\tilde{L}(z)=\hat{u}^{c}(z)\overline{\tilde{u}^{c}(z)}$. Since $\hat{L}(z)$
is given, due to the assumption that $\theta_{i}$, $0\leq i\leq s-1$, are
given, $\hat{u}^{c}(z)$ and $\tilde{u}^{c}(z)$ are also a solution to (52) and
(64). Next, we find the expressions for $\hat{u}^{c}(z),\tilde{u}^{c}(z)$ from
which we infer $g_{k}^{b}$. The definitions in (II-A) give

$\displaystyle\hat{u}^{c}(z)$
$\displaystyle=-\sum_{l=0}^{s-1}\overline{g^{a}_{l}t_{l}(z)}=-\sum_{l=0}^{s-1}\overline{g^{a}_{l}}\prod_{\begin{subarray}{c}i=0\\\
i\neq l\end{subarray}}^{s-1}(\bar{z}\overline{\theta_{i}}-1)$
$\displaystyle=-\sum_{l=0}^{s-1}\overline{g^{a}_{l}}\prod_{\begin{subarray}{c}i=0\\\
i\neq l\end{subarray}}^{s-1}\left(\frac{1}{z\theta_{i}}-1\right)$ (76)
$\displaystyle=-\frac{(-1)^{s-1}}{z^{s-1}}\sum_{l=0}^{s-1}\left(\overline{g^{a}_{l}}\prod_{\begin{subarray}{c}i=0\\\
i\neq
l\end{subarray}}^{s-1}\overline{\theta_{i}}\right)\prod_{\begin{subarray}{c}i=0\\\
i\neq l\end{subarray}}^{s-1}\left(z\theta_{i}-1\right)$ (77)
$\displaystyle=-\frac{(-1)^{s-1}}{z^{s-1}}\sum_{l=0}^{s-1}\left(\overline{g^{a}_{l}}\theta_{l}^{-n}\prod_{\begin{subarray}{c}i=0\\\
i\neq
l\end{subarray}}^{s-1}\overline{\theta_{i}}\right)\theta_{l}^{n}t_{l}(z).$

Equality (76) holds because $z,\theta_{i}\in\mathbb{T}$ while (77) follows
from extracting the term $\frac{-1}{z\theta_{i}}$ from the product in (76) and
using the fact that $\frac{1}{\theta_{i}}=\overline{\theta_{i}}$. Similarly,
from (II-A), we have

$\displaystyle\tilde{u}^{c}(z)$
$\displaystyle=\sum_{l=0}^{s-1}\overline{g^{a}_{l}}\theta_{l}^{-n}\overline{t_{l}(z)}=\sum_{l=0}^{s-1}\overline{g^{a}_{l}}\theta_{l}^{-n}\prod_{\begin{subarray}{c}i=0\\\
i\neq l\end{subarray}}^{s-1}(\bar{z}\overline{\theta_{i}}-1)$
$\displaystyle=\sum_{l=0}^{s-1}\overline{g^{a}_{l}}\theta_{l}^{-n}\prod_{\begin{subarray}{c}i=0\\\
i\neq l\end{subarray}}^{s-1}\left(\frac{1}{z\theta_{i}}-1\right)$
$\displaystyle=\frac{(-1)^{s-1}}{z^{s-1}}\sum_{l=0}^{s-1}\left(\overline{g^{a}_{l}}\theta_{l}^{-n}\prod_{\begin{subarray}{c}i=0\\\
i\neq
l\end{subarray}}^{s-1}\overline{\theta_{i}}\right)\prod_{\begin{subarray}{c}i=0\\\
i\neq l\end{subarray}}^{s-1}\left(z\theta_{i}-1\right)$
$\displaystyle=\frac{(-1)^{s-1}}{z^{s-1}}\sum_{l=0}^{s-1}\left(\overline{g^{a}_{l}}\theta_{l}^{-n}\prod_{\begin{subarray}{c}i=0\\\
i\neq l\end{subarray}}^{s-1}\overline{\theta_{i}}\right)t_{l}(z).$

Hence, $\hat{u}^{c}(z)$ and $\tilde{u}^{c}(z)$ are of the same form as
$\frac{(-1)^{s-1}}{z^{s-1}}\hat{u}(z)$ and
$\frac{(-1)^{s-1}}{z^{s-1}}\tilde{u}(z)$, respectively, from where it is also
seen that $\hat{u}^{c}(z)$ and $\tilde{u}^{c}(z)$ produce the same Laurent
polynomials as $\hat{u}(z),\tilde{u}(z)$ when $z\in\mathbb{T}$. Thus,
replacing $g^{a}_{l}$ in $\hat{u}(z)$ and $\tilde{u}(z)$ with
$g^{b}_{l}=\left(\overline{g^{a}_{l}}\theta_{l}^{-n}\prod_{\begin{subarray}{c}i=0\\\
i\neq l\end{subarray}}^{s-1}\overline{\theta_{i}}\right)$ results in the same
Laurent polynomials. ∎

From the expression of $\bm{g}^{b}$ in Theorem 6, we get

$\displaystyle\overline{g_{l}^{b}}\theta_{l}^{-n}\left(\prod_{\begin{subarray}{c}i=0\\\
i\neq
l\end{subarray}}^{s-1}\bar{\theta}_{i}\right)=g_{l}^{a}\theta_{l}^{n}\theta_{l}^{-n}\left(\prod_{\begin{subarray}{c}i=0\\\
i\neq
l\end{subarray}}^{s-1}\theta_{i}\right)\left(\prod_{\begin{subarray}{c}i=0\\\
i\neq l\end{subarray}}^{s-1}\overline{\theta_{i}}\right)=g_{l}^{a}$ (78)

where the second equality holds due to $\theta_{i}\in\mathbb{T}$, $0\leq i\leq
s-1$. Hence, $\bm{g}^{a}$ can be obtained from $\bm{g}^{b}$ with the same
transform from Theorem 6 that produces $\bm{g}^{b}$ from $\bm{g}^{a}$.
Therefore, we will refer to $\bm{g}^{b}$ as being the dual of $\bm{g}^{a}$ and
vice versa. Clearly, the solutions to (52) and (64) can be paired as
$(\bm{g}^{a},\bm{g}^{b})$, with $\bm{g}^{a}$ and $\bm{g}^{b}$ being duals of
each other.

We start by providing a solution to (III-A). As described, we first
investigate the minimum number of samples needed to recover the support
$\bm{\theta}^{+}$.

######  Theorem 7.

Assume that $n\geq 4s-1$, $z_{j}^{n}=e^{i\gamma}$ for $0\leq j\leq 4s-2$,
$g_{k}^{+}\neq 0$, $\theta_{k}^{+}\in\mathbb{T}$ and $(\theta_{k}^{+})^{n}\neq
e^{-i\gamma}$ for $0\leq k\leq s-1$. Solving (64) (or equivalently (III-A))
results in
$\hat{L}(z)=c\hat{L}^{+}(z),L(z)+z^{n}\tilde{L}(z)+z^{-n}\overline{\tilde{L}}(z)=c\left(L^{+}(z)+z^{n}\tilde{L}^{+}(z)+z^{-n}\overline{\tilde{L}^{+}}(z)\right)$
for some real scalar $c$.

###### Proof.

Inserting
$y^{+}(z_{j})=\left(L^{+}(z_{j})+z_{j}^{n}\tilde{L}^{+}(z_{j})+z_{j}^{-n}\overline{\tilde{L}^{+}(z_{j})}\right)/\hat{L}^{+}(z_{j})$
into (52) and multiplying both sides of the equation with $\hat{L}^{+}(z_{j})$
results in the equation system

$\displaystyle\left(L^{+}(z_{j})+z_{j}^{n}\tilde{L}^{+}(z_{j})+z_{j}^{-n}\overline{\tilde{L}^{+}(z_{j})}\right)\hat{L}(z_{j})-\left(L(z_{j})+z_{j}^{n}\tilde{L}(z_{j})+z_{j}^{-n}\overline{\tilde{L}(z_{j})}\right)\hat{L}^{+}(z_{j})$
$\displaystyle=0$ (79)

which holds for all $z_{j}$ such that $\hat{L}^{+}(z_{j})\neq 0$. Since
$z_{j}^{n}=e^{i\gamma}$ and $(\theta_{k}^{+})^{n}\neq e^{-i\gamma}$ for $0\leq
j\leq 4s-2$, $0\leq k\leq s-1$, it follows that $\hat{L}^{+}(z_{j})\neq 0$ for
$0\leq j\leq 4s-2$. Using the fact that $z_{j}^{n}=e^{i\gamma}$ for $0\leq
j\leq 4s-2$, (79) gives

$\displaystyle\left(L^{+}(z_{j})+e^{i\gamma}\tilde{L}^{+}(z_{j})+e^{-i\gamma}\overline{\tilde{L}^{+}(z_{j})}\right)\hat{L}(z_{j})$
$\displaystyle-(L(z_{j})+e^{i\gamma}\tilde{L}(z_{j})+e^{-i\gamma}\overline{\tilde{L}(z_{j})})\hat{L}^{+}(z_{j})=0,\quad
0\leq j\leq 4s-2.$ (80)

Since $z_{j}\neq 0$, it follows from Proposition 3 that equation (III-A) is
equivalent to a polynomial equation in $z_{j}$ of degree $4s-2$. Applying the
identity theorem for polynomials to (III-A) gives that the relation

$\displaystyle\left(L^{+}(z)+e^{i\gamma}\tilde{L}^{+}(z)+e^{-i\gamma}\overline{\tilde{L}^{+}(z)}\right)\hat{L}(z)=\left(L(z)+e^{i\gamma}\tilde{L}(z)+e^{-i\gamma}\overline{\tilde{L}(z)}\right)\hat{L}^{+}(z)$
(81)

holds for all $z$. Proposition 4 shows that the roots of $\hat{L}^{+}(z)$ are
$1/\theta_{k}^{+}$, $0\leq k\leq s-1$, where each one of them has multiplicity
two. These roots are not roots of
$L^{+}(z)+e^{i\gamma}\tilde{L}^{+}(z)+e^{-i\gamma}\overline{\tilde{L}^{+}(z)}$.
Assume the contrary. Noting that
$L^{+}(z)+e^{i\gamma}\tilde{L}^{+}(z)+e^{-i\gamma}\overline{\tilde{L}^{+}(z)}=\left|e^{i\gamma}\hat{u}^{+}(z)+\tilde{u}^{+}(z)\right|^{2}$
implies that the roots of $\hat{L}^{+}(z)$ must be roots of
$e^{i\gamma}\hat{u}^{+}(z)+\tilde{u}^{+}(z)$. It follows from (II-A) that

$e^{i\gamma}\hat{u}^{+}\left(\frac{1}{\theta_{k}^{+}}\right)+\tilde{u}^{+}\left(\frac{1}{\theta_{k}^{+}}\right)=g_{k}^{+}t_{k}^{+}\left(\frac{1}{\theta_{k}^{+}}\right)\left(e^{i\gamma}(\theta_{k}^{+})^{n}-1\right).$

Since $(\theta_{k}^{+})^{n}\neq e^{-i\gamma}$ and $g_{k}\neq 0$, it follows
that $t_{k}^{+}\left(\frac{1}{\theta_{k}^{+}}\right)\neq 0$ and thus
$g_{k}^{+}t_{k}^{+}\left(\frac{1}{\theta_{k}^{+}}\right)\neq 0$, which implies
that the roots of $\hat{L}^{+}(z)$ are not roots of
$L^{+}(z)+e^{i\gamma}\tilde{L}^{+}(z)+e^{-i\gamma}\overline{\tilde{L}^{+}(z)}$,
a contradiction. Hence, the roots of $\hat{L}^{+}(z)$ are roots of
$\hat{L}(z)$ which implies that $\hat{L}(z)=c\hat{L}^{+}(z)$ for some real
scalar $c$.

Inserting $\hat{L}(z)=c\hat{L}^{+}(z)$ into (81) and dividing out
$\hat{L}^{+}(z)$, we conclude that

$c\left(L^{+}(z)+e^{i\gamma}\tilde{L}^{+}(z)+e^{-i\gamma}\overline{\tilde{L}^{+}(z)}\right)=L(z)+e^{i\gamma}\tilde{L}(z)+e^{-i\gamma}\overline{\tilde{L}(z)}.$

∎

From the result in Theorem 7, it follows that when $g_{k}\neq 0$, $0\leq k\leq
s-1$, we can find $\hat{L}^{+}(z)$ (up to a real scalar $c$) after $4s-1$
measurements with shifted harmonic samples from which we can recover the
support $\bm{\theta}^{+}$ (below we will deal with the case when $S$, $S\leq
s$, of $g_{0},\ldots,g_{s-1}$ are non-zero). Furthermore, from these
measurements, we also obtain
$L^{+}(z)+e^{i\gamma}\tilde{L}^{+}(z)+e^{-i\gamma}\overline{\tilde{L}^{+}(z)}=\left|e^{i\gamma}\hat{u}^{+}(z)+\tilde{u}^{+}(z)\right|^{2}$
(up to real scalar $c$). For arbitrary $\bm{\theta}^{+}\in\mathbb{T}^{s}$ and
$\bm{g}^{+}$, the assumptions of Theorem 7 are satisfied w.p.1, implying that
$\tilde{\bm{G}}(\bm{z})$ in (III-A) has a one-dimensional null space w.p.1
when the $m\geq 4s-1$ samples in $\bm{z}$ are shifted harmonics. Moreover, for
$\bm{g}^{+}$ with no zero elements and any $\bm{\theta}^{+}\in\mathbb{T}^{s}$,
it is evidently easy to choose a $\bm{z}$ such that the assumptions of Theorem
7 are satisfied.

From Theorem 7 one can find $|\bm{g}^{+}|$ up to a real scalar.

######  Proposition 5.

Given the assumptions in Theorem 7, $|\bm{g}^{+}|$ can be found, up to a real
scalar $c$, with $O(\mathrm{poly}(s))$ complexity.

###### Proof.

From Theorem 7 one obtains $\bm{\theta}^{+}$ and

$q(z)=L(z)+z^{n}\tilde{L}(z)+z^{-n}\overline{\tilde{L}(z)}=c\left(L^{+}(z)+z^{n}\tilde{L}^{+}(z)+z^{-n}\overline{\tilde{L}^{+}(z)}\right)=c|e^{i\gamma}\hat{u}^{+}(z)+\tilde{u}^{+}(z)|^{2}$

for some unknown real scalar $c$. Since
$(\theta_{k}^{+})^{-1}=\overline{\theta_{k}^{+}}$ and

$q(\overline{\theta_{k}^{+}})=c\left|e^{i\gamma}\hat{u}^{+}(\overline{\theta_{k}^{+}})+\tilde{u}^{+}(\overline{\theta_{k}^{+}})\right|^{2}=c\left|g_{k}^{+}t_{k}^{+}\left(\overline{\theta_{k}^{+}}\right)\left(e^{i\gamma}(\theta_{k}^{+})^{n}-1\right)\right|^{2},$

it follows that

$c|g_{k}^{+}|^{2}=\frac{q(\overline{\theta_{k}^{+}})}{\left|t_{k}^{+}\left(\overline{\theta_{k}^{+}}\right)\left(e^{i\gamma}(\theta_{k}^{+})^{n}-1\right)\right|^{2}}.$

Clearly, the operations leading to $c|\bm{g}^{+}|$ are of
$O(\mathrm{poly}(s))$ complexity, showing that $c|\bm{g}^{+}|$ can be found
with $O(\mathrm{poly}(s))$ complexity. ∎

A direct application of this proposition is that we can find the relative
magnitudes of the different frequency components with $4s-1$ measurements and
$O(\mathrm{poly}(s))$ complexity.

When it comes to recovering $\bm{g}^{+}$, the situation is trickier. Given
$\left|e^{i\gamma}\hat{u}^{+}(z)+\tilde{u}^{+}(z)\right|^{2}$, up to a real
scalar, Proposition 4 implies that the roots of
$\left|e^{i\gamma}\hat{u}^{+}(z)+\tilde{u}^{+}(z)\right|^{2}$ appear as
conjugate pairs: if $r$ is a root of
$\left|e^{i\gamma}\hat{u}^{+}(z)+\tilde{u}^{+}(z)\right|^{2}$ then so is
$1/\overline{r}$ as well. Let
$(r_{0},1/\overline{r_{0}}),\ldots,(r_{s-2},1/\overline{r_{s-2}})$ be the
$s-1$ conjugate root pairs. As seen from (II-A),
$e^{i\gamma}\hat{u}^{+}(z)+\tilde{u}^{+}(z)$ is of the same form as
$\tilde{u}^{+}(z)$ in (II-A) but with coefficients
$(e^{i\gamma}(\theta_{k}^{+})^{n}-1)g_{k}^{+}$, $0\leq k\leq s-1$. Hence, for
an arbitrary $\bm{g}^{+}$, Lemma 2 implies that
$(r_{0},1/\overline{r_{0}}),\ldots,(r_{s-2},1/\overline{r_{s-2}})$ are
distinct w.p.1. Only one member from each pair can be a root of
$e^{i\gamma}\hat{u}^{+}(z)+\tilde{u}^{+}(z)$. Since we do not know beforehand
which member of each pair is the actual root of
$e^{i\gamma}\hat{u}^{+}(z)+\tilde{u}^{+}(z)$, we need to go through all the
$2^{s-1}$ possibilities. Denote by $\bm{q}^{k}=[q_{0}^{k},\ldots,q_{s-2}^{k}]$
the $k$:th possibility, meaning that $q_{j}^{k}$ is either $r_{j}$ or
$1/\overline{r_{j}}$ for $0\leq j\leq s-2$, $0\leq k\leq 2^{s-1}-1$. For each
$\bm{q}^{k}$, we get an equation system in $\bm{g}$

$\displaystyle
e^{i\gamma}\tilde{u}(q_{j}^{k})+\hat{u}(q_{j}^{k})=\sum_{\begin{subarray}{c}l=0\end{subarray}}^{s-1}g_{l}t_{l}^{+}(q_{j}^{k})(e^{i\gamma}(\theta_{l}^{+})^{n}-1)=0,\quad
0\leq j\leq s-2.$ (82)

The equation system in (82) with $s-1$ equations represents the zeros of a
polynomial of degree $s-1$ and thus it has a unique solution $\bm{g}^{k}$ up
to a complex scalar. Each $\bm{g}^{k}$, with proper normalization, is a
solution to (64), i.e.,
$\bm{y}^{+}=|\bm{V}^{T}(\bm{z})\bm{V}(\bm{\theta}^{+})\bm{g}^{k}|^{2}$.
Therefore, after $4s-1$ measurements as in R4 (with shifted harmonic samples)
there are $2^{s-1}$ solutions to R4. To resolve the ambiguity of which
$\bm{g}^{k}$ equals $\bm{g}^{+}$, we need an additional measurement - this is
the purpose of the arbitrary vector $\bm{a}$ in R5. W.p.1, only one of
$\bm{g}^{k}$, $0\leq k\leq 2^{s-1}-1$, will satisfy
$y_{m}^{+}=|\bm{a}^{T}\bm{V}(\bm{\theta}^{+})\bm{g}^{k}|^{2}$, resolving the
ambiguity and thus giving us $\bm{g}^{+}$.

In the derivations above, we assumed that $g_{k}\neq 0$, $0\leq k\leq s-1$;
i.e., we assumed exact knowledge of the number of non-zero elements. If
$g_{k}=0$ for some $k$, $\tilde{\bm{G}}(\bm{z})$ in (III-A) can have a null
space of more than one dimension. Assume that $S$, $S\leq s$, of the
$g_{0},\ldots,g_{s-1}$ are non-zero. One can now iterate through different
realizations of $\tilde{\bm{G}}(\bm{z})$, reducing the value of $s$ in each
step by 1 until a one-dimensional null space is encountered. Since for each
such realization of $\tilde{\bm{G}}(\bm{z})$, $\bm{w}^{+}$ in (75) belongs to
its null space, the solution will not be missed and the iteration terminates
after at most $s-S+1$ steps (when $s$ is reduced to $S$, giving the unique
solution $c\bm{w}^{+}$ as implied by the derivations above).

Thus, an algorithm for producing the $2^{s-1}$ solutions to R4 and that gives
the exponential recovery complexity in R5 is

R5-Alg-Harmonic:

  1. Input: $s,n\geq 4s-1,\bm{y}^{+}$, an arbitrary vector $\bm{a}$ and samples $\bm{z}\in\mathbb{T}^{m}$ where $4s-1\leq m<\min\\{n+1,8s-3\\}$.

  2. Assumptions: All samples in $\bm{z}$ are shifted harmonics.

  3. Output: $\bm{\theta}^{+},e^{i\alpha}\bm{g}^{+}$ for some unknown $\alpha$

  4. 1.

Construct the $m\times 4s$ matrix $\tilde{\bm{G}}(\bm{z})$ in (III-A).

  5. 2.

Keep reducing the value of $s$ (and thus the dimension of
$\tilde{\bm{G}}(\bm{z})$) by one until $\tilde{\bm{G}}(\bm{z})$ has a one-
dimensional null space (assume $s=S$ when this occurs for the fist time). Let
$\bm{w}$ denote a vector from its one-dimensional null space.

  6. 3.

Construct the Laurent polynomial $\hat{L}(z)=\sum_{k=0}^{2S}w_{k}z^{k-S}$. The
conjugate of the unique roots of $\hat{L}(z)$ equal $\bm{\theta}^{+}$.
$|\bm{g}^{+}|$ can now be computed (up to a global real scalar) from
Proposition 5.

  7. 4.

Construct the Laurent polynomial $Q(z)=\sum_{k=0}^{2S-2}w_{k+2S+1}z^{k-S+1}$.

  8. 5.

The $2S-2$ roots of $Q(z)$ can be paired as
$(r_{0},1/\overline{r_{0}}),\ldots,(r_{S-2},1/\overline{r_{S-2}})$. From each
pair, choose one element. This choice can be done in $2^{S-1}$ different ways
and let $\bm{q}^{k}$, $0\leq k\leq 2^{S-1}-1$, denote the $k$:th choice.

  9. 6.

For each $\bm{q}^{k}$, solve (82) with $s=S$ and let $\bm{g}^{k}$ denote the
solution (up to a complex scalar). Normalize each $\bm{g}^{k}$ so that
$\bm{y}^{+}=|\bm{V}^{T}(\bm{z})\bm{V}(\bm{\theta}^{+})\bm{g}^{k}|^{2}$

  10. 7.

Let
$\bm{g}^{+}=\arg\min_{\bm{g}\in\\{\bm{g}^{0},\ldots,\bm{g}^{2^{S-1}-1}\\}}|y_{m}-|\bm{a}^{T}\bm{V}(\bm{\theta}^{+})\bm{g}|^{2}|$.
Output $\bm{g}^{+}$ and $\bm{\theta}^{+}$.

From the theoretical results above, it follows that R5-Alg-Harmonic recovers
$\bm{\theta}^{+}$ and $\bm{g}^{+}$ w.p.1. Steps 1-6 produce all solutions to
R4 while step 7 finds the correct $\bm{g}^{+}$ and thereby gives the result in
R5. The complexity of R5-Alg-Harmonic is $O(\mathrm{poly}(s)2^{s})$. Up until
step 6, the total complexity is $O(\mathrm{poly}(s))$. Solving (82) has
complexity $O(\mathrm{poly}(s))$ but it is repeated $2^{s-1}$ times in step 6.
It is possible to utilize Theorem 6 to directly generate the dual solution for
each $\bm{g}^{k}$ and thereby reduce the number of times (82) is solved by
half, but this still results in complexity $O(\mathrm{poly}(s)2^{s})$.
Fortunately, it turns out that by taking $4s-2$ additional measurements, the
number of possible solutions can be reduced from $2^{s-1}$ to only two.

######  Theorem 8.

Given the assumptions in Theorem 7, further assume that
$z_{j}^{n}=e^{i\omega}$, $\hat{L}^{+}(z_{j})\neq 0$ for $4s-1\leq j\leq 6s-3$
and some $\omega$ satisfying $e^{i\gamma}\neq e^{i\omega}$ (with $\gamma$ from
Theorem 7); $z_{j}^{n}=e^{i\phi}$, $\hat{L}^{+}(z_{j})\neq 0$ for $6s-2\leq
j\leq 8s-4$ and some $\phi$ satisfying $e^{i\phi}\neq e^{i\gamma}$ and
$(e^{i\phi}-e^{i\gamma})(e^{-i\omega}-e^{-i\gamma})\neq(e^{-i\phi}-e^{-i\gamma})$.
Solving (52) (or equivalently (III-A)) results in
$\hat{L}(z)=c\hat{L}^{+}(z),L(z)=cL^{+}(z),\tilde{L}(z)=c\tilde{L}^{+}(z)$ for
some real scalar $c$.

###### Proof.

After the $4s-1$ measurements $z_{j}$, $0\leq j\leq 4s-2$, that satisfy
$z_{j}^{n}=e^{i\gamma}$ we conclude from Theorem 7 that
$\hat{L}(z)=c\hat{L}^{+}(z)$. Inserting this relation into (81) and dividing
out $\hat{L}^{+}(z)$, we get

$\displaystyle
cL^{+}(z)-L(z)=e^{i\gamma}(\tilde{L}(z)-c\tilde{L}^{+}(z))+e^{-i\gamma}(\overline{\tilde{L}(z)}-c\overline{\tilde{L}^{+}(z)}).$
(83)

Further using $\hat{L}(z)=c\hat{L}^{+}(z)$ in (79) and dividing out
$\hat{L}^{+}(z_{j})$ gives that

$\displaystyle
cL^{+}(z_{j})-L(z_{j})=z_{j}^{n}(\tilde{L}(z_{j})-c\tilde{L}^{+}(z_{j}))+z_{j}^{-n}(\overline{\tilde{L}(z_{j})}-c\overline{\tilde{L}^{+}(z_{j})})$
(84)

holds for all $z_{j}$ for which $\hat{L}^{+}(z_{j})\neq 0$. Subtracting (83),
evaluated at $z=z_{j}$, from (84) we get

$\displaystyle(z_{j}^{n}-e^{i\gamma})(\tilde{L}(z_{j})-c\tilde{L}^{+}(z_{j}))+(z_{j}^{-n}-e^{-i\gamma})(\overline{\tilde{L}(z_{j})}-c\overline{\tilde{L}^{+}(z_{j})})$
$\displaystyle=0.$ (85)

Proposition 3 shows that (85) is equivalent to a polynomial equation in
$z_{j}$ of degree $2s-2$. Hence, with the $2s-1$ samples $z_{j}$, $4s-1\leq
j\leq 6s-3$, for which $z_{j}^{n}=e^{i\omega}$, (85) becomes

$\displaystyle(e^{i\omega}-e^{i\gamma})(\tilde{L}(z_{j})-c\tilde{L}^{+}(z_{j}))+(e^{-i\omega}-e^{-i\gamma})(\overline{\tilde{L}(z_{j})}-c\overline{\tilde{L}^{+}(z_{j})})$
$\displaystyle=0,\quad 4s-1\leq j\leq 6s-3.$ (86)

Applying the identity theorem for polynomials to the equations in (86) implies
that

$\displaystyle(e^{i\omega}-e^{i\gamma})(\tilde{L}(z)-c\tilde{L}^{+}(z))+(e^{-i\omega}-e^{-i\gamma})(\overline{\tilde{L}(z)}-c\overline{\tilde{L}^{+}(z)})$
$\displaystyle=0,\quad\forall z$ (87)

Using the $2s-1$ samples $z_{j}=e^{i2\pi j/n}e^{i\phi}$, $6s-2\leq j\leq
8s-4$, in (85), we obtain

$\displaystyle(e^{i\phi}-e^{i\gamma})(\tilde{L}(z_{j})-c\tilde{L}^{+}(z_{j}))+(e^{-i\phi}-e^{-i\gamma})(\overline{\tilde{L}(z_{j})}-c\overline{\tilde{L}^{+}(z_{j})})$
$\displaystyle=0,\quad 6s-2\leq j\leq 8s-4.$ (88)

Applying again the identity theorem for polynomials, the equations in (88)
imply that

$\displaystyle(e^{i\phi}-e^{i\gamma})(\tilde{L}(z)-c\tilde{L}^{+}(z))+(e^{-i\phi}-e^{-i\gamma})(\overline{\tilde{L}(z)}-c\overline{\tilde{L}^{+}(z)})$
$\displaystyle=0,\quad\forall z.$ (89)

Multiplying (87) with
$\frac{(e^{-i\phi}-e^{-i\gamma})}{(e^{-i\omega}-e^{-i\gamma})}$, which is a
finite number due to the assumptions on $\phi,\omega,\gamma$, and subtracting
it from (89) we get

$\displaystyle\left((e^{i\phi}-e^{i\gamma})-\frac{(e^{-i\phi}-e^{-i\gamma})}{(e^{-i\omega}-e^{-i\gamma})}\right)(\tilde{L}(z)-c\tilde{L}^{+}(z))=0,\quad\forall
z.$ (90)

The assumptions on $\phi,\omega,\gamma$ imply that the scalar in front of
$\tilde{L}(z)-c\tilde{L}^{+}(z)$ in (90) is non-zero. Thus, the identity
theorem for polynomials implies that (90) given
$\tilde{L}(z)=c\tilde{L}^{+}(z)$. Inserting this relation into (83) we also
obtain $L(z)=cL^{+}(z)$, finishing the proof. ∎

Clearly, Theorem 8 holds even if $m\geq 8s-3$ as long as a subset of $8s-3$
measurements satisfies the assumptions of Theorem 8. Arbitrary
$\bm{\theta}^{+}\in\mathbb{T}^{s}$ and $\bm{g}^{+}$ satisfy the assumptions of
Theorem 8 w.p.1, so that $\bm{G}(\bm{z})$ in (III-A) has a one-dimensional
null space w.p.1. Moreover, for any $\bm{\theta}^{+}\in\mathbb{T}^{s}$ and
$\bm{g}^{+}$ with no zero elements, it is evidently easy to choose a $\bm{z}$
satisfying the assumptions of Theorem 8. Similarly to Theorem 5, we can relax
the requirements on $\bm{z}$ from Theorem 8.

######  Theorem 9.

Assume that $n\geq 4s-1$. Given any $\bm{\theta}^{+}\in\mathbb{T}^{s}$ and any
$\bm{g}^{+}\in\mathbb{C}^{s}$ with $s$ non-zero elements, solving (52) (or
equivalently (III-A)) for an arbitrary $\bm{z}\in\mathbb{T}^{8s-3}$ results in
$\hat{L}(z)=c\hat{L}^{+}(z),L(z)=cL^{+}(z),\tilde{L}(z)=c\tilde{L}^{+}(z)$
(for some real scalar $c$) w.p.1.

###### Proof.

For the given $\bm{\theta}^{+}$, $\bm{g}^{+}$, there clearly exists a vector
$\hat{\bm{z}}=[\hat{z}_{0},\ldots,\hat{z}_{8s-4}]$ of $8s-3$ samples that
satisfies the assumptions of Theorem 8. Thus, $\bm{G}(\hat{\bm{z}})$ in
(III-A) has a one-dimensional null space. Since $\bm{G}(\hat{\bm{z}})$ is a
$(8s-3)\times(8s-2)$ matrix, a one-dimensional null space of
$\bm{G}(\hat{\bm{z}})$ implies that there is a $(8s-3)\times(8s-3)$ sub-matrix
$\bm{M}(\hat{\bm{z}})$ of $\bm{G}(\hat{\bm{z}})$ with full rank $8s-3$.
$\bm{M}(\hat{\bm{z}})$ is obtained from $\bm{G}(\hat{\bm{z}})$ by removing one
column from $\bm{G}(\hat{\bm{z}})$. We will show that we can replace
$\hat{\bm{z}}$ with the arbitrary $\bm{z}$ and obtain that $\bm{M}(\bm{z})$ is
full rank w.p.1.

From (51), it follows that $y_{j}^{+}$ is a ratio of two Laurent polynomials
in the variable $z_{j}$. Thus, the determinant $\left|\bm{M}(\bm{z})\right|$
is a ratio $P(\bm{z})/Q(\bm{z})$ of multivariate Laurent polynomials
$P(\bm{z})$ and $Q(\bm{z})$. Since $\bm{M}(\hat{\bm{z}})$ is invertible,
$P(\hat{\bm{z}})\neq 0$, so $P(\bm{z})$ is not zero. The set
$\mathbb{T}^{8s-3}$ is homeomorphic to the $(8s-3)$-dimensional cube
$\mathbb{U}^{8s-3}$, where $\mathbb{U}=[-1,1)$ is the unit interval. Hence,
the zero set $\mathcal{S}=\\{\bm{z}|\bm{z}\in\mathbb{T}^{8s-3},P(\bm{z})=0\\}$
is either empty or is a set of points in $\mathbb{U}^{8s-3}$ taken from a
hypersurface of dimension at most $8s-4$. An arbitrary
$\bm{z}\in\mathbb{T}^{8s-3}$, corresponding to an arbitrary point in
$\mathbb{U}^{8s-3}$, is inside $\mathcal{S}$ w.p.0, proving the theorem. ∎

Theorem 9 is true for arbitrary $\bm{z}\in\mathbb{T}^{m}$ with $m\geq 8s-3$ as
well, since any subset of $8s-3$ samples from $\bm{z}$ ($8s-3$ rows from
$\bm{G}(\bm{z})$) satisfies Theorem 9 w.p.1. Similar to Proposition 5 we have

######  Proposition 6.

Given the assumptions in Theorem 9, $|\bm{g}^{+}|$ can be found, up to a real
scalar $c$, with $O(\mathrm{poly}(s))$ complexity.

###### Proof.

From Theorem 9 we obtain $\bm{\theta}^{+}$ and $L(z)=cL^{+}(z)$. Using
$(\theta_{k}^{+})^{-1}=\overline{\theta_{k}^{+}}$ and

$L(\overline{\theta_{k}^{+}})=cL^{+}(\overline{\theta_{k}^{+}})=2c|g_{k}^{+}|^{2}\left|t_{k}^{+}\left(\overline{\theta_{k}^{+}}\right)\right|^{2}$

it follows that

$c|g_{k}^{+}|^{2}=\frac{L(\overline{\theta_{k}^{+}})}{2\left|t_{k}^{+}\left(\overline{\theta_{k}^{+}}\right)\right|^{2}}.$

Since the operations leading to $c|\bm{g}^{+}|$ are of $O(\mathrm{poly}(s))$
complexity, it follows that $c|\bm{g}^{+}|$ can be found with
$O(\mathrm{poly}(s))$ complexity. ∎

Next, we show that given $L(z),\tilde{L}(z)$ in (III-A), we can solve for
$|\tilde{u}(z)|^{2},|\hat{u}(z)|^{2}$.

######  Lemma 4.

Given $L(z),\tilde{L}(z)$ in (III-A), we have that

$\displaystyle|\hat{u}(z)|^{2}=\frac{L(z)\pm\sqrt{L^{2}(z)-4K(z)}}{2}$
$\displaystyle|\tilde{u}(z)|^{2}=\frac{L(z)\mp\sqrt{L^{2}(z)-4K(z)}}{2}$

where $K(z)=|\tilde{L}(z)|^{2}$.

###### Proof.

$L(z)=|\tilde{u}(z)|^{2}+|\hat{u}(z)|^{2}$ implies
$|\tilde{u}(z)|^{2}=L(z)-|\hat{u}(z)|^{2}$. Furthermore,
$\tilde{L}(z)=\hat{u}(z)\overline{\tilde{u}(z)}$ gives
$K(z)=|\tilde{L}(z)|^{2}=|\hat{u}(z)|^{2}|\tilde{u}(z)|^{2}$. Hence,
$K(z)=(L(z)-|\hat{u}(z)|^{2})|\hat{u}(z)|^{2}=L(z)|\hat{u}(z)|^{2}-|\hat{u}(z)|^{4}$,
which is a quadratic equation in $|\hat{u}(z)|^{2}$. Solving the quadratic, we
obtain

$\displaystyle|\hat{u}(z)|^{2}=\frac{L(z)\pm\sqrt{L^{2}(z)-4K(z)}}{2}$ (91)

which together with $L(z)=|\tilde{u}(z)|^{2}+|\hat{u}(z)|^{2}$ gives

$\displaystyle|\tilde{u}(z)|^{2}=\frac{L(z)\mp\sqrt{L^{2}(z)-4K(z)}}{2}.$ (92)

∎

We have the option to choose sign in (92), resulting in the opposite sign in
(91). First, assume that $-$ is chosen in (92). Denote the roots of the
obtained $|\tilde{u}^{+}(z)|^{2}$, which appear as conjugate pairs, as
$(r_{0},1/\overline{r_{0}}),\ldots,(r_{s-2},1/\overline{r_{s-2}})$. For
arbitrary $\bm{g}^{+}$, since
$\tilde{L}^{+}(z)=\hat{u}^{+}(z)\overline{\tilde{u}^{+}(z)}$ and due to Lemma
2 and Theorem 1, $s-1$ roots of $|\tilde{u}^{+}(z)|^{2}$ appear as roots of
$\tilde{L}^{+}(z)$ w.p.1. Assume WLOG that the roots $r_{0},\ldots,r_{s-2}$
appear as roots of $\tilde{L}^{+}(z)$. Hence, we construct $\tilde{u}^{+}(z)$
(up to a complex scalar) as $\hat{u}^{+}(z)=\prod_{j=0}^{s-2}(z-r_{j})$ (the
other $s-1$ roots of $\tilde{L}^{+}(z)$ produce $\hat{u}^{+}(z)$). Solving the
equations

$\displaystyle\sum_{k=0}^{s-1}g_{k}t_{k}^{+}(r_{j})=0,\quad 0\leq j\leq s-2$
(93)

for $\bm{g}$ gives a unique $\bm{g}^{m}$ up to a complex scalar. Instead,
choosing $+$ as sign in (91) produces a different $|\tilde{u}^{+}(z)|^{2}$ and
a different solution $\bm{g}^{p}$ (up to a complex scalar) to the equations
(93). Since these are the only possible values for $\bm{g}$, one of them
equals $\bm{g}^{+}$ (up to a complex scalar). If no prior information about
$\bm{g}^{+}$ is at hand, then as before, in order to resolve the ambiguity, we
take an additional measurement $y_{m}$ with an arbitrary measurement vector.
Interestingly, Theorem 6 shows that the two possible solutions $\bm{g}^{m}$
and $\bm{g}^{p}$ are duals.

An important special case is when $L^{2}(z)=4K(z)$, which implies that
$|\hat{u}(z)|^{2}=|\tilde{u}(z)|^{2}=L(z)/2$. Theorem 1 shows that this occurs
for arbitrary $\bm{g}$ if and only if $\theta_{k}^{n}$, $0\leq k\leq s-1$, are
all equal. Since $\bm{\theta}\in\mathbb{T}^{s}$, this means that
$\theta_{k},0\leq k\leq s-1$, are shifted harmonics and thus $\bm{x}$ is
$s$-sparse in a DFT basis. In this case, the situation is the same as in
R5-Alg-Harmonic since there are $2^{s-1}$ different ways to choose the roots
for $\tilde{u}(z)$ (each one of them producing a solution to the equations
after normalization); thus, the recovery complexity suddenly becomes
exponential in $s$. Hence, the implication of this observation is that when
$\bm{x}$ is $s$-sparse in a DFT basis, the recovery complexity of $\bm{x}$ is
$O(\mathrm{poly}(s)2^{s})$ (as mentioned in 2.b in Section I). Nevertheless,
due to Theorem 1, as soon as not all $\theta_{j}^{n}$, $0\leq j\leq s-1$, are
equal (i.e., at least one frequency is not a shifted harmonic), $L^{2}(z)\neq
4K(z)$ w.p.1 since then the roots of $\tilde{u}(z)$ and $\hat{u}(z)$ are
different w.p.1.

If $S$ of the $g_{0},\ldots,g_{s-1}$ are non-zero, it is readily seen that
$\bm{w}^{+}$ in (63) belongs to the null space of $\bm{G}(\bm{z})$ for $s\geq
S$; thus the same iterative procedure as in step 2 of R5-Alg-Harmonic can be
used to find $S$. Hence, we have the following algorithm that recovers
$\bm{g}^{+},\bm{\theta}^{+}$

R5-Alg-General:

  1. Input: $s,n\geq 4s-1,\bm{y}^{+}$, an arbitrary vector $\bm{a}$ and arbitrary $\bm{z}\in\mathbb{T}^{m}$ where $m\geq 8s-3$.

  2. Output: $\bm{\theta}^{+},e^{i\alpha}\bm{g}^{+}$ for some unknown $\alpha$

  3. 1.

Construct the $m\times(8s-1)$ matrix $\bm{G}(\bm{z})$ in (III-A).

  4. 2.

Keep reducing the value of $s$ (and thus the dimension of $\bm{G}(\bm{z})$) by
one until $\bm{G}(\bm{z})$ has a one-dimensional null space (assume $s=S$ when
this occurs for the first time). Let $\bm{w}$ denote a vector from its one-
dimensional null space.

  5. 3.

Construct the Laurent polynomials
$\hat{L}(z)=\sum_{k=0}^{2S}w_{k}z^{k-S},\tilde{L}(z)=\sum_{k=0}^{2S-2}w_{k+2S+1}z^{k-S+1},L(z)=\sum_{k=0}^{2S-2}w_{k+4S}z^{k-S+1}$,
respectively. The conjugate of the unique roots of $\hat{L}(z)$ equal
$\bm{\theta}^{+}$. $|\bm{g}^{+}|$ can now be computed (up to a global real
scalar) from Proposition 6.

  6. 4.

Let $K(z)=|\tilde{L}(z)|^{2}$ and

$Q(z)=\frac{L(z)+\sqrt{L^{2}(z)-4K(z)}}{2}.$

We have two cases

    1. a.

$L^{2}(z)\neq 4K(z)$ ($(\theta_{j}^{+})^{n}$, $0\leq j\leq S-1$, are not all
equal):

       * •

Denote by $q_{0},\ldots,q_{S-2}$ the $S-1$ roots of $Q(z)$ that appear as
roots of $\tilde{L}(z)$ and solve for $\bm{g}$ the $S-1$ polynomial equations

$\sum_{k=0}^{S-1}g_{k}t_{k}(q_{j})=0,\quad 0\leq j\leq S-2.$

Denote by $\bm{g}^{a}$ the obtained solution which is unique up to a complex
scalar and construct its dual $\bm{g}^{b}$ as shown in Theorem 6. Let
$\mathcal{G}=\\{\bm{g}^{a},\bm{g}^{b}\\}$. Normalize each
$\bm{g}\in\mathcal{G}$ so that
$\bm{y}^{+}=|\bm{V}^{T}(\bm{z})\bm{V}(\bm{\theta}^{+})\bm{g}|^{2}$.

    2. b.

$L^{2}(z)=4K(z)$ ($(\theta_{j}^{+})^{n}$, $0\leq j\leq S-1$, are all equal):

       * •

In this case, $L(z)/2=|\hat{u}(z)|^{2}=|\tilde{u}(z)|^{2}$. The roots of
$L(z)/2$ can be paired as
$(l_{0},1/\overline{l_{0}}),\ldots,(l_{S-2},1/\overline{l_{S-2}})$. From each
pair, choose one element. This choice can be done in $2^{S-1}$ different ways
and let $\bm{l}^{k}$, $0\leq k\leq 2^{S-1}-1$, denote the $k$:th choice.

       * •

For each $\bm{l}^{k}$, $0\leq k\leq 2^{S-1}-1$, solve for $\bm{g}$ the $S-1$
polynomial equations

$\sum_{k=0}^{S-1}g_{k}t_{k}(l_{j})=0,\quad 0\leq j\leq S-2.$

Denote by $\bm{g}^{k}$ the obtained solution which is unique up to a complex
scalar. Let $\mathcal{G}=\\{\bm{g}^{0},\ldots,\bm{g}^{2^{S-1}-1}\\}$.
Normalize each $\bm{g}\in\mathcal{G}$ so that
$\bm{y}^{+}=|\bm{V}^{T}(\bm{z})\bm{V}(\bm{\theta}^{+})\bm{g}|^{2}$.

  7. 5.

Let
$\bm{g}^{+}=\arg\min_{\bm{g}\in\mathcal{G}}|y_{m}-|\bm{a}^{T}\bm{V}(\bm{\theta}^{+})\bm{g}|^{2}$.
Output $\bm{g}^{+}$ and $\bm{\theta}^{+}$.

###  III-B Derivation of R3

From the derivation of R4 and R5, we can easily construct an $\bm{A}$ and a
corresponding recovery algorithm that produces R3. Choose
$\bm{z}\in\mathbb{T}^{m}$ and $\bm{\theta}\in\mathbb{T}^{n}$ such that
$\bm{z}$ and any subset of $s$ samples from $\bm{\theta}$ satisfy the
assumptions of Theorem 7 or Theorem 8; it is clear from the assumptions of
these theorems that such $\bm{z}$ and $\bm{\theta}$ are easily found. Let
$\bm{a}\in\mathbb{C}^{n}$ be an arbitrary vector. Construct $\bm{A}$ in R3 as
$\bm{A}=[\bm{V}(\bm{z})^{T}\bm{V}(\bm{\theta});\bm{a}^{T}]$. Denote by
$p_{0},\ldots,p_{s-1}$ the unknown positions of the $s$ non-zero elements in
$\bm{x}$ and let $\bm{g}=[g_{0},\ldots,g_{s-1}]$ denote these elements. Hence,
$\bm{V}(\bm{z})^{T}\bm{V}(\bm{\theta})\bm{x}=\bm{V}(\bm{z})^{T}\bm{V}(\bm{\theta}_{\bm{p}})\bm{g}$,
where $\bm{\theta}_{\bm{p}}=[\theta_{p_{0}},\ldots,\theta_{p_{s-1}}]$.
Clearly, $\bm{\theta}_{\bm{p}}$ can be seen as the unknown $\bm{\theta}$ in
R4, meaning that the first $m$ measurements in R3 correspond to R4. From these
measurements, one obtains $\bm{\theta}_{\bm{p}}$ and the different $\bm{g}$
satisfying (45). The additional measurement with $\bm{a}$ enables application
of R5-Alg-Harmonic or R5-Alg-General which finally produce the results in R3.
Thus, a recovery algorithm for R3 is

R3-Alg:

  1. Input: $s,n\geq 4s-1,\bm{y}^{+}$, an arbitrary vector $\bm{a}$ and $\bm{A}=[\bm{V}(\bm{z})^{T}\bm{V}(\bm{\theta});\bm{a}^{T}]$. Sample vector $\bm{z}\in\mathbb{T}^{m}$, where $m\geq 4s-1$, and $\bm{\theta}\in\mathbb{T}^{n}$.

  2. Assumptions: If $4s-1\leq m<\min\\{n+1,8s-3\\}$, then all samples in $\bm{z}$ are shifted harmonics and $\bm{\theta}$ is chosen such that $\bm{z}$ and any subset of $s$ samples from $\bm{\theta}$ satisfy the assumptions of Theorem 7.
If $m\geq 8s-3$, $\bm{z}$ and $\bm{\theta}$ are chosen such that $\bm{z}$ and
any subset of $s$ samples from $\bm{\theta}$ satisfy the assumptions of
Theorem 8.

  3. Output: $\bm{x}$ (up to a global phase)

  4. 1.

If $4s-1\leq m\leq\min\\{n+1,8s-3\\}$, run algorithm R5-Alg-Harmonic until
step 3 from which we obtain $\bm{\theta}^{+}$. We have that

$p_{j}=\arg\min_{k}|\theta_{j}^{+}-\theta_{k}|,\quad 0\leq j\leq S-1$

are the positions of the $S$ non-zero elements in $\bm{x}$. Let
$\bm{g}=[x_{p_{0}},\ldots,x_{p_{S-1}}]^{T}$ be the vector of the unknown non-
zero elements from $\bm{x}$. Continue running R5-Alg-Harmonic from step 3 to
obtain $\bm{g}$ as $\bm{g}=\bm{g}^{+}$ in the end.

  5. 2.

If $m\geq 8s-3$, run algorithm R5-Alg-General until step 3 from which we
obtain $\bm{\theta}^{+}$. We have that

$p_{j}=\arg\min_{k}|\theta_{j}^{+}-\theta_{k}|,\quad 0\leq j\leq S-1$

are the positions of the $S$ non-zero elements in $\bm{x}$. Let
$\bm{g}=[x_{p_{0}},\ldots,x_{p_{S-1}}]^{T}$ be the vector of the unknown non-
zero elements from $\bm{x}$. Continue running R5-Alg-General from step 3 to
obtain $\bm{g}$ as $\bm{g}=\bm{g}^{+}$ in the end.

R3-Alg shows that there is a unique $\bm{x}$ (up to a global phase) satisfying
the constraint in (4) for
$\bm{A}=[\bm{V}(\bm{z})^{T}\bm{V}(\bm{\theta});\bm{a}^{T}]$ with $m=O(s)$.
This observation is what solves, and significantly improves upon, open problem
3.5 in [18].

##  IV Acknowledgements

I thank Dr. Zhibin Yu for proofreading the work.

## References

  * [1] E. J. Candes, J. Romberg, T. Tao, “Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information,” _IEEE Trans. Inf. Theory_ , vol. 52, no. 2, pp. 489–509, Feb. 2006.
  * [2] D. L. Donoho, “Compressed sensing,” _IEEE Trans. Inf. Theory_ , vol. 52, no. 4, pp. 1289–1306, Apr. 2006.
  * [3] E. J. Candes, T. Tao, “Near-optimal signal recovery from random projections: Universal encoding strategies?” _IEEE Trans. Inf. Theory_ , vol. 52, no. 12, pp. 5406–5425, Dec. 2006.
  * [4] E. J. Candes, T. Tao, “Decoding by Linear Programming,” _IEEE Trans. Inf. Theory_ , vol. 51, no. 12, pp. 4203–4215, Dec. 2005.
  * [5] M. Rani, S. B. Dhok, R. B. Deshmukh, “A Systematic Review of Compressive Sensing: Concepts, Implementations and Applications,” _IEEE Access_ , vol. 6, pp. 4875 - 4894, Jan. 2018.
  * [6] M. Huang, Z. Xu, “No existence of linear algorithm for Fourier phase retrieval,” available on https://arxiv.org/abs/2209.05673
  * [7] J. Bourgain, S. Dilworth, K. Ford, S. Konyagin, D. Kutzarova, “Explicit constructions of RIP matrices and related problems,” Duke Math. J., vol. 159, no. 1, pp. 145–185, 2011.
  * [8] Y.M. Lu, M. Vetterli, “Sparse spectral factorization: Unicity and reconstruction algorithms,” ICASSP 2011.
  * [9] Y. Shechtman, A. Beck, Y. C. Eldar, “GESPAR: Efficient Phase Retrieval of Sparse Signals,” IEEE Trans. Sig. Proc., vol. 64, no. 4, pp. 928-938, Feb. 2014.
  * [10] K. Jaganathan, S. Oymak, B. Hassibi, “Sparse phase retrieval: Convex algorithms and limitations,” IEEE Int. Symp. Inf. Theory, Istanbul, Turkey, Jul. 2013.
  * [11] M. Moravec, J. Romberg, R. Baraniuk, “Compressive phase retrieval,” in Proc. SPIE 6701, Wavelets XII, 670120 (27 September 2007).
  * [12] Y. Xia, Z. Xu, “Sparse Phase Retrieval Via PhaseLiftOff,” IEEE Trasn. Sig. Proc., vol. 64, no. 4, pp. 928-938, Feb. 2014.
  * [13] R. Pedarsani, D. Yin, K. Lee, K. Ramchandran, “PhaseCode: Fast and Efficient Compressive Phase Retrieval Based on Sparse-Graph Codes,” _IEEE Trans. Inf. Theory_ , vol. 63, no. 6, pp. 3663 - 3691, Apr. 2017.
  * [14] Roger A. Horn, Charles R. Johnson (1991), Topics in matrix analysis, Cambridge University Press.
  * [15] A. Conca, D. Edidin, M. Hering, C. Vinzant, “An algebraic characterization of injectivity in phase retrieval,” Appl. Comput. Harmon. Anal., vol. 38, no. 2, pp. 346–356, 2015.
  * [16] Y. Xia, Z. Xu, “The recovery of complex sparse signals from few phaseless measurements,” available on https://arxiv.org/abs/1911.11301
  * [17] J. Dong, L. Valzania, A. Maillard, T. Pham, S. Gigan, M. Unser, “Phase Retrieval: From Computational Imaging to Machine Learning: A tutorial,” IEEE Sig. Proc. Mag., vol. 40, no. 1, pp. 45 - 57, Jan. 2023.
  * [18] Z. Xu, “The Minimal Measurement Number Problem in Phase Retrieval: A Review of Recent Developments,” _Journal of Mathematical Research with Applications_ , vol. 37, no. 1, pp. 40–46, Jan. 2017.
  * [19] G. Tang, B. N. Bhaskar, P. Shah, B. Recht, “Compressed Sensing Off the Grid,” _IEEE Trans. Inf. Theory_ , vol. 59, no. 11, Nov. 2013.
  * [20] D. Kapetanović, “Sparse Channel Estimation From Discrete-Time Fourier Transform Beam Measurements,” IEEE Trans. Wire. Comm., vol. 22, no. 9, pp. 6356 - 6368, Sept. 2023.