1

Low-complexity implementation of convex
optimization-based phase retrieval
Sercan Ö. Arık, and Joseph M. Kahn, Fellow, IEEE

Abstract— Phase retrieval has important applications in
optical imaging, communications and sensing. Lifting the
dimensionality of the problem allows phase retrieval to be
approximated as a convex optimization problem in a higherdimensional space. Convex optimization-based phase retrieval
has been shown to yield high accuracy, yet its low-complexity
implementation has not been explored. In this paper, we study
three fundamental approaches for its low-complexity
implementation: the projected gradient method, the Nesterov
accelerated gradient method, and the alternating direction
method of multipliers (ADMM). We derive the corresponding
estimation algorithms and evaluate their complexities. We
compare their performance in the application area of directdetection mode-division multiplexing. We demonstrate that they
yield small estimation penalties (less than 0.2 dB for transmitter
processing and less than 0.6 dB for receiver equalization) while
yielding low computational cost, as their implementation
complexities all scale quadratically in the number of unknown
parameters. Among the three methods, ADMM achieves
convergence after the fewest iterations and the fewest
computational operations.
Index Terms—Phase retrieval, convex optimization, alternative
direction method of multipliers, optical communications, modedivision multiplexing.

P

I. INTRODUCTION

HASE retrieval corresponds to the extraction of phase
information about signals solely from their measured
intensities, possibly with the aid of various transformations.
Phase retrieval problems are commonly encountered in
applications where phase measurements are impossible or
impractical. One such field is optics, due to the fundamental
difficulty of measuring phase at carrier frequencies of
hundreds of THz [1]. The applications of phase retrieval in
optics include imaging (crystallographic imaging [2-4],
speckle imaging [5], coherent diffractive imaging [6,7],
astronomical imaging [8,9]), sensing (wavefront sensing [11]),
and
communications
(direct-detection
mode-division
multiplexing, MDM [10]). There are also important
applications of phase retrieval outside of optics, such as multiinput multi-output (MIMO) wireless communications [12],
This manuscript was submitted in September 2017. This work was
supported in part by Huawei Technologies Co., Ltd.
S. Ö. Arık is with Baidu Silicon Valley Artificial Intelligence Lab,
Sunnyvale, CA 94089 (e-mail: sercanarik@baidu.com).
J. M. Kahn is with the E. L. Ginzton Laboratory, Department of Electrical
Engineering, Stanford University, Stanford, CA 94305 USA (e-mail:
jmk@ee.stanford.edu).

speech recognition [13] and compressed sensing [14-16].
The traditional setting of phase retrieval is in the recovery
of signals from measured intensities in the space/spatial
frequency or time/temporal frequency domains. This setting is
encountered in many applications, because many systems
(e.g., optical systems such as graded-index media, lenses, or
simply free-space propagation) transform signals from the
spatial domain into the spatial frequency domain, and many
detection systems measure only intensities. The approaches
typically used for phase retrieval in this traditional problem
setting are of the form of sequential gradient descent. For
example, the Gerchberg-Saxton algorithm is based on
iteratively transforming a candidate solution back and forth
between the two domains (e.g., space and spatial frequency),
imposing constraints in one domain before transforming to the
other [17,18]. Despite the good performance of such heuristic
approaches, there is no algorithm for phase retrieval that
guarantees achievement of a globally optimal solution at a low
(non-deterministic polynomial-time hard) computational
complexity. A fundamental challenge is the non-convexity of
the intensity constraint after an orthogonal transformation
[19].
Furthermore, in many applications, the problem setting for
phase retrieval is more general than that addressed by
traditional algorithms. An important generalization is
encountered in recovery of unknown signals from the
measured intensities of their inner products with known
signals. (When the known signals are complex exponentials,
we obtain the special case described in the previous
paragraph.) The inner product operation occurs in linear
systems when measurements are performed in the appropriate
domain. For example, for a linear MIMO optical channel at a
given frequency, the known signals may correspond to
(baseband envelopes of) known training sequence vectors, and
the signals to be recovered may correspond to rows of the
channel transfer matrix [10]. In such settings, generalization of
low-complexity heuristic approaches based on sequential
gradient descent is challenging, as there might not be an
orthogonal transformation between the known signals (e.g.,
Fourier
transformation)
with
a
low-complexity
implementation (e.g., fast Fourier transform algorithm). The
concept of sequential algorithms may be generalized in the
form of alternating minimization, based on alternation
between the missing phase information and the candidate
solutions [20], but such methods do not reliably achieve good
performance [2,10].
Addressing phase retrieval problems in a general setting,

2
Candes et al. recently proposed an approach of lifting the
optimization problem to a high-dimensional space, where a
convex optimization problem can efficiently approximate the
original phase retrieval problem [21]. Convex optimizationbased phase retrieval has been demonstrated to achieve high
performance in imaging [21-23] and in communications [10].
Implementation complexity is critical in many applications,
yet the low-complexity implementation of convex
optimization-based phase retrieval has not yet been explored.
For example, in short-reach optical communications [10],
transceiver power consumption is an important concern, and
can determine the practical feasibility of any proposed
technology. Specialized low-complexity algorithms are
required for such applications. Methods employed in generalpurpose optimization software, such as interior-point methods,
are typically not suitable for such applications.
In this paper, we study three fundamental methods for lowcomplexity implementation of convex optimization-based
phase retrieval: the projected gradient method in Section III,
the Nesterov accelerated gradient method in Section IV, and
the alternating direction method of multipliers (ADMM) in
Section V. For each method, we describe adaptive algorithms
and present the associated computational complexities. In
comparing the performance of these methods, we focus on an
application area in optical communications: direct-detection
MDM, whose architecture is described in detail in Section VI.
We demonstrate that all these methods yield excellent
performance at low implementation complexity in either
transmitter- or receiver-based MIMO processing.

We focus on the problem of extracting unknown signal
vectors from the magnitudes of their inner products with
known signal vectors. We assume that D-dimensional known
complex vectors x (n ) are used for estimation of Ddimensional unknown complex vectors h l . With additive

(n ) , the available measurements are
h l H x (n )

2

+ n l (n ) = d l (n ),

N

min

rank(Cl ) + λ

∑ x( )

n H

Cl x (n ) − d l

(n )

2

,

n =1

(1)

where 1 ≤ n ≤ N and 1 ≤ l ≤ L . We might have L = D , e.g.,
in direct-detection MDM for a D × D MIMO channel, or we
might have L >> D , e.g., in diffractive imaging when
multiple illuminations are used.
In the absence of noise, the formulation (1) can be
expressed equivalently as L different rank-minimization
problems. Defining D × D complex Hermitian matrix
parameters C l , we have
min rank(Cl )
. (2)
(n )
s.t. Cl ≥ 0 and x(n )H Cl x(n ) = dl , for 1 ≤ n ≤ N
Optimal solutions to (2) are obtained for the matrix parameters

ˆ = h H h . To consider the effect of noise, the equality
C
l
l
l
constraints in (2) can be expressed as an additive penalty term
in a maximum-likelihood formulation of signal retrieval with
additive noise . In this paper, we focus on spatially white zero-

(3)

s.t. Cl ≥ 0
where λ is a tradeoff constant and ≥ denotes that the matrix
is positive semidefinite. Generalization of (3) to other noise
types is possible by considering the maximum-likelihood
formulation for the corresponding noise statistics, which might
or might not be convex in x (n ) .
As rank-minimization problems are non-deterministic
polynomial-time hard, it is desired to apply a convex
relaxation that can approximate the optimization problem (3)
efficiently, while yielding a global optimal solution. An
efficient convex relaxation is based on substituting the trace
function [21] in place of the rank function, which yields
N

min

g (Cl ) = tr (Cl ) + λ

∑ x( )

n H

Cl x ( n ) − d l

(n )

n=1

2

.

(4)

s.t. Cl ≥ 0

After computing the optimal matrix parameter Ĉ l , its first

ˆ = ∑D χ
ˆ ˆ ˆ H
ˆ 1,1 u
ˆ 1,1 , where C
χ
l
k =1 i ,k u i ,k u i ,k

principal component

, yields the estimate for h l . The overall dimensionality of the

problem is squared in exchange for the benefit of convex
approximation.
We will next focus on three techniques that aim to
approximately solve the increased-dimensionality problem
efficiently with low computational complexity. In estimating

() ()

II. CONVEX OPTIMIZATION-BASED PHASE RETRIEVAL

noise n l

mean real-valued Gaussian noise, for which the rankminimization problem can be expressed as

*

computational complexity, we assume ℜ , ℑ
and ( )
operations have negligible cost. We quantify the
computational complexity of other operations, which are
generally performed on complex numbers, in terms of the
number of floating point operations (flops) required. We
assume optimal ordering of algebraic operations to minimize
complexity. For example, multiplication of a matrix, vector
and scalar can be done with lower complexity by first
multiplying the matrix and the vector and then the scalar, or
by first multiplying the vector and the scalar and then the
matrix, as opposed to first multiplying the scalar and the
matrix and then the vector. We also assume that any computed
result can be reused at multiple steps without being
recomputed, by utilizing a working memory that can store
intermediate results. The structure of the objective function in
(4) forms the basis for the reduced-complexity algorithms.
III. PROJECTED GRADIENT METHOD
Projected gradient methods represent an efficient first-order
approach for constrained convex optimization problems. The
principle is based on iteratively searching for optimal values
of the variables along a direction determined by gradient
values, while imposing the constraints at each step.
Starting with an initial guess Cl 0 , the projected gradient
method applies the following update equation recursively at
iteration step k + 1 :
(5)
Cl k +1 ← P Cl k − t k ∇g Cl k .

[]

[

]

( [ ] [ ] ( [ ]))

3

( ) denotes the gradient operator corresponding to
objective function and P( ) denotes the projection

In (5), ∇g

the
operator that is required to satisfy the constraints.
We will next derive the various operators and analyze their
computational complexities. The step size can be updated by a
search algorithm. A common approach is a back-tracking line
search algorithm:
while ( g C l k < g Cl k − t k ∇g Cl k ):

( [ ])

( [ ] [ ] ( [ ]))

do t [k ] = t [k ] ⋅ τ
(6)
where 0 < τ < 1 . The computational complexity of backtracking line search is dominated by the complexity of the

(

)

objective function computation, which is O D 2 LN flops per
iteration step.
Computation of ∇g Cl , the gradient with respect to C l ,
requires computing the following derivatives for 1 ≤ i ≤ D and
i +1 ≤ k ≤ D :
2
∂g (C)
⎛N
⎞
= 1 + 2λ ⎜ ∑ x i(n ) ℜ x (n )H C l x (n ) ⎟
∂ℜ(C lii )
⎝ n =1
⎠
,
(7)
N
2
⎛
(n ) (n ) ⎞
− 2λ ⎜ ∑ ℜ d l
xi
⎟
⎝ n =1
⎠
N
∂g (C)
⎛
⎞
= 2λℜ⎜ ∑ x i(n ) x (kn )* x (n )H C l x (n ) ⎟
∂ℜ(C lik )
⎝ n =1
⎠
,
(8)
⎛ N (n ) (n ) (n )* ⎞
− 2λℜ⎜ ∑ d l x i x k ⎟
⎝ n =1
⎠
∂g (C)
⎛ N (n ) (n )* (n )H
⎞
= 2λℑ⎜ ∑ xi x k x Cl x(n ) ⎟
∂ℑ(Clik )
⎝ n =1
⎠
.
(9)
N
⎛
(n ) (n ) (n )* ⎞
− 2λℑ⎜ ∑ dl xi x k ⎟
⎝ n =1
⎠
2
(
n)
Computation of the x
and x (n )* x (n ) terms require

( )

(

)

(

)

( )

(

)

i

i

k

( )
2

O ND flops in total. For all expressions in (7)-(9), the
x (n )H Cl x (n ) terms are computed once with a complexity of
2

( )

O D2 L flops. Given

(n )* (n )
x i(n ) , x i x k , and x (n )H Cl x (n )

terms,

of

computation

N
(n ) (n )* (n )H C x (n )
∑ xi x k x
l

n =1
N

(n ) 2 ℜ x (n )H C x (n ) require
l

∑ xi

n =1

(

)

(

)

or

2 N − 1 flops for each i, k

(

)

and l, yielding a total complexity of O D2 LN flops. To
construct all the derivatives, performing all the additions and
subtractions in (7)-(9) requires O D2 L flops. The total

( )

(

)

complexity of computing all the derivatives is O D2 LN flops.
Imposing the constraint is based on the projection operator
onto the positive-semidefinite cone:
⎛D
⎞ D
P⎜⎜ ∑ λiuiui H ⎟⎟ = ∑ max(0, λi )uiui H .
(10)
⎝ i =1
⎠ i =1
Computation of the projection operator requires eigenvalue
decomposition. Hardware-efficient eigenvalue decomposition
techniques depend on the size of the matrix. For D = 3,
techniques based on the analytical expressions of the

eigenvalues and eigenvectors are fastest [1]. For D > 25,
divide-and-conquer is the fastest algorithm [2]. For 3 < D <
25, (the range of interest for practical applications, e.g., in
MDM systems), tridiagonal QR iteration is the fastest

( )

algorithm, with a complexity of O D3 flops [24]. After
replacing any negative eigenvalues by 0, the reconstruction of

( )

the projection requires O D3

flops. Overall, applying

( )

projection L times requires in total O D 3 L flops per
iteration,
For the projected gradient method, the overall
computational cost is dominated by the computation of
gradient
and
projection
operations,
and
is
flops per iteration. The total
max O D 2 LN , O D 3 L
computational complexity is also proportional to the number
of iterations. Hence, minimization of the total number of
iterations is also crucial.

((

) ( ))

IV. NESTEROV ACCELERATED GRADIENT METHOD
Optimal first-order methods are derived from classical
projected gradient methods, with the goal of reducing the
number of iterations required to achieve convergence.
One modification of the projected gradient method,
pioneered by Nesterov [25], is the accelerated gradient method
[26], which has update equations:

Cl [k + 1] ← P(K l [k ] − t[k + 1]∇g (K l [k ])),
γ[k + 1] ←

1
⎛1
⎜ + 1+ 1
⎜2
4 γ[k ]2
⎝

⎞
⎟
⎟
⎠

,

(11)
(12)

γ[k + 1](1 − γ[k ])
(C l [k + 1]− C l [k ]) . (13)
γ[k ]
For determining the step size in (11), as above, the backtracking line search algorithm (6) can be used. The
fundamental distinction of the Nesterov accelerated gradient
method is the adjustment of the step size using the γ k
K l [k + 1] ← C l [k + 1] +

[]

[]

parameter. The evolution of γ k as a function of the iteration
number is shown in Fig. 1.
The Nesterov accelerated gradient method has a slightly
higher complexity than the projected gradient method. The
computational cost for the update step (12) is O 1 flops and

( )

()

2

for (13) is O D L flops. The major computational cost is step

((

) ( ))

(11), which requires max O D 2 LN , O D 3 L flops per
iteration, similar to the projected gradient algorithm. Again,
the total computational complexity is proportional to the
number of iterations. The modification for acceleration
reduces the number of steps required for convergence, despite
a slight increase in the computational cost per iteration step.
V. ALTERNATING DIRECTION METHOD OF MULTIPLIERS
A general approach to address constrained convex
optimization problems with complicated objective functions is

4
1

N

1

( N ) = λ x (n ) x (n )* x (n ) x (n )* .
Tikab
∑ i k b a

(22)

0.8

0.6

0.6

0.4

0.4

0.2

0.2

γ[k]

0.8

0

0

5

10

15

20

(γ[k + 1](1 − γ[k]))/γ[k]

n =1

0

25

Iteration number
Fig. 1

γ[k ]

vs. iteration number for the Nesterov accelerated projected

The minimization problem (17) reduces to a set of D 2 linear
equations. Consider the matrix representation of the set of
linear equations (20)
⎡ C l11 [k + 1] ⎤
⎢ C [k + 1] ⎥
l12
⎢
⎥
~
⎥ = b l(N ) [k ],
(23)
X (N ) + ρI ⎢
⎢
⎥
⎢C lD(D −1) [k + 1]⎥
⎢ C lDD[k +1] ⎥
⎣
⎦

(

)

( )
where b l N [k ] is a vector of size D 2 × 1 such that
(N ) − (E [k ] / 2 + ρD [k ]) − δ / 2 ,
b l(n(()i −1)D + k ) [k ] = U lik
lik
lik
i ,k
N
T
~ (N )
2
2
X = λ ∑ X(n )X(n ) is a matrix of size D × D , and
n =1

gradient method.

to decompose them into smaller problems, each of which is
then easier to handle. The ADMM technique combines the
benefits of the dual decomposition and augmented Lagrangian
methods for constrained optimization,
We start by observing that the optimization problem can be
equivalently expressed in equality constraint form:

min

g (Cl ) + I ≥0 (Dl )

s.t. Cl = Dl

,

(14)

( )

2
)
(n )* (n )
X((in−1
)D + k = xi x k is a vector of size D × 1 . For the case of

(~

)

invertible X(N ) + ρI , the solution of (23) is given by

⎡ Cl11[k + 1] ⎤
⎢ C [k + 1] ⎥
⎢ l12
⎥ ~
−1
⎢
⎥ = X(N ) + ρI bl(N )[k ] = F (N )bl(N )[k ]. (24)
⎢
⎥
⎢ClD(D −1)[k + 1]⎥
⎢ ClDD[k +1] ⎥
⎣
⎦

(

)

where I ≥0 D l is the indicator function such that
⎧ 0, D l ≥ 0
.
(15)
I ≥0 (D l ) = ⎨
⎩∞, otherwise
The corresponding augmented Lagrangian is in the form
L(C l , D l , E l ) = g (C l ) + I ≥0 (D l )
.
(16)
2
+ El , Dl − Cl + ρ Dl − Cl F
ADMM is based on replacing minimization of (16) by
alternating minimization, where each step optimizes the
matrix parameters C l , D l and E l separately:

Inverting an ordinary matrix of size D 2 × D 2 requires O D 6
flops using naïve algorithms such as Gauss-Jordan

C l [k + 1] ← argmin L(C l , D l [k ], E l [k ]) ,

(17)

D l [k + 1] ← arg min L(C l [k + 1], D l , E l [k ]) ,

(18)

N −1
T
T
⎛
⎞
= ⎜ λX (N ) X (N ) + λ ∑ X (n ) X (n ) + ρI ⎟
n =1
⎝
⎠

E l [k + 1] ← E l [k ]+ ρ (Dl [k + 1]− Cl [k + 1]).

(19)

( )

(

~
F (N ) = X (N ) + ρI

(

Cl

Dl

We now analyze the steps (17)-(19).
For step (17), we find the analytical solution by setting the
first partial derivatives with respect to Clik to zero (since the
second derivatives are always positive), which yield
D

(

)

(N ) =
ρClik [k + 1] + ∑ Clab [k + 1]Tikab
a ,b =1

( N ) − ⎛ Elik [k ] + ρD [k ]⎞ − δ i , k
Ulik
⎜
⎟
lik
2
⎝ 2
⎠
where δ i ,k is a Kronecker delta function, and
N

( N ) = λ d (n ) x (n ) x (n )* ,
U lik
∑ l
i
k
n =1

,

(20)

)

elimination, and requires O D4.746 flops with optimized
Coppersmith-Winograd-like construction methods [27].
Iterative computation of the matrix inverse term F ( N ) is also
possible using the Sherman–Morrison–Woodbury lemma,
exploiting its structure:

= F (N −1) −

−1

)

−1

.

(25)

F (N −1) X (N ) X ( N )T F (N −1)
1 / λ + X (N )T F (N −1) X (N )

~

Note that we have the initial conditions X (0) = 0 and
F (0) = I / ρ . Computation of X (N ) requires O D 2 flops.
4
Given X (N ) and F (N −1) , the update of F ( N ) requires

( )

( )

OD

flops. As the update is repeated for all n , the total complexity
of computing F ( N ) becomes O D 4 N flops. Computation of

(

)

(

)

( N ) terms requires an additional
all U lik
O D2 NL flops. Given

( )

( N ), computation of all b (N ) [k ] vectors requires
U lik
O D2 L
l
flops.
(21)

Overall,

((

2

computation

) (

4

))

of

Cl [k + 1]

max O D LN , O D N flops at each iteration.

requires

5
For step (18), an analytical solution can be obtained by
equating the gradient to zero. We have the gradient
2 ⎞
⎛ 1
∇ D[k ]⎜
El [k ] + (Dl [k ] − Cl [k + 1]) ⎟ =
⎜ 2ρ
⎟
F⎠
⎝
.
(26)
⎛ 1
⎞
2⎜⎜
El [k ] + (Dl [k ] − Cl [k + 1])⎟⎟
2
ρ
⎝
⎠
1
If Cl [k + 1] −
El [k ] ≥ 0 , then the optimal solution is
2ρ
1
D l [k + 1] = C l [k + 1] −
E l [k ] . For the general case, the
2ρ

⎛
⎞
1
optimal D l k + 1 is the projection P⎜⎜ C l [k + 1] −
E l [k ]⎟⎟ .
2
ρ
⎝
⎠
The overall update equation becomes
⎛
⎞
1
(27)
D l [k + 1] ← P⎜⎜ C l [k + 1] −
E l [k ]⎟⎟ .
2ρ
⎝
⎠

[

modulus) of each modal field amplitude is available at the
receiver. This has important consequences for MIMO signal
(a)

(b)

]

Computation of the

( )

1
C l [k + 1] −
E l [k ] terms requires
2ρ

(c)

O D2 L flops in total. The projection operator is applied L
times and dominates the computational complexity, requiring
O D 3 L flops per iteration.

( )

( )
) ( ))

Lastly, step (19) requires O D2 L flops at each step. The
overall
computational
complexity
per
step
is
2
4
3
flops
per
iteration,
and
is
max O D LN , O D N , O D L
again proportional to the total number of iterations.
The computational complexities of the three methods are
summarized in Table I.

((

) (

TABLE I
COMPUTATIONAL COMPLEXITIES OF PHASE RETRIEVAL METHODS
Computational complexity per iteration

( (
( (

) (
) (

))
))

Projected gradient

max O D 2 LN , O D 3 L

Nesterov

max O D 2 LN , O D 3 L

ADMM

max O D 2 LN , O D 4 N , O D 3 L

((

) (

) (

))

VI. APPLICATION EXAMPLE: DIRECT-DETECTION MODEDIVISION MULTIPLEXING

In this section, we study application of the phase retrieval
techniques described above in direct-detection MDM [10].
MDM aims to increase the data throughput per fiber by
employing a fiber supporting D propagating modes. At the
transmitter, D independent data signals are multiplexed onto
these modes, ideally providing a D-fold increase in
throughput. During transmission, the data signals are
intermixed by random coupling between modes, which
corresponds to a D × D channel transfer matrix that has offdiagonal entries. We assume that group velocity dispersion
and modal dispersion are negligible, such that memoryless
signal processing and channel estimation are sufficient.
Using direct detection simplifies receiver implementation,
making the technique potentially suitable for short-reach
communications. However, only the intensity (squared

Fig. 2 Direct detection MDM systems using: (a) digital MIMO precoding at
the transmitter, (b) optical MIMO precoding at the transmitter, and (c)
optical MIMO equalization at the receiver. Mux: Multiplexer, mod.:
modulator, det.: detector, I: in-phase, Q: quadrature.

processing, channel transfer matrix estimation and information
decoding. First, compensating for modal crosstalk requires
complex-valued D × D MIMO signal processing. As shown
in Fig. 2, direct-detection MDM can be implemented using a
digital MIMO precoder or optical MIMO precoder at the
transmitter, or using an optical MIMO equalizer at the receiver
(additional details are given in [10]). Second, adjustment of
the MIMO precoder or equalizer requires estimation of the
D × D channel transfer matrix using phase retrieval. Third,
information must be demodulated from the optical carrier
using noncoherent or differentially coherent detection.
Here, we are concerned with the second consequence,
namely, how to perform channel estimation using phase
retrieval. Estimation is facilitated using a known training
sequence, which is composed of N D-dimensional signals. For
a given D , the required size of the training signals can be
found empirically. If the number of training vectors N is not
large enough, high estimation error is obtained due to
insufficient measurement diversity. On the other hand,
increasing N beyond a certain threshold does not improve
estimation accuracy significantly, and it is desired to set the
value of N around this threshold. Note that for optical
precoding or equalization methods, the training vectors can be
transmitted at a low symbol rate, and high-speed digital signal
processing can be avoided. For all three methods (digital or
optical MIMO precoding or optical equalization), channel
estimation corresponds to the phase retrieval problems
described in Section 2.

6
For a MIMO precoder, the channel transfer matrix needs to
be estimated up to a constant diagonal unitary matrix [10]. We
have L = D , since for each transmitted sequence, D different

the receiver is thermal noise-limited, so the received modal
intensities are corrupted by additive, spatially white real
Gaussian noise. Realizations of the channel transfer matrix are

Group 1

(a)

2×2
MIMO

MMF
Spatial
Demux

Group 2

SNR = 12 dB

-1

a2

10

SNR = 14 dB

b1
4×4
MIMO

Crosstalk-free
Interior point-based
Projected gradient
Nesterov accelerated gradient
ADMM

SNR = 8 dB

a1

b2

-2

10

b3

(b)

Polarization
Beam
Spatial Splitter
Mode 1

4 × 4 MIMO Processing using Mach-Zehnder Array

Average BER

b4
-3

10

-4

10

-5

10

Spatial
Mode 2

-6

10

-7

10
Polarization
Converter

Fig. 3 (a) Demultiplexing using receiver MIMO processing in mode group
subspaces, (b) implementation of optical MIMO processing using MachZehnder array.

inner products are generated by the D rows of the channel
transfer matrix, i.e., h l in (1) corresponds to the l th row of the
channel transfer matrix.
For a MIMO equalizer, the channel transfer matrix needs
to be estimated up to a complex exponential, hence an
additional optimization procedure is used to estimate the
unknown phases of each row of the channel transfer matrix,
yielding L = D + 1 [10].
Considering L = D or L = D + 1 with Table I, for a MIMO
precoder or a MIMO equalizer, we estimate the computational
complexity per iteration to be max O D 3 N , O D 4 for the
projected gradient and Nesterov accelerated gradient methods,
and O D 4 N for the ADMM. In other words, they are of the
order of the dimensionality of the high-dimensional space,
which is the square of the number of unknowns.
In short-reach systems, the end-to-end link may be
optimized to have negligible modal crosstalk between mode
groups [28]. In this case, MIMO processing can be
implemented within each mode group separately. Fig. 3
illustrates processing in mode group subspaces in a system
using receiver MIMO equalization implemented by an array of
Mach-Zehnder interferometers [29]. Processing in mode group
subspaces can reduce the total hardware complexity of MIMO
processing, as the required number of Mach-Zehnder
interferometers scales with the square of the number of modes
processed. In addition, mode group subspace processing can
significantly reduce the computational complexity of channel
estimation, as the estimation algorithms can be implemented
separately for each mode group (e.g., of sizes D = 2 and
D = 4 for the system shown in Fig. 3).

((

(

) ( ))

)

VII. NUMERICAL RESULTS
In the simulations presented here, we consider a directdetection MDM system employing on-off keying. We assume

0

10

20

30
40
Number of iterations

50

60

70

Fig. 4 Average BER vs. number of iterations for a direct-detection MDM
system using on-off keying, transmitter MIMO precoding, D = 6 modes and
N = 300 training signals. Horizontal lines show the performance achieved in
the two reference cases: for the ideal crosstalk-free channel and when an
interior-point technique [28] is used to implement convex estimation-based
phase retrieval. ADMM: alternating direction method of multipliers.

chosen from an ensemble of fully random complex-valued
D × D unitary matrices to model worst-case modal crosstalk.
We assume D = 6 modes and, accordingly, we choose
N = 300 training signals to obtain low estimation error. We
assume that these N training signals are generated randomly,
with each vector element an independent, identically
distributed circularly symmetric Gaussian random variable.
The convex optimization formulation tradeoff constant is
chosen as λ = 10, as in [10], as this choice yielded the lowest
estimation error in the original convex optimization problem
formulation (4). For projected gradient and Nesterov
accelerated gradient methods, the backtracking line search
parameter is chosen as τ = 0.3, which is at an optimal point
when the tradeoff between the number of backtracking line
search iterations and the direction of search accuracy is
considered (when the value of τ is close to 1, convergence of
the backtracking line search would be slow and when the
value of τ value is too small, the estimation error would be
high when moving in the direction of the step size).
The system electrical signal-to-noise ratio (SNR) is defined
as

⎡
SNR = ⎢
⎣
where

D

∑i =1 y i

(n )

2

(D ⋅ σ )⎤⎥

2

,

(28)

⎦

denotes an average over training signals. Since (28) is

proportional to the square of the optical power per mode, a 1-dB
change of the optical power per mode corresponds to a 2-dB
change of 10 log10 (SNR ). A range of SNR values is

considered that corresponds to practical short-reach optical
links. Link performance is quantified by the average bit-error
ratio (BER). Note that short-reach communication systems
typically use forward error-correction codes with low
complexity, necessitating BER targets around 10-3 to 10-4.

7
We compare the performance of the three low-complexity
convex optimization-based phase retrieval techniques to the
performance achieved in two reference cases. The first

Average BER

10-2

Crosstalk-free
Interior point-based
Projected gradient
Nesterov
ADMM

10-3

10-4

10-5

7

8

9

10

11

12

13

14

SNR (dB)

Fig. 5 Average BER vs. SNR for a direct-detection MDM system using onoff keying, transmitter MIMO precoding, D = 6 modes and N = 300 training
signals. ADMM: alternating direction method of multipliers.

Average BER

10-2

Crosstalk-free
Interior point-based
Projected gradient
Nesterov
ADMM

10-3

10-4

10-5
7

8

9

10
11
SNR (dB)

12

13

14

Fig. 6 Average BER vs. SNR for a direct-detection MDM system using onoff keying, receiver MIMO equalization, D = 6 modes and N = 300 training
signals. ADMM: alternating direction method of multipliers.

reference case is the ideal crosstalk-free channel, where
additive noise is the only impairment. The second reference
case includes modal crosstalk and uses CVX [30], a wellknown interior-point-based package, for convex estimationbased phase retrieval to determine the MIMO precoder or
equalizer.
Fig. 4 shows the average BER vs. number of iterations
using transmitter MIMO precoding. The convergence speed is
rather insensitive to the SNR value. The projected gradient
method converges after ~50 iterations, the Nesterov
accelerated gradient method converges after ~10 iterations,
and the ADMM converges after just a few iterations. In the
remainder of the simulations, we fix the total number of
iterations to 70 for the projected gradient method, 30 for the
Nesterov accelerated gradient method and 5 for the ADMM,
to ensure convergence beyond the respective knees of the
curves shown in Fig. 4.
Table II provides estimates of the number of operations per
second required for channel estimation based on the
assumptions described in Sections III-V. We note that the
mapping of an algorithm to hardware can be optimized given
tradeoffs between the costs of computation, memory, and
communication on a chip. Hence, the numerical values in
Table II should be considered as rough order-of-magnitude

estimates for straightforward hardware mapping on a singlecore processor without parallelization of the computation.
The operation counts in Table II are computed assuming
channel estimation at intervals of 1, 10 and 100 ms. Based on
studies of channel dynamics [31], we estimate that short-reach
multimode links with D~6 modes will change on time scales
of ms to tens of ms [31]. Even if a direct-detection MDM
system must update channel estimates at 1 ms intervals, the
operation count for ADMM is only ~5 Gflop/s. This is about
two orders of magnitude smaller than the operation count for a
long-haul coherent MDM receiver [28].
TABLE II
OPERATIONS PER SECOND FOR CONVEX OPTIMIZATION-BASED PHASE
RETRIEVAL IN A DIRECT-DETECTION MDM SYSTEM USING TRANSMITTER
MIMO PRECODING, ASSUMING L = D = 6, N = 300.

Update interval
Projected gradient
Nesterov
ADMM

Operations per second required for channel
estimation for update intervals indicated
1 ms
10 ms
100 ms
~50 Glops/s
~5 Gflops/s
~500 Mflops/s
~20 Gflops/s
~2 Gflops/s
~200 Mflops/s
~5 Gflops/s
~0.5 Gflops/s
~50 Mflops/s

Fig. 5 shows the average BER vs. SNR using transmitter
MIMO precoding. All three techniques exhibit performance
very close to CVX, incurring SNR penalties (with respect to
crosstalk-free case) less than 0.3 dB for target average BERs
of 10-3 to 10-4.
Fig. 6 shows the average BER vs. SNR using receiver
MIMO
equalization.
Receiver
MIMO
equalization
necessitates estimation procedures to determine the transmitter
MIMO precoder, as well as additional estimation to
disambiguate common phase factors [10]. Consequently, there
are more opportunities for estimation errors than in transmitter
precoding. The projected gradient and Nesterov accelerated
gradient methods yield similar performance, incurring SNR
penalties (with respect to crosstalk-free case) less than 0.6 dB
for target average BERs of 10-3 to 10-4. For ADMM, the SNR
penalty is slightly higher, but still less than 0.8 dB for the
same target average BER range.
VIII. CONCLUSIONS
We have studied three low-complexity approaches for
implementing convex optimization-based phase retrieval: the
projected gradient method, the Nesterov accelerated gradient
method and ADMM. The projected gradient method is a firstorder technique, which aims to impose the constraint at each
iteration step. The Nesterov accelerated gradient method
achieves an acceleration of the projected gradient method by
using an adaptive step size parameter. The ADMM is based on
redefining the convex optimization objective for uncoupled
solutions of the separable objective functions using alternating
minimization. We have demonstrated application of the three
techniques to direct-detection MDM. All three methods
accurately approximate the solution to the original phase
retrieval problem, incurring small estimation penalties. The
computational complexities of all three methods scale with the
square of the number of unknowns per iteration. Among the

8
three, ADMM yields convergence after the smallest number of
iterations.
ACKNOWLEDGEMENTS
Discussions with Prof. Emmanuel J. Candes are gratefully
acknowledged.
REFERENCES
[1] A. Walther, “The question of phase retrieval in optics,” J. Mod. Opt. vol.
10, no. 1, pp. 41-49, 1963.
[2] E. J. Candes, X. Li, and M. Soltanolkotabi, “Phase retrieval via Wirtinger
flow: Theory and algorithms,” IEEE Trans. Inf. Theory, vol. 61, no. 4, pp.
1985-2007, Feb. 2015.
[3] F. Fogel, I. Waldspurger, A. d’Aspremont, “Phase retrieval for imaging
problems,” Math. Program. Comp. vol. 8, no. 3, pp. 311-355, Mar. 2016.
[4] R. P. Millane, “Phase retrieval in crystallography and optics,” J. Opt. Soc.
Am. A vol. 7, no. 3, pp. 394-411, 1990.
[5] R. Horisaki, R. Egami, J. Tanida, “Single-shot phase imaging with
randomized light (SPIRaL),” Opt. Expr. vol. 24, no. 4, pp. 3765-3773, Feb.
2016.
[6] M. H. Seaberg, A. d'Aspremont, and J. J. Turner, “Coherent diffractive
imaging using randomly coded masks,” App. Phys. Lett. vol. 107, no. 23,
231103, Nov. 2015.
[7] J. A. Rodriguez, R. Xu, C. C. Chen, Y. Zou, and J. Miao, “Oversampling
smoothness: an effective algorithm for phase retrieval of noisy diffraction
intensities,” J. App. Crystallogr. vol. 46, pp. 312-318, 2013.
[8] J. C. Dainty, and J. R. Fienup, “Phase retrieval and image reconstruction
for astronomy,” in Image Recovery: Theory and Applications, H. Stark,
(Academic Press Inc., 1987), pp. 231-275.
[9] F. Roddier, Adaptive Optics in Astronomy (Cambridge University Press,
1999).
[10] S. O. Arik and J. M. Kahn, “Direct-detection mode-division multiplexing
in modal basis using phase retrieval,” Opt. Lett. vol. 41, no. 18, pp. 42654268, Sep. 2016.
[11] R. A. Gonsalves and R. Chidlaw, “Wavefront sensing by phase retrieval,”
Proc. SPIE 0207, Applications of Digital Image Processing III no. 32, Dec.
1979.
[12] G. Wunder, H. Boche, T. Strohmer, and P. Jung, “Sparse signal
processing concepts for efficient 5G system design,” IEEE Access vol. 3, pp.
195-208 Feb. 2015.
[13] R. Balan, P. Casazza, and D. Edidin, “On signal reconstruction without
phase,” App. and Comp. Harmonic Anal. vol. 20, no. 3, pp. 345-356, May
2006.
[14] M. B. McCoy and J. A. Tropp, “The achievable performance of convex
demixing,” [Online]. Available: https://arxiv.org/abs/1309.7478.
[15] C. Hegde and R. B. Baraniuk, “Signal recovery on incoherent
manifolds,” IEEE Trans. on Info. Theory vol. 58, no. 12, pp/ 7204-7214, Jul.
2012.
[16] P. Schniter and S. Rangan, “Compressive phase retrieval via generalized
approximate message passing,” IEEE Trans. Signal Process. vol. 63, no. 4,
pp. 1043-1055, Feb. 2015.
[17] R. W. Gerchberg and W. O. Saxton, “A practical algorithm for the
determination of the phase from image and diffraction plane pictures,” Optik
vol. 35, pp. 237-246, 1972.
[18] J. R. Fienup, “Phase retrieval algorithms: a comparison,” App. Opt. vol.
21, no. 15, pp. 2758-2769, Aug. 1982.
[19] V. Elser, “Phase retrieval by iterated projections,” J. Opt. Soc. Am. A vol.
20, no. 1, pp. 40-55, 2003.
[20] P. Netrapalli, P. Jain, and S. Sanghavi, “Phase retrieval using alternating
minimization,” IEEE Trans. Signal Process. vol. 63, no. 18, pp. 4814-4826,
Jun. 2015.
[21] E. J. Candes, T. Strohmer, and V. Voroninski, “Phaselift: Exact and
stable signal recovery from magnitude measurements via convex

programming,” Commun. Pure and App. Math. vol. 66, no. 8, pp. 1241-1274,
Aug. 2013.
[22] E. J. Candes, Y. C. Eldar, T. Strohmer, and V. Voroninski, “Phase
retrieval via matrix completion,” SIAM J. Imaging Sci. vol. 57, no. 2, pp. 225251, 2015.
[23] L. H. Yeh, J. Dong, J. Zhong, L. Tian, M. Chen, G. Tang, M.
Soltanolkotabi, and L. Waller, “Experimental robustness of Fourier
ptychography phase retrieval algorithms,” Opt. Expr. vol. 23, no. 26, pp.
33214-33240, Dec. 2015.
[24] J. W. Demmel, Applied Numerical Linear Algebra (Siam, 1997).
[25] Y. Nesterov, Introductory Lectures on Convex Optimization, (Kluwer
Academic Publishers, 2004).
[26] A. Beck and M. Teboulle, “A fast iterative shrinkage-thresholding
algorithm for linear inverse problems,” SIAM J. Imaging Sci. vol. 2, no. 1, pp.
183–202, 2009.
[27] W. Williams, “Multiplying matrices faster than Coppersmith-Winograd,”
in Proceedings of ACM Symposium on Theory of Computing 2012, pp. 887898.
[28] S. O. Arık, K. P. Ho, J. M. Kahn, “Group delay management and multiinput multi-output signal processing in mode-division multiplexing systems,"
J. Lightw. Techn. vol. 34, no. 11, pp. 2867-2880, Jun. 2016.
[29] D. A. B. Miller, “Self-aligning universal beam coupler,” Opt. Expr. vol
21, no. 5, pp. 6360-6370, Mar. 2013.
[30] M. Grant and S. Boyd, "CVX: Matlab Software for Disciplined Convex
Programming," [Online]. Available: http://cvxr.com/cvx.
[31] K. Choutagunta and J. M. Kahn, “Dynamic Channel Modeling for ModeDivision Multiplexing", J. Lightw. Techn. vol. 35, no. 12, pp. 2451-2463, Jun.
2017.

