Phase retrieval and the topology of zero sets

Introduction

In optics, detectors measure intensity: they record the squared modulus $|E|^2$ of an electromagnetic field $E$, but not its phase. This is a fundamental limitation. The intensity tells you how much light arrives at each point, but the phase encodes the relative timing of the wave, and with it, much of the structural information about the object that produced the field. Recovering this lost phase from intensity measurements alone is known as the phase retrieval problem.

In many optical setups, the relationship between the object and the measured field is well described by Fourier optics: the field in the far field or at the focal plane of a lens is the Fourier transform of the field at the object plane. This means that the problem has a clean mathematical formulation. Suppose you have a compactly supported function $f(k)$ and you form its Fourier transform $F(x) = \int f(k) \, e^{-ikx} \, dk$. You observe $|F(x)|^2$ for all real $x$, and you want to recover $\operatorname{arg} F(x)$. The question is: does the intensity uniquely determine the phase, and if so, how do you compute it?

The approach taken in this post is based on the zeros of the complex-analytic extension of $F$. Since $F(z)$ is entire, it is determined — up to an exponential factor — by its zeros. The phase of $F(x)$ on the real line is then expressible entirely in terms of these zeros, and the phase retrieval problem reduces to the problem of locating them.

We begin with a review of the Weierstrass and Hadamard factorisation theorems, which describe how entire functions decompose into products over their zeros. We then specialise to Fourier transforms of compactly supported functions, derive a closed-form expression for the derivative of the phase in terms of zeros, and show that in 1D, the phase cannot be uniquely recovered from intensity alone. In 2D, however, the situation changes: the zero set acquires a topological structure that resolves the ambiguity.

Fundamental theorem of algebra, stretched to its limits

It is a commonly accepted and embraced miracle that complex polynomials are uniquely determined by their zeros (counted with multiplicity). That is, any polynomial function $f(z) \in \mathbb{C}[z]$ uniquely factors as

where $C$ is a constant and $\{a_n\}$ are the zeros of $f(z)$.

What about general complex-analytic functions? Such functions can have infinitely many zeros. For instance, the zeros of the function $f(z) = \sin \pi z$ occur exactly at the integers, suggesting a factorisation

Alas, this won’t converge — but it was the insight of Euler that, if we arrange the terms a bit, we get something convergent:

So perhaps the fundamental theorem of algebra can be extended as follows.

Fake Theorem 1. Let $f(z)$ be a complex-analytic function with zeros $\{a_n\}$, counted with multiplicity. Then $f(z)$ decomposes as

Unfortunately, this cannot be true either. In fact, our theorem is false for two reasons.

1. The expression fails to converge. For instance, it will fail if our zeros are $\{1, 2, \ldots\}$, not to mention if we have a zero at $0$.
2. There are entire functions without zeros. In fact, such functions are easy to produce: for any entire function $g(z)$, the function $e^{g(z)}$ is entire and never vanishes.

One man’s bug is another man’s feature, however. Since $e^{g(z)}$ has no zeros, that means that we can multiply any factor $(1 - w)$ for $w = z / a_n$ by an exponential without perturbing the zero set. Can we use that freedom to drag each multiplicative factor in our product close to $1$, so that the product converges?

For $p \geq 0$, we define the Weierstrass elementary factor $E_p(w)$ as

Each $E_p$ is entire with a single simple zero at $w = 1$ and no others. The point of the construction is visible the moment you write it in exponential form:

We’ve surgically removed the first $p$ terms of the log series, and $E_p(w)$ now hugs $1$ to order $p+1$:

or more precisely

Returning to a general entire function $f(z)$, let ${a_n}$ be its nonzero zeros, ordered so that $|a_1| \leq |a_2| \leq \cdots$. Actually, let’s take $\{a_n\}$ to be the nonzero zeros, and if $f(0) = 0$, let’s let $m$ be the order of that zero, since we’ll treat it separately. Consider the product

Fix a disk $|z| \leq R$. For $n$ large enough so that $|a_n| > R$, we have $|z / a_n| \leq 1$, so we know that

Now choose the exponents $p_n$ large enough so that

for every $R$. If anything, you can just pick $p_n = n$. Don’t worry about what $p_n$ to pick just yet. Just know that it can be done.

Now let’s consider the ratio $h(z) = f(z) / P(z)$. By construction, every zero cancels, with matching orders, so $h(z)$ is entire and never zero. Now we already knew such functions existed: it was our second point of failure of Fake Theorem 1: for any entire function $g(z)$, the exponential $e^{g(z)}$ is entire without zeros. As a matter of fact, this characterises them:

Lemma 2. If $h(z)$ is an entire function without zeros, then it is the exponential $e^{g(z)}$ of another entire function $g(z)$.

It’s easy to construct $g(z)$: it’s the antiderivative of $h'(z) / h(z)$. Actually, there is some nontriviality in the existence of this antiderivative, which hinges on $\mathbb{C}$ being simply connected. For instance, $h(z) = z$ is entire and without zeros on $\mathbb{C} \setminus \{0\}$, but $h' / h = 1/z$ has no primitive there because

In any case, going back to our construction, the ratio $h(z) = f(z) / P(z)$ is entire and without zeros, so by Lemma 2, it’s the exponential of another entire function. This brings us to:

Theorem 3. Let $f(z)$ be an entire function, and let $\{a_n\}$ be the nonzero zeros of $f$, counted with multiplicity. If $f(z)$ has a zero at $z = 0$, let $m$ be the order of that zero. Then there exists an entire function $g$ and a sequence of integers $\{p_n\}$ such that

Example 4. Turning back to $\sin \pi z$, we set $m = 1$ and $e^{g(z)} = \pi$, and write

Notice what happens when we pair up the positive and negative integers:

Thus the Weierstrass product simplifies to

Taming the factorisation

The Weierstrass factorisation theorem gives us a factorisation in great generality; but vagueness is the price we pay for that generality. Indeed, notice how loose the factorisation is:

Bad News 5. The exponents $p_n$ can be whatever, so long as the product converges.

Bad News 6. $g(z)$ is any entire function — an infinite-dimensional unknown.

So a factorisation exists, but you’re left with an infinite amount of slack in how to write out the factorisation. What we’ll find out in this section is that the moment you bound how fast $f(z)$ grows, all of that slack collapses to finitely many numbers. Both $p_n$ and $g(z)$ become controlled by a single parameter.

For an entire function $f(z)$, define

and define the order of $f(z)$ to be

That’s a scary-looking expression, but we’re effectively just saying that $M(r)$ grows like $e^{r^{\rho}}$ and we want to know $\rho$.

• If $f(z)$ is a polynomial, $\rho = 0$.
• If $f(z) = e^z$, $\rho = 1$, and if $f(z) = e^{g(z)}$ for some polynomial of degree $d$, then $\rho = d$.
• The order can be infinite, for instance for $e^{e^z}$.

Order is just a measure of growth rate, and we’ll find that it controls both the density of zeros, and the nature of the exponential $g(z)$ in the Weierstrass factorisation, thereby addressing Bad News 5 and 6.

Let’s start with Bad News 5. We will be working inside the disk of radius $r$. For a point $a$ inside the disk, define the Blaschke factor $\varphi_a(z)$ as

This function has the following properties:

• It has one zero, at $a$, inside the disk.
• It has one pole, at $r^2 / \overline{a}$, outside the disk.
• On the boundary circle, the modulus is exactly $1$.

Now consider an analytic function $f(z)$. Assume for convenience that $f(0) \neq 0$. Strip off the zeros of $f(z)$ by setting

Then $h(z)$ is analytic and zero-free in the disk, and on the boundary circle, $|h(z)| = |f(z)|$. The maximum modulus principle implies that

Evaluating $h(0)$, we conclude that

What does this mean intuitively? If we add a zero to $f$ within a given disk, then $f$ has to swing correspondingly large elsewhere on that disk. Consequently, the order of $f$ is directly related to the density of zeros of $f$.

If $f$ has order $\rho$, then $\log M(r) \leq r^{\rho + \varepsilon}$, hence

Since every zero with $|a_n| < r / e$ contributes at least $\log e = 1$ to the sum, the zero count $n(r)$ grows at most as $O(r^{\rho + \varepsilon})$. Consequently,

This addresses Bad News 5: If $f(z)$ has order $\rho$, then a single fixed $p = \lfloor \rho \rfloor$ makes $\sum_n 1/|a_n|^{p + 1}$ converge, so one and the same elementary factor $E_p$ works for every zero.

Now let’s turn to Bad News 6. Let $f(z)$ be an entire function of order $\rho$. Then we know that we can write

where now the “$p$” in $E_p$ is fixed (and set to $\lfloor \rho \rfloor$). What can we say about $g(z)$? It’s easy to see that $e^{g(z)}$ must have order $\rho$, too. Hence

which implies that

So the real part of $g(z)$ grows at most polynomially. But if the real part grows polynomially, then so does the imaginary part. In fact, this is another instance of the maximum modulus principle. If we set $A = \sup_{|z| \leq r} \operatorname{Re} g(z)$, then $g(z)$ takes values in the half-plane $\mathcal{H}_A = \{w \in \mathbb{C} : \operatorname{Re} w < A\}$. A half-plane maps conformally to a unit disk by a Mobius transformation, and on this unit disk, the maximum modulus principle kicks in to bound the imaginary part.

Concretely, the Mobius transformation

maps $\mathcal{H}_A$ to the unit disk. Now assume without loss of generality that $g(0) = 0$ (otherwise take $g(z) - g(0)$) and consider the composition $\varphi(z) = T(g(z))$. It maps the disk of radius $r$ to the unit disk, and it has a single zero, at $z = 0$. Now look at $\varphi(z) / z$. It’s analytic, has no zeros, and it’s bounded by $1/r$ on $|z| = r$, so by the maximum modulus principle,

Now apply the inverted Mobius map,

By definition,

Take moduli, use the reverse triangle inequality to find that

Now pick $r = 2|z|$ to find

The full modulus on a circle is controlled by the real part on a circle twice as big. So if $\operatorname{Re} g(z)$ grows polynomially, then so do $\operatorname{Im} g(z)$ and $|g(z)|$. And since every complex analytic function has a Taylor series, we can conclude that $g(z)$ is a polynomial. This brings us to:

Theorem 7 (Hadamard Factorisation Theorem). Let $f(z)$ be an entire function of order $\rho$. If it has a zero at $0$, let $m$ be the order of that zero. Let $a_n$ be the other zeros, counted with multiplicity. Then

where $p = \lfloor \rho \rfloor$ and $g(z)$ is a polynomial of degree $\leq \rho$.

Not tame enough

We have seen that the main obstacle to a clean factorisation of an analytic function is its growth rate: if the function grows too fast, the desired products fail to converge.

Let’s take a look at the Hadamard factorisation at order $\rho = 1$:

Under what conditions can we get rid of the exponential factors? And can we determine $\alpha z + \beta$?

An entire function $f(z)$ is said to be of exponential type if the inequality

holds for some positive $A$ and $B$. A useful condition to impose on functions of exponential type is the convergence of the integral

Such functions are called functions of class $C$. The distribution of zeros is particularly well-behaved for such functions, and accordingly, the Hadamard factorisation can be further tamed.

Theorem 8. Let $f(z)$ be a function of class $C$. If it has a zero at $z = 0$, let $m$ be its order. Let $\{a_n\}$ be the remaining zeros, counted with multiplicity. Then $f(z)$ factorises as

where $\alpha \in \mathbb{R}$. The data $(C,m,\alpha,\{a_n\})$ uniquely determine $f(z)$, and the real number $\alpha$ is the asymmetry of vertical growth,

where

Back to Fourier transforms

Consider a function $f(k) \in L^2(\mathbb{R})$ with compact support $S$. Let $a = \operatorname{ess} \operatorname{inf} S$ and $b = \operatorname{ess} \operatorname{sup} S$. Form its Fourier transform

It is not hard to see that this function is of class $C$. Moreover, one can verify that

As such, we can factorise $F(z)$ as in Theorem 8, with $\alpha = (a + b) / 2$.

Theorem 9. Adapt the notation as above. Let $\{a_n\} = \{\xi_n + i \eta_n\}$ be the zeros of $F(z)$, counted with multiplicity. Then

for all $x \in \mathbb{R}$.

Notice that the individual summands are effectively just bump functions, which get progressively sharper when a zero approaches the real line.

To see why this is true, simply expand

and insert the factorisation of $F(x)$:

In words, then, we have a closed-form expression for the derivative of the phase of $F(x)$, depending entirely on its zeros.

Example 10. It is reasonably straightforward to verify the claim numerically. In the example below, I generated a random smooth complex spectrum $f(k)$ on the interval $[1, 3.5]$ consisting of a few random Gaussian bumps. The Fourier transform $F(z)$ is computed numerically, after which the zeros are detected and extracted by counting the winding number

The true phase slope is extracted directly from the computed Fourier transform, while the predicted phase slope is obtained by applying Theorem 9, using the detected zeros.

We can immediately observe a few things:

• We’re seeing a near-perfect match between the true and predicted phase slope. The deviations near the edges of the $x$-range are merely because there are additional zeros outside our search range.
• Most zeros are very close to the real axis. This is typical of functions of class $C$.

1D phase extraction

We now turn to the main phase retrieval problem in 1D. Suppose that you only know $|F(x)|^2$ for real $x$. To what extent can you recover the zeros of the complex-analytic extension $F(z)$, and thereby recover the derivative of the phase?

Since $|F(x)|^2 = F(x) \, \overline{F(x)}$, it’s really a matter of finding the analytic continuation $z \mapsto F(z) \, \overline{F(\overline{z})}$ that extends the known values of the real line. Here’s how we do that. First, note that

So $|F(x)|^2$ is the Fourier transform of the function

Conversely then, $g(\tau)$ must be the inverse Fourier transform of $|F(x)|^2$:

Notice that this is something we can explicitly compute with the data that we got. As such, we are also able to directly compute, for any $z \in \mathbb{C}$, the values of

But now simply note that

so we’ve got our expression for $F(z) \, \overline{F(\overline{z})}$.

Now, if $G(z) = 0$ at some complex number $z_0$, then obviously $G(z) = 0$ at $\overline{z_0}$ as well. This raises the question: which of the two is the zero of $F(z)$? Is it $z_0$ or $\overline{z_0}$? As it turns out, it is fundamentally impossible to tell. To see this, notice that, if $z_0$ is a zero of $F(z)$, then we can define

which has the exact same intensity on the real axis, but with a zero $\overline{z_0}$ instead of $z_0$.

What are the implications for the phase? One can verify that

which is a $2\pi$ step located at $x = \xi_0$, with $\eta_0$ (i.e. the distance of the zero to the real axis) governing the sharpness of the $2\pi$ step.

In fact, this agrees with the expression in Theorem 9. If we flip a zero, $z_0 \mapsto \overline{z_0}$, then the summand corresponding to that zero reverses sign:

which is precisely the derivative of the arctangent expression above.

Moving to 2D

In the previous section, we saw that we are fundamentally unable to determine the phase of $F(x)$ from the intensity alone, and we were able to pinpoint the exact source of ambiguity.

Moving to 2D, however, the situation dramatically changes. In fact, it is a known theorem of Hayes that the multidimensional phase retrieval problem has a unique solution, up to trivial ambiguities. Therefore, if we know $|F(x,y)|^2$ for real $x,y$, we should, in theory, be able to uniquely recover $\operatorname{arg} F(x,y)$.

Suppose that we know $|F(x,y)|^2$ on $\mathbb{R}^2$. Now set $y = 0$. By the projection–slice theorem, $F(x,0)$ is the (real) Fourier transform of the projection of $f(k_x,k_y)$ onto the first axis. We can now play the exact same game that we did in the previous section, and infer the zeros of $z \mapsto F(z,0)$ up to the aforementioned fundamental ambiguity of not being able to tell a zero from its conjugate.

In fact, we can apply this exact same game to the restriction of $|F(x,y)|^2$ to any line in $\mathbb{R}^2$. And every time we do so, we encounter the same fundamental ambiguity, that we cannot tell a zero from its conjugate.

Yet the theorem of Hayes tells us that there must be a way to resolve the ambiguity when all of the data is taken together. How do we do that? What additional information do we have that is lost when looking only at individual slices?

The key insight is that the zeros, viewed through any slice, form scattered points in $\mathbb{C}$, but when taken altogether, they form a well-defined, coherent complex curve in $\mathbb{C}^2$. Mathematically, the structure of $V = F^{-1}(0)$ for a holomorphic function $F \colon \mathbb{C}^2 \to \mathbb{C}$ is that of an analytic hypersurface. If $F$ is irreducible, then such a hypersurface is topologically connected.

To illustrate how this works, imagine that the zeros of $F(z,w) \colon \mathbb{C}^2 \to \mathbb{C}^2$ form a 1D curve as illustrated below.

Intersecting the curve with a line yields isolated points, and for each such point, we cannot distinguish it from its conjugate. But taken together, the zeros trace out a connected curve. The only remaining ambiguity is therefore a single global conjugation, rather than independent conjugations of individual points.

Attempt on simulated data

tbc...