# Quantization-Aware Phase Retrieval

We address the problem of phase retrieval (PR) from quantized measurements. The goal is to reconstruct a signal from quadratic measurements encoded with a finite precision, which is indeed the case in many practical applications. We develop a rank-1 projection algorithm that recovers the signal subject to ensuring consistency with the measurement, that is, the recovered signal when encoded must yield the same set of measurements that one started with. The rank-1 projection stems from the idea of lifting, originally proposed in the context of PhaseLift. The consistency criterion is enforced using a one-sided quadratic cost. We also determine the probability with which different vectors lead to the same set of quantized measurements, which makes it impossible to resolve them. Naturally, this probability depends on how correlated such vectors are, and how coarsely/finely the measurements get quantized. The proposed algorithm is also capable of incorporating a sparsity constraint on the signal. An analysis of the cost function reveals that it is bounded, both above and below, by functions that are dependent on how well correlated the estimate is with the ground truth. We also derive the Cramér-Rao lower bound (CRB) on the achievable reconstruction accuracy. A comparison with the state-of-the- art algorithms shows that the proposed algorithm has a higher reconstruction accuracy and is about 2 to 3 dB away from the CRB. The edge, in terms of the reconstruction signal-to-noise ratio, over the competing algorithms is higher (about 5 to 6 dB) when the quantization is coarse.

There are no comments yet.

## Authors

• 6 publications
• 15 publications
• ### Taking the edge off quantization: projected back projection in dithered compressive sensing

Quantized compressive sensing (QCS) deals with the problem of representi...
05/11/2018 ∙ by Chunlei Xu, et al. ∙ 0

• ### Phase retrieval with Bregman divergences: Application to audio signal recovery

Phase retrieval aims to recover a signal from magnitude or power spectra...
11/25/2020 ∙ by Pierre-Hugo Vial, et al. ∙ 0

• ### Quantized Compressive Sensing with RIP Matrices: The Benefit of Dithering

In Compressive Sensing theory and its applications, quantization of sign...
01/17/2018 ∙ by Chunlei Xu, et al. ∙ 0

• ### Robust Wirtinger Flow for Phase Retrieval with Arbitrary Corruption

We consider the phase retrieval problem of recovering the unknown signal...
04/20/2017 ∙ by Jinghui Chen, et al. ∙ 0

• ### DeepFPC: Deep Unfolding of a Fixed-Point Continuation Algorithm for Sparse Signal Recovery from Quantized Measurements

We present DeepFPC, a novel deep neural network designed by unfolding th...
12/02/2019 ∙ by Peng Xiao, et al. ∙ 0

• ### Total least squares phase retrieval

We address the phase retrieval problem with errors in the sensing vector...
02/01/2021 ∙ by Sidharth Gupta, et al. ∙ 0

• ### Estimation from Quantized Gaussian Measurements: When and How to Use Dither

Subtractive dither is a powerful method for removing the signal dependen...
11/16/2018 ∙ by Joshua Rapp, et al. ∙ 0

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## I Introduction

The objective of phase retrieval (PR) is to estimate a signal from intensity measurements given as

 bi=∣∣a⊤ix∗∣∣2,i=1:m, (1)

where are known sampling vectors. For a complex vector , the notation denotes its Hermitian transpose. In addition, the measurements may get corrupted by noise. The measurement model considered in (1) arises in a number of imaging applications, such as X-ray crystallography [1], holography [2], electron microscopy [3]

, etc. For example, the diffraction patterns of objects to be imaged using X-ray crystallography closely approximate their Fourier transforms. The sensors can only record the intensities of the complex wave-field; and the phase, which contains critical structural information about the object, is not directly measured. Thus, it becomes imperative to recover the phase, starting from the Fourier magnitude/intensity measurements, in order to reconstruct the object accurately. The fundamental objective of PR is to solve this otherwise ill-posed inverse problem by taking into account prior information about the underlying signal, such as non-negativity, compact support, sparsity, etc. One could also resolve the phase ambiguity by considering oversampled measurements exceeding the signal dimension (

). In the special case where the vectors correspond to the discrete Fourier transform (DFT) basis vectors, the PR problem reduces to classical Fourier PR, wherein one seeks to reconstruct a signal starting from its Fourier magnitude or intensity. The generalized setting, which involves projections onto random sampling vectors is the one that we shall consider. Before giving a formal statement of the problem considered in this paper, we provide a concise review of the existing PR literature to put our contributions into perspective.

### I-a A Survey of Phase Retrieval Literature

The PR problem has its origin in optics and astronomy. The initial contributions were due to Fienup [4, 5], and Gerchberg and Saxton [6], who proposed iterative error-reduction algorithms that bounce estimates back and forth between the object and the measurement domains, and incorporate respective priors. The most widely used priors in the signal domain are causality, non-negativity, compact support, sparsity, etc. Fienup’s algorithm has been the most popular technique for PR in the optics community, and works reasonably well for a wide class of imaging problems. A comprehensive overview of the Fienup algorithm and several of its variants can be found in [7] and the references therein. A notable variant of the Fienup algorithm was developed by Quatieri et al. [8] for reconstructing minimum-phase signals from their DFT magnitude measurements, wherein one iteratively imposes the causality constraint in the signal domain, and combines the measured magnitude spectrum with the current estimate of phase in the frequency domain. Apart from the iterative algorithms, there exist non-iterative techniques [9], which rely on the Hilbert transform relationship between the log-magnitude and the phase of the Fourier transform of minimum-phase signals, in order to reconstruct them from magnitude-/phase-only measurements. The two-dimensional (2-D) counterpart of such results and exact reconstruction guarantees were proposed in [10] in the context of digital holography. We recently developed a non-iterative algorithm [11, 12]

to solve the PR problem for a class of 2-D parametric models, by extending the concept of minimum-phase signals in 1-D. An exact PR methodology for signals lying in shift-invariant spaces was developed in

[13]. We also recently constructed generalized minimum-phase signals and developed corresponding 2-D Hilbert integral equations [14].
Moravec et al. addressed the problem of PR within the realm of sparsity and magnitude-only compressive measurements [15]. The compressive PR problem received considerable attention because of its wide applicability. Signals encountered in a number of applications indeed admit a sparse representation in an appropriately chosen basis. Yu and Vetterli proposed a sparse spectral factorization technique [16], and established uniqueness guarantees, where the objective was to recover a sparse signal from its autocorrelation sequence. A greedy local search-based algorithm for sparse PR, referred to as GESPAR, was proposed by Schechtman et al. [17]. The scalability and accuracy of GESPAR has been established in [17] by extensive simulations in a variety of practical settings. Netrapalli et al. [18] developed an analytical convergence guarantee for the well-known alternating minimization (Alt. Min.) framework for PR, with and without the constraint of sparsity, referred to as AltMinPhase and SparseAltMinPhase, respectively. Their work is the first one in the literature to establish the correctness of Alt. Min. for PR, subject to the so-called spectral initialization [18]. Vaswani et al. [19] recently proposed an Alt. Min. technique for low-rank PR, where the task is to recover a low-rank matrix from the quadratic measurements corresponding to projections with each of its columns. Other notable contributions for sparsity regularized PR include dictionary learning for PR (DOLPHIn) [20], cone programming [21], compressive PR via generalized message passing [22], simulated annealing for sparse Boolean signals [23], majorization-minimization for recovery from under-sampled measurements [24], etc. Fogel et al. [25] have recently shown that incorporating signal priors, such as sparsity and positivity lead to a significant speed-up of iterative reconstruction techniques.
A distinctive contribution in PR is the PhaseLift framework pioneered by Candès et al. [26, 27]. PhaseLift relies on the idea of lifting a vector to a matrix such that the quadratic measurements of get converted to an equivalent set of linear measurements of . Reconstruction is achieved by solving a tractable semi-definite program (SDP). One encounters two scenarios within this framework: (i) the absence of a signal prior, together with a set of oversampled measurements; or (ii) reconstruction subject to a signal prior such as sparsity. Ohlsson et al.’s compressive PR with lifting (CPRL) [28] technique falls in the second category, where sparsity is enforced via an penalty. Schechtman et al. [29] also developed a sparse PR technique in the context of sub-wavelength imaging with partially incoherent light by employing the idea of PhaseLift, wherein the sparsity constraint on the underlying image is imposed via a log-det penalty. The PhaseLift technique requires spectral decomposition of an matrix within every iteration, where is the dimension of the signal to be recovered, and, therefore, its complexity per iteration scales as

. However, since one needs to compute only the top few eigenvectors of a symmetric matrix, one could employ

power iterations [30, Chapter 7] to significantly speed-up the computation. Gradient-descent approaches for PR without lifting include the Wirtinger flow (WF) method [31] and its truncated version (TWF) [32]. These algorithms lead to stable reconstruction, and have convergence guarantees provided that the starting point is accurate, which is typically achieved by using spectral initialization. The WF and TWF algorithms are also scalable with respect to the signal dimension. Waldspurger et al. developed the PhaseCut technique [33], where the PR problem is formulated as a non-convex quadratic program and solved using a provable block-coordinate-descent approach. Although the number of variables in the resulting SDP in PhaseCut is larger than that in PhaseLift, it has been shown in [33] that the proposed algorithm for PhaseCut has a per-iteration complexity comparable with that of the iterative error-reduction type algorithms.
The issue of noise robustness of several PR algorithms has been considered, but the effect of quantization, which is ubiquitous in practical acquisition systems, has not been addressed in the PR literature. In contrast, a significant amount of research has gone into developing algorithms for compressive sensing (CS) reconstruction with quantized measurements starting from the work of Zymnis et al. [34], who developed two reconstruction algorithms for quantized CS based on -regularized maximum likelihood (ML) and least-squares estimation. Laska et al. [35] studied the effect of quantization and saturation on linear compressive measurements and a bit-precision analysis for CS was carried out by Ardestanizadeh et al. [36]. The extreme case of binary/one-bit quantization, wherein one retains only the sign of the measurements, has been investigated in detail in the context of CS [37, 38, 39, 40, 41]. However, the issue of measurement quantization and its ramifications on the reconstruction algorithms have not been considered in the PR literature. Our attempt in this paper is to fill this void.

### I-B Our Contributions

We introduce the problem of QPR, that is phase retrieval from quantized measurements (Section II) and address issues related to distinguishability (Section III) of two distinct signals in terms of their quantized quadratic measurements. Issues related to noise robustness and quantizer design are also part of this development. Subsequently, the optimization framework for QPR is developed by combining the principles of consistent recovery and lifting, and we propose two gradient-based iterative projection algorithms for signal reconstruction (Section IV). Our algorithms are amenable to incorporating the sparsity prior as well, which boosts the reconstruction signal-to-noise ratio (SNR). The descent property of the projected gradient approach for QPR is established in Section IX-D of the supplementary material. An analysis of the cost function shows that it can be bounded both above and below, in a probabilistic setting, as a function of the distance between the ground-truth and the estimate, and parameters that measure the precision of the quantizer (Section IX-B, supplementary material). Since no explicit quantization-aware PR algorithms exist in the literature, we consider the quantization noise to be additive as far as the implementation of the state-of-the-art PR algorithms (such as PhaseLift, truncated Wirtinger flow (TWF), AltMinPhase, GESPAR, compressive PR with lifting (CPRL)) are considered, and make performance comparisons for the QPR problem first without a sparsity prior (Section V) and subsequently with the sparsity prior incorporated (Section VII). For benchmarking the performance, we derive the Cramér-Rao bounds (CRB) assuming white Gaussian noise contamination prior to quantization (Section IX-C, supplementary material). The comparisons show that the proposed algorithm achieves a reconstruction mean-squared error (MSE) that is within dB of the CRB and about dB better than the best performing technique from the state-of-the-art (Section VI). Our recent work on PR from binary measurements [44] follows as a special case of the formalism developed in this paper.

## Ii Measurement Model for QPR

Recall that the objective of PR is to reconstruct a signal from intensity measurements of the form (1), which encompasses the special case of signal reconstruction from DFT intensity measurements. In this paper, we consider the scenario where the intensity measurements are acquired with a finite precision, that is, they are quantized using a finite codebook , containing distinct symbols. The acquired measurements take the form

 (2)

where the encoding map is determined using a set of thresholds , such that , whenever , for . A codebook of size requires bits for quantization. The measurement model (2) is more practical than (1) as any real-world measurement device would have a finite bit precision. In principle, for the consistent recovery framework pursued in this work, the encoding symbols could be arbitrary as long as the association between and is known. However, for the sake of simplicity, and for making comparisons with the state-of-the-art techniques that use the encoded measurement values, we assume that

 τj−1

meaning that the encoding symbol for an interval is chosen to be a point falling inside that interval. Since the state-of-the-art PR algorithms aim to minimize an appropriate loss function in the measurement domain, selecting the encoding symbols

from the corresponding intervals is a reasonable choice. As we shall see, the encoding symbols have no bearing on the reconstruction performance of the proposed approach. However, it would affect the performance of the state-of-the-art techniques.

### Ii-a The Principle of Consistent Recovery

Since quantization is a non-invertible mapping, it would not be possible to determine exactly. In other words, there could be two candidate signals and such that their quantized intensity measurements match exactly and one needs to design an optimization objective that does not distinguish between them. If one uses the conventional squared-error loss , wherein the intensity measurements of and are compared with the acquired measurement , one might get different values of the error corresponding to and , although the measurement process does not distinguish between them. Therefore, instead of minimizing the traditional squared-error loss, we seek a solution that is consistent with the measurements. The idea of consistency has earlier been used in the context of reconstruction from quantized CS samples [39]. To elucidate the idea of consistent recovery in the PR context, let

 Hj={ai,1≤i≤m\,|\,% yi=sj},j=1:k, (4)

denote the collection of sampling signals such that the corresponding measurements are encoded as . A solution is said to be consistent with the measurements if

 τj−1≤∣∣a⊤i^x∣∣2<τj,\,\,whenever\,\,ai∈Hj. (5)

Essentially, the consistency condition in (5) ensures that the reconstructed vector, when passed through the same acquisition process, matches the measurements that one started with. One can impose the constraint of sparsity depending on the application. Accordingly, we have two problem statements as specified below.

1. Quantized PR (QPR): Find such that

2. Sparse QPR (SQPR): Find satisfying such that

### Ii-B Performance Metrics

If is a consistent solution to any of the PR problems posed above, so is . In order to factor out the effect of the global sign, an appropriate metric to quantify the accuracy of reconstruction vis-à-vis the ground truth would be the global-sign-invariant reconstruction SNR [27] defined as111For a complex ground-truth signal , the global-phase-invariant reconstruction SNR can be defined as .:

 Reconstruction SNR=maxα∈{−1,+1}∥x∗∥22∥α^x−x∗∥22. (6)

The MSE of reconstruction is defined as the reciprocal of the SNR metric. The second metric that would be relevant in the context of QPR is the consistency (denoted as ) of the reconstruction with the measurements, defined as

 (7)

where is the indicator of the event . Essentially, the consistency metric quantifies the fraction of measurements correctly explained by . Naturally, , and is the best that one could hope to achieve.

## Iii Effect of Quantization on the Measurements

We next address the issues of distinguishability, noise-robustness of the quantized measurements, and the design of the quantizer. Akin to [32], one could consider two sampling models: (i) the real model, where ; and (ii) the complex model, where the test vectors can be expressed as , , with and drawn independently drawn from . The subscripts re and im denote the real and imaginary parts, respectively, whereas denotes the identity matrix. The signal is assumed to have unit norm in both models, without loss of generality. We consider the real signal model throughout our development and the analysis carried out in this section is only applicable to the real model. The notation denotes the unit sphere in . In the following proposition, we derive the distribution of the full-precision quadratic measurement.

###### Proposition 1

The quadratic measurement is

distributed with one degree of freedom for any

and , that is .

Proof: Observe that

 ∣∣a⊤x∗∣∣2=a⊤x∗x∗⊤a=a⊤UΛU⊤a=~a21,

where , and follows the same distribution as , since is orthonormal. The proposition is now a direct consequence of the fact that is the square of a random variable.

The cumulative distribution function (c.d.f.) of

is given by

 Fb(b)=γ(b2,12)√π, (8)

where denotes the lower incomplete gamma function.

### Iii-a Distinguishability of the Quantized Measurements

Let and be two linearly independent signals in , meaning that . We analyze the probability of and being mapped to the same set of quantized measurements, which renders them indiscernible. For any reconstruction algorithm to succeed, it is necessary to collect a sufficient number of measurements to keep low. To begin with, we place an upper-bound on the error probability in the case where only one measurement is acquired. Let and . It would be impossible for any algorithm to differentiate from from their quantized measurements if

 Q(b1)=Q(b2). (9)

Assuming noise-free measurements, there are two events that could possibly lead to (9): (i) where both and fall within the same quantization bin; and (ii) where both and are saturated, that is, . The probability of the first event is upper-bounded by that of the event , where

 δ=max1≤j≤k−1\,\,(τj−τj−1), (10)

is the the precision of the quantizer, indicating that and are closer apart than the precision of the quantizer. The probability of the second event is denoted as , which arises in the case of measurement saturation. The probability satisfies . In order to place an upper-bound on , we bound the error probabilities and separately.

#### Iii-A1 An Upper Bound on P(E1)

Consider the separation between the quadratic measurements of and , given by

 ψ = |b1−b2|=∣∣∣∣a⊤x1∣∣2−∣∣a⊤x2∣∣2∣∣ (11) = |a⊤(x1x⊤1−x2x⊤2)Va|.

An important property of is established in Proposition 2, eventually leading to a bound on .

###### Proposition 2

Let and

be two linearly independent vectors in

. The matrix

has two nonzero eigenvalues of equal magnitude

, such that , where is the coefficient of correlation.

Proof: Observe that has two nonzero eigenvalues, since, for any vector in the orthogonal complement of , we have . Therefore, the eigenvectors corresponding to the nonzero eigenvalues must be of the form , for scalars and that are not simultaneously equal to zero. If the corresponding eigenvalue is , then we have

 (x1x⊤1−x2x⊤2)(αx1+βx2)=ν1(αx1+βx2). (12)

Since and are linearly independent and of unit norm, by comparing terms in (12), we get that

 α+β(x⊤1x2)=ν1α,\,\,and\,\,α(x⊤2x1)+β=−ν1β.

Writing and substituting in , and using the fact that , we have

 β(1+ν1−ρ21−ν1)=0,

leading to , since .
Therefore, considering the separation in (11) and invoking the eigenvalue decomposition (EVD) , we get

 ψ = |a⊤Va|=∣∣a⊤RDR⊤a∣∣=|~a⊤D~a|, (13)

where . The top two diagonal entries of have identical magnitudes , and the remaining diagonal entries are zero. Taking this property into account, the separation may be bounded from below as

 ψ≥√1−ρ2|~a21−~a22|, (14)

where and are respectively the first and second entries of , which is a Gaussian random vector since is orthonormal and . Hence, the probability of can be upper-bounded as follows:

 P(E1) = P(ψ≤δ)(a)≤P(√1−ρ2|~a21−~a22|≤δ), ⇒P(E1) ≤ P(∣∣~a21−~a22∣∣<δ√1−ρ2), (15)

where the inequality (a) in (15) follows from the fact that the event implies that happens, as a consequence of (14). The random variables and are independent and follow the distribution (cf. Proposition 1

). The moment generating function (m.g.f.) of the

random variable is

 Mχ21(t)=(1−2t)−1/2,\,\,for\,\,t<12.

Let , which is the difference between two such independent random variables. The m.g.f. of turns out to be

 Mu(t)=(1−4t2)−12=⎛⎝1414−t2⎞⎠12,|t|<12. (16)

Comparing (16) with the m.g.f. of the variance-gamma distribution [45] with parameters , , , and , given by

 M(t)=eμot(α2−β2α2−(β+t)2)λ0, (17)

we obtain an equivalence of (17) and (16) with , , , and

. The probability density function (p.d.f.) corresponding to the m.g.f. in (

16) is given by

 fu(u)=12πK0(|u|2), (18)

where is the modified Bessel-function of the second kind and zeroth order. Since the p.d.f. in (18) is symmetric, the upper-bound on in (15) reduces to

 P(E1)≤2π∫δ√1−ρ20K0(u)dug(δ′), (19)

where . In order to obtain a more readily interpretable upper bound, we approximate the integral in (19) as . Figure 1(a) shows that is a reliable and accurate approximation to . Consequently,

 P(E1)≤1−exp(−1.6δ√1−ρ2). (20)

#### Iii-A2 Computing an Upper Bound on P(E2)

Since , which corresponds to both and going into saturation, where and are not independent. The event is clearly a subset of the event , and using Proposition 1, the probability of can be bounded as

 P(E2)≤P(b1≥τk−1)=1−γ(τk−12,12)√π. (21)

Finally, we combine the upper bounds on the error events and to obtain a bound on , which dictates the minimum number of measurements required to ensure that two linearly independent signals can be discerned from their quantized intensity measurements.

#### Iii-A3 An Upper Bound on Pe

Combining (20) and (21) results in an upper bound on , given by

 Pe≤Pmaxe=2−γ(τk−12,12)√π−exp(−1.6δ√1−ρ2), (22)

which, we recall, is the probability with which the measurements would be indistinguishable. The upper bound, in turn, must be less than unity for it to be meaningful. The quantizer parameters and may be chosen accordingly.

#### Iii-A4 A Numerical Example

Consider independent measurements corresponding to and . The overall probability that the measurements will be indistinguishable is upper-bounded as If we desire to keep this value below a certain , the minimal number of measurements would be . Consider the case where and a four-bit quantizer (), with equiprobable intervals, meaning that, , for . Such a quantizer ensures that roughly the same number of measurements fall in each interval, for large enough number of measurements . The illustrative values chosen for and are indicative of a fair degree of correlation between and and coarse quantization. The values of and for this particular choice turn out to be and . The corresponding upper bound on in (22) evaluates to . For independent measurements, we have . If we desire to achieve , we must have .
The minimum number of measurements required to maintain is shown in Figure 1(b) as a function of , for various values of . The value of is kept fixed at , which corresponds to , thereby guaranteeing that the probability of measurement saturation is no more than . For a fixed , the number of required measurements increases with increasing (coarse quantization). Moreover, for a given quantizer precision , one needs to collect more measurements in order to maintain distinguishability as the correlation increases. The upper bound on in (22) is not tight and the minimum number of measurements required, as given by the preceding analysis, is independent of the reconstruction algorithm, and is necessary, but not sufficient.

### Iii-B The Issue of Quantizer Design

As shown in Section III, a quadratic measurement of the form of a signal follows the c.d.f. , plotted in Figure 2(a). The quantizer may be designed using the Lloyd-Max algorithm [46], which jointly optimizes for the thresholds and the codebook, such that the quantization SNR, defined as

 SNRquant=∑mi=1b2i∑mi=1(bi−yi)2, (23)

is maximized, where and are as defined in (1) and (2), respectively. As emphasized in Section II, the proposed formalism based on consistent recovery is independent of the choice of encoding symbols. This aspect will also become clear as we develop the optimization framework for QPR in Section IV. Consequently, there is no guarantee that maximization of would lead to a superior phase retrieval performance. In the present context, it would be more appropriate to optimize the thresholds such that the expected estimation error, given by

 (τ∗1,⋯,τ∗k−1)=argminτ1,⋯,τk−1\,\,E∥^x(τ1,⋯,τk−1)−x∗∥22, (24)

where the expectation is calculated over the distribution of the sampling vectors , is minimized. Unfortunately, solving (24) iteratively or in closed-form is mathematically intractable. Therefore, we adopt an approach similar to that proposed by Zymnis et al. [34] in the context of quantized CS, wherein the quantizer is designed to have equiprobable intervals (along the lines of companding [47]). The design of such a quantizer is facilitated by the knowledge of the c.d.f. as illustrated in Figure 2(a) for levels. The thresholds are marked as with and . For the purpose of illustration, we show in Figure 2(b), as a function of the number of bits for equiprobable-interval quantizers. The encoding symbols for the levels are taken as the midpoints of the corresponding intervals, whereas the encoding symbol for the level is set to , where is as defined in (10). We observe from Figure 2(b) that tends to increase almost linearly as the number of bits increases, and attains a value of approximately dB when the number of quantization bits is . However, for coarse quantization, is about dB or lower, which is far too low for the existing PR algorithms, considering that one models the quantization noise as additive. This underscores the need for developing a quantization-aware PR algorithm.

### Iii-C Noise Robustness of Quantized Measurements

In practice, due to noise, the acquired measurements take the form

 (25)

where are independent and identically distributed (i.i.d.) additive noise samples. The quantization process is inherently noise-robust as long as the perturbation due to noise does not alter the output symbol. We shall illustrate this inherent robustness using Monte-Carlo simulations.
Let . We define the robustness factor corresponding to a representative quadratic measurement as

 pr=P(Q(b)=Q(b+ξ)),\,\,where\,\,b=∣∣a⊤x∗∣∣2, (26)

which measures the probability that the noise does not alter the quantized measurement. Recall that the intensity measurement follows the distribution. Since it is cumbersome to obtain an analytical expression for , for the purpose of illustration, we adopt a Monte-Carlo approach to estimate as

 ^pr=1NtrialNtrial∑t=11{Q(bt)=Q(bt+ξt)}, (27)

where denotes the indicator function, and and are drawn independently from the and distributions, respectively. The number of trials is taken as . The variation of as a function of , for different number of quantization levels , is shown in Figure 3. An equiprobable interval quantizer is considered. We observe from Figure 3 that the measurements get increasingly robust to noise as the quantization gets coarser. As one would expect, drops monotonically as the noise variance increases. For binary measurements (), the noise does not alter the measurements with probability approximately , even for relatively high noise variance of . The trade-off, however, is that such a coarse quantization will compromise distinguishability of the measurements.

## Iv The Optimization Framework and Reconstruction Algorithms for QPR

### Iv-a The QPR Optimization Framework

We combine the requirement of consistent recovery with the principle of lifting [26, 27] and formulate an appropriate cost function. The central idea behind lifting is to write the quadratic expression as

 ∣∣a⊤ix∣∣2=Tr(AiX),

where , , and denotes the trace operator. The two representations are equivalent, but the advantage of the lifted version is that it enables one to express the quadratic measurements in 1-D as a set of linear measurements in 2-D. We assume that measurements are encoded with the symbol and denote the sampling vectors in as , for . Since the matrix is a rank-1 and positive semi-definite (PSD) matrix by construction, it would be imperative to enforce these conditions. Effectively, we seek such that , , and

 τj−1≤Tr(AjiX)≤τj;\,\,i=1:mj,\,\,j=1:k, (28)

where . We incorporate the inequality constraints in (28) arising out of the consistency criterion in the optimization objective by using the one-sided quadratic function :

 f(u)={12u2,\,\,if\,\,u≤0,\,and0,\,\,otherwise.

Therefore, the QPR problem may be formulated as

 ^X=argminX⪰0\,F(X)\,\,\,subject to\,\,\,rank(X)=1, (29)

where the optimization objective is given by

 F(X)=k∑j=1mj∑i=1f(τj−Tr(AjiX))+f(%Tr(AjiX)−τj−1). (30)

Although the QPR objective function is convex in , the minimization posed in (29) is non-convex, because of the rank constraint.

### Iv-B Reconstruction Algorithms for QPR

We develop two projected gradient-based algorithms to solve (29), wherein one retains the best rank-1 approximation of the estimate obtained following a gradient-based update. The algorithms could be terminated whenever the measurement consistency requirement in (28) is met or when a maximum number of iterations have elapsed. The proposed algorithmic framework is amenable to accommodating constraints such as positivity, sparsity, etc., which are relevant in many practical imaging modalities. The sparsity prior may be incorporated by hard-thresholding the estimate obtained subsequent to applying the rank-1 constraint, akin to the approach we developed in [43].

#### Iv-B1 Projected Gradient-Descent (PGD) for Quantized PR

The first algorithm employs a simple projected gradient-descent technique, wherein one iteratively computes an update of the form

 Xt+1=Prank−1(Xt−ηtGt), (31)

starting with an initial estimate , where is the step-size parameter, and is the gradient matrix of evaluated at . We shall refer to (31) as the QPR update. The rank-1 projection operator applied on a symmetric matrix of size is defined as follows:

 Prank−1(Y)=max{λmax,0}v1v⊤1,

where is the largest eigenvalue of , having as the associated eigenvector. An estimate of the underlying ground-truth signal can be obtained as .

Cai et al. developed a similar singular-value thresholding algorithm, albeit with a soft-threshold, for solving the problem of low-rank matrix recovery from linear measurements

[48]. Calculating the gradient of requires the gradient of functions of the form , which is given by

 ∇h(X)=f′(Tr(AX)−τ)A⊤,

where, for any , is given by

 f′(u)={u,\,\,if\,\,u≤0,\,and0,\,\,otherwise.

To incorporate sparsity of , the QPR update is subjected to a hard-thresholding operation of the form

 xt+1←Ps(xt+1), (32)

where returns the best -sparse approximation of its argument and is obtained by picking the top entries in magnitude. We refer to (32) as the SQPR update. To determine , we adopt the exact line-search strategy by solving

 ηt=argminη>0\,\,F(Xt−ηGt), (33)

using a grid search over a chosen interval. The steps involved in QPR and SQPR are listed in Algorithm 1. The PGD algorithm for QPR possesses the descent property, that is, the updates generated by (31) satisfy . A proof of this claim is provided in Section IX-D of the supplementary material.

#### Iv-B2 Accelerated Projected Gradient-Descent (APGD)

Although the cost function in (30) is convex, the rank-1 projection step in (31) is not. In general, an acceleration of the PGD using Nesterov’s scheme [49] is not guaranteed in this setting. Nonetheless, motivated by the accelerated singular-value hard-thresholding strategy adopted in [50] for low-rank matrix completion problems, we go ahead with incorporating a momentum factor in the QPR and SQPR algorithms and investigate empirically if it would result in acceleration. It turns out, from the simulation results, that incorporating the momentum factor indeed results in accelerated convergence (Section V-A contains the simulation results). The steps of the QPR and SQPR algorithms with acceleration, referred to as QPR-A and SQPR-A, respectively, are summarized in Algorithm 2, where the step-sizes are chosen following (33). The updates for the momentum terms and in Algorithm 2 are based on the recommendations given in [27].

## V Numerical Experiments Without the Sparsity Constraint

In this section, we demonstrate the following via numerical simulations: (i) effect of Nesterov’s acceleration scheme; (ii) choice of quantizer design — Lloyd-Max quantizer (LMQ) versus equiprobable quantizer (EQ); and (iii) comparison of the proposed quantized PR algorithm with state-of-the-art PR techniques in the absence of any external additive noise, and without the assumption of sparsity. An assessment of the robustness to noise and reconstruction with the incorporation of a sparsity prior will be presented in Sections VI and VII, respectively. The state-of-the-art techniques used for comparison are PhaseLift, TWF, and AltMinPhase.
The PhaseLift approach is implemented using the PGD and the APGD algorithms, and referred to as PL and PL-A, respectively. The Matlab implementation of the TWF algorithm is taken from the authors’ website. The spectral initialization technique of [18] is employed to initialize PhaseLift, TWF, and AltMinPhase, wherein one sets to be equal to the eigenvector corresponding to the largest eigenvalue of the matrix , normalized to have unity -norm. The spectral initialization depends on the measurement vector , which is a function of the encoding symbols. The QPR-A algorithm, however, is initialized with an all-zero vector, thereby avoiding any dependence of the reconstructed signal on the choice of the codebook .

Experiments are conducted for the real signal model, where is drawn uniformly at random on and the measurement vectors , with . The number of measurements is taken as . The step-size parameter is chosen by an exhaustive search over the range , with a grid spacing of .

### V-a Effect of the Momentum Factor: QPR Versus QPR-A

Since the rank-1 projection operator in the update rule (31) is not convex, it is not obvious a priori that incorporating the momentum factor, explained in Section IV-B2, would necessarily lead to fast convergence or give any performance gains. A similar dilemma was encountered in the context of PhaseLift as well. Therefore, we compare the performances of QPR and QPR-A to determine which of them results in superior reconstruction, and also compare the results to PhaseLift with and without acceleration. We consider measurements quantized using levels, EQ for QPR and QPR-A, and LMQ for PL and PL-A. These are the optimal quantizer settings for the respective algorithms as will be demonstrated in the following subsection.

The average reconstruction SNRs and their standard deviations for QPR and QPR-A, calculated over

independent trials, are shown in Figure 4(a) with respect to iterations. The same metrics for PL and PL-A are shown in Figure 4(b). We observe from Figure 4(a) that QPR-A results in an improvement of approximately dB over QPR after iterations. The variation of the reconstruction SNR around its average value is also found to be slightly smaller for QPR-A. On the other hand, PL-A leads to faster convergence than PL, as can be inferred from Figure 4(b), although the final SNR after iterations settles to more or less the same value for both of them. This trend was found to be consistent for different values of and . The reconstructed signals obtained in a random trial using QPR-A and PL-A are compared against the ground-truth in Figures 4(c) and 4(d), respectively. QPR-A yields an improvement of approximately dB in the reconstruction SNR over PL-A. For further comparisons, we consider only QPR-A and PL-A owing to their superiority over QPR and PL, respectively.

### V-B Quantizer Design: Lloyd-Max Versus Equiprobable

In LMQ, the thresholds and the encoding symbols are jointly optimized to maximize the quantization SNR. On the other hand, the thresholds in EQ are chosen such that each interval has equal probability. The encoding symbols in EQ are taken as , for , and , where is as defined in (10). Since the consistency criterion is enforced in QPR-A, the specific choice of the codebook has no bearing on the performance of QPR algorithms, as explained in Section III-B. In other words, the estimated signals obtained using QPR-A corresponding to two different quantizers having the same intervals, but different codewords, would be the same. The variations of reconstruction SNR versus iterations for and , averaged over random trials, are shown in Figure 5. Comparing Figures 5(a) with 5(b) and 5(c) with 5(d), we observe that the performances of PL-A and TWF improve significantly when LMQ is used for measurement quantization; whereas the performance of AltMinPhase remains approximately the same under LMQ and EQ. The reconstruction SNR of QPR-A corresponding to LMQ initially increases with iterations, but drops as the number of iterations exceeds . However, when the EQ is used for measurement quantization, we observe that QPR-A leads to a steady increase of the reconstruction SNR as the iterations progress. The experiment indicates that the EQ is a better choice than LMQ for QPR-A, whereas for the competing algorithms, the LMQ is better. The superior performance of the competing algorithms with LMQ is not too surprising since they are not quantization-aware. Any quantization noise would only be treated as additive noise and their performance would be the best when the quantization noise variance is the least, which is what the LMQ guarantees. Therefore, in order to facilitate a fair comparison, in the sequel, we report the performance of QPR-A with the EQ and its competitors with the LMQ.

### V-C Comparison of QPR-A With the State-of-the-Art

The comparative performances of QPR-A and the competing algorithms in terms of the reconstruction SNR and measurement consistency are shown in Figures 6 and 7, respectively, corresponding to various quantization levels . We observe from Figure 6 that QPR-A ultimately results in superior reconstruction SNR despite having a suboptimal initialization, and this trend is consistent for all values of . We observe from Figure 6 that the improvement in reconstruction SNR obtained using QPR-A over TWF, the best performing competing technique, is approximately dB for , and reduces to nearly dB as increases to . Moreover, unlike the competing techniques, the reconstruction SNR of QPR-A does not seem to saturate fast. Naturally, the reconstruction SNR increases with for all algorithms because the quantization gets finer. As far as measurement consistency is concerned, we observe from Figure 7 that QPR-A steadily improves with iterations and eventually performs on par with its competitors. In other words, as the iterations progress, the reconstruction produced by QPR-A provides an accurate explanation of the acquired measurements.

To summarize, the accelerated algorithms QPR-A and PL-A offer superior reconstruction performance than QPR and PL, respectively. An EQ is a better choice for QPR-A and LMQ for the competing algorithms. In the absence of noise, the QPR-A technique is at least dB better than the best performing PR technique, namely, the TWF. Before concluding this section, we show an application of QPR-A to the reconstruction problem in frequency-domain optical coherence tomography (FDOCT). An example of natural image reconstruction using QPR-A is shown in Section IX-A of the supplementary document.

### V-D Application of QPR-A to Frequency-Domain Optical Coherence Tomography (FDOCT)

We consider signal reconstruction in FDOCT, a non-invasive imaging technique used for obtaining structural details of biological specimens. A detailed description of the acquisition setup and the signal model in FDOCT that is relevant to the present discussion can be found in [42]. The interference pattern formed by the reflected signals from the object and reference arms approximates the Fourier transform of the object wave and is recorded by the spectrometer. The key challenge is to reconstruct the reflected wave from the object arm, which carries structural information about the specimen, from the intensity recordings of the spectrometer. Since the reflected wave exhibits a strong peak only when there is a significant change of refractive index in the specimen, the assumption of sparsity is appropriate in this context. However, the QPR-A algorithm leads to a fairly accurate reconstruction of the tomograms even without the sparsity assumption, as we shall show next.
FDOCT reconstruction of the glass specimen333The FDOCT data is the courtesy of Prof. R. A. Leitgeb, Medical University of Vienna, Austria. produced by the QPR-A algorithm is shown in Figure 8 corresponding to three different quantization levels, namely , , and . Considering the max- reconstruction [43] as the ground-truth, finite-precision measurements of the form (2) are collected corresponding to every scan-line. We consider a downsampled version (by a factor of four) of the back-scattered wave along each scan-line, leading to a signal dimension of , for the purpose of illustration. We observe that the QPR-A algorithm can reconstruct the back-scattered wave and recover the structural details in the specimen reliably. Imposing the sparsity constraint iteratively (more on the effect of sparsity in Section VII), with a sparsity level , helps eliminate the background noise significantly, as one can observe from Figures 8(b), 8(d), and 8(f). These results also show that increasing the quantizer precision leads to a higher accuracy in tomogram reconstruction.

## Vi Noise Robustness: MSE vis-à-vis the CRB

We now consider the effect of additive white Gaussian noise, prior to quantization, giving rise to quantized measurements

 (34)

where the noise samples are drawn independently from the distribution. For illustration, the ground-truth signal is taken as a sum of two sinusoids:

 x∗ℓ=Cx∗[1.5sin(4π(ℓ−1)n)+2.5cos(14π(ℓ−1)n)], (35)

where and the constant is chosen such that . The derivation of the CRB for a -level quantizer is given in Section IX-C of the supplementary material. Since the CRB is used as a theoretical benchmark, we use the reconstruction MSE as a performance metric. The reconstruction MSEs of different algorithms are compared against the CRB for three different quantization levels, namely and . The MSEs are computed according to (6), and averaged over noise realizations corresponding to a fixed measurement matrix , whose entries are i.i.d. and follow the distribution. Since the measurement matrix is random, we need one more level of averaging, which is performed over different measurement matrices. The results are shown in Figure 9 as a function of the input SNR defined as . The legends CRB-EQ and CRB-LMQ in Figure 9 denote the Cramér-Rao bounds corresponding to quantizers EQ and LMQ, respectively. The Fisher information matrix corresponding to the LMQ was found to be nearly rank-deficient for and therefore we omitted CRB-LMQ in Figure 9(a).
We observe that QPR-A attains reconstruction MSEs within dB of the corresponding CRB, whereas the other algorithms do not follow the CRB with increasing input SNR, especially for coarse quantization (). For input SNR greater than dB, QPR-A has a reconstruction MSE closer to the CRB than other techniques. At low input SNRs (below dB), the additive noise leads to a violation of the consistency condition in (28), which forms the basis of QPR-A. As a result, the QPR-A algorithm does not lead to a significant improvement over TWF and PL-A. We believe that a relaxed consistency condition in (28) to account for noise might result in more accurate recovery at low input SNRs. This aspect requires a separate investigation.

## Vii Numerical Experiments on QPR With Sparsity

We now take into account sparsity of the signal and analyze the reconstruction capability of SQPR-A vis-à-vis the state-of-the-art algorithms for sparse PR. As the experimental results show, when the ground-truth is indeed sparse, incorporating that prior actually helps improve the reconstruction performance. This point is illustrated by comparing the reconstructed signals using SQPR-A and QPR-A for binary quantization () in Figure 10. We observe that imposition of the sparsity prior leads to an improvement of about dB in the reconstruction SNR.

### Vii-a Impact of the Sparsity Prior: SQPR-A versus QPR-A

We compare SQPR-A with QPR-A for different values of and the relative sparsity level , the fraction of nonzero entries in the ground-truth. The results are averaged over trials and presented in Figures 11(a) and 11(b). The figures show a distinct improvement coming from the sparsity prior. For example, for sparsity level , and , SQPR-A is about dB better than QPR-A. The measurement consistency index is also higher and exhibits a faster convergence with iterations. As the underlying signal gets denser, the improvement one can achieve by enforcing sparsity diminishes (cf. Figures 11(c) and 11(d)).

### Vii-B SQPR-A Versus State-of-the-Art Sparse PR Algorithms

We now compare SQPR-A with three state-of-the-art sparse PR techniques: (i) CPRL [28]; (ii) GESPAR [17]; and (iii) SparseAltMinPhase [18], the sparse counterpart of AltMinPhase. Recall that in CPRL, sparsity is enforced by incorporating an penalty term, thereby leading to the following semi-definite program:

 argminX⪰0\,Tr(X)+λ1∥X∥1\,\,s.t.\,\,m∑i=1(%Tr