Phase recovery using MaxCut

In this example, we relax the phase retrieval problem similar to the classical MaxCut semidefinite program and recover the phase of the signal given the magnitude of the linear measurements.

Phase recovery has wide applications such as in X-ray and crystallography imaging, diffraction imaging or microscopy and audio signal processing. In all these applications, the detectors cannot measure the phase of the incoming wave and only record its amplitude i.e complex measurements of a signal $x \in \mathbb{C}^p$ are obtained from a linear injective operator $A$, but we can only measure the magnitude vector $Ax$, not the phase of $Ax$.

Recovering the phase of $Ax$ from $|Ax|$ is a nonconvex optimization problem. Using results from this paper, the problem can be relaxed to a (complex) semidefinite program (complex SDP).

The original representation of the problem is as follows:

\[\begin{array}{ll} \text{find} & x \in \mathbb{C}^p \\ \text{subject to} & |Ax| = b \end{array}\]

where $A \in \mathbb{C}^{n \times p}$ and $b \in \mathbb{R}^n$.

In this example, the problem is to find the phase of $Ax$ given the value $|Ax|$. Given a linear operator $A$ and a vector $b= |Ax|$ of measured amplitudes, in the noiseless case, we can write $Ax = \text{diag}(b)u$ where $u \in \mathbb{C}^n$ is a phase vector, satisfying $|\mathbb{u}_i| = 1$ for $i = 1,\ldots, n$.

We relax this problem as Complex Semidefinite Programming.

Relaxed Problem similar to MaxCut

Define the positive semidefinite hermitian matrix $M = \text{diag}(b) (I - A A^*) \text{diag}(b)$. The problem is:

\[\begin{array}{ll} \text{minimize} & \langle U, M \rangle \\ \text{subject to} & \text{diag}(U) = 1\\ & U \succeq 0 \end{array}\]

Here the variable $U$ must be hermitian ($U \in \mathbb{H}_n $), and we have a solution to the phase recovery problem if $U = u u^*$ has rank one. Otherwise, the leading singular vector of $U$ can be used to approximate the solution.

using Convex, SCS, LinearAlgebra

n = 20
p = 2
A = rand(n, p) + im * randn(n, p)
x = rand(p) + im * randn(p)
b = abs.(A * x) + rand(n)

M = diagm(b) * (I(n) - A * A') * diagm(b)
U = ComplexVariable(n, n)
objective = inner_product(U, M)
c1 = diag(U) == 1
c2 = isposdef(U)
p = minimize(objective, c1, c2)
solve!(p, SCS.Optimizer; silent = true)
Problem statistics
  problem is DCP         : true
  number of variables    : 1 (800 scalar elements)
  number of constraints  : 2 (1_640 scalar elements)
  number of coefficients : 1_221
  number of atoms        : 20

Solution summary
  termination status : OPTIMAL
  primal status      : FEASIBLE_POINT
  dual status        : FEASIBLE_POINT
  objective value    : -2696.977

Expression graph
  minimize
   └─ real (affine; real)
      └─ sum (affine; complex)
         └─ diag (affine; complex)
            └─ …
  subject to
   ├─ == constraint (affine)
   │  └─ + (affine; complex)
   │     ├─ diag (affine; complex)
   │     │  └─ …
   │     └─ Convex.NegateAtom (constant; negative)
   │        └─ …
   └─ PSD constraint (convex)
      └─ vcat (affine; real)
         ├─ hcat (affine; real)
         │  ├─ …
         │  └─ …
         └─ hcat (affine; real)
            ├─ …
            └─ …
evaluate(U)
20×20 Matrix{ComplexF64}:
       1.0-1.72242e-18im  …    0.908392-0.4182im
 0.0843738+0.996467im          0.493334+0.869839im
 -0.701167+0.713019im         -0.338727+0.940865im
 -0.438539+0.898749im        -0.0225063+0.999746im
 -0.393766+0.919246im         0.0267331+0.999642im
 -0.355226+0.934814im     …   0.0682504+0.997666im
 -0.874913+0.484347im          -0.59217+0.805812im
  0.999311-0.0378387im         0.891881-0.452254im
  0.789284+0.614079im          0.973721+0.22773im
  0.356749+0.934236im          0.714717+0.699412im
 -0.642287+0.7665im       …    -0.26288+0.964822im
  0.377471+0.926055im          0.730119+0.683316im
 -0.952214+0.305534im         -0.737159+0.675715im
   0.99397+0.109946im          0.948829-0.315783im
 -0.568588+0.822661im         -0.172452+0.985016im
 -0.469565+0.882935im     …  -0.0573014+0.998355im
 -0.191165+0.981592im          0.236834+0.971549im
 -0.405244+0.914245im         0.0142165+0.999898im
  0.867785+0.497007im           0.99607+0.088563im
  0.908392+0.4182im                 1.0+1.00568e-17im

Verify if the rank of $U$ is 1:

B, C = eigen(evaluate(U));
length([e for e in B if (abs(real(e)) > 1e-4)])
1

Decompose $U = uu^*$ where $u$ is the phase of $Ax$

u = C[:, 1];
for i in 1:n
    u[i] = u[i] / abs(u[i])
end
u
20-element Vector{ComplexF64}:
                  1.0 + 0.0im
 -0.08587617156793281 - 0.9963058180884196im
   0.7007444287759066 - 0.7134123951400961im
   0.4389476384188856 - 0.8985126436096953im
  0.39554091208698205 - 0.9184483583007802im
   0.3591401372349993 - 0.9332836448941051im
   0.8732753348467414 - 0.4872270410684444im
  -0.9986936614787048 + 0.05109765671983609im
  -0.7826418178857034 - 0.6224723165704331im
 -0.35306666434319034 - 0.9355981672329061im
   0.6376163361964394 - 0.7703540795085265im
 -0.37182170163690464 - 0.9283041646959452im
   0.9512405083789296 - 0.30845015029822154im
  -0.9943767183319128 - 0.10590062341391507im
   0.5687889470923485 - 0.8224835157409404im
  0.47214939201546746 - 0.8815185486530753im
  0.18909979675222313 - 0.9819578742839572im
   0.4049479605176449 - 0.9143397340554549im
  -0.8654893733197792 - 0.5009272848134108im
  -0.9098300101257855 - 0.4149811473723996im

This page was generated using Literate.jl.