\section{Introduction}

Random projections are at the heart of many algorithms in machine learning, signal processing and numerical linear algebra. Recent developments ranging from classification with random features \cite{rahimi2008random}, kernel approximation \cite{yang2017randomized} and sketching for matrix optimization \cite{tropp2017practical, yurtsever2017sketchy}, to sublinear-complexity transforms \cite{yu2016orthogonal} and randomized linear algebra are all enabled by random projections. Computing random projections for realistic signals such as images, videos, and modern big data streams is computation- and memory-intensive. Thus, from a practical point of view, any increase in the size and speed at which one can do the required processing is highly desirable.

This fact has motivated work on using dedicated hardware based on physics rather than traditional CPU and GPU computation to obtain random projections. A notable example is \finalrev{the} scattering of light in random media (Figure \ref{fig:opu_setup} (left)) with an optical processing unit (OPU). The OPU enables rapid (20 kHz) projections of high-dimensional data such as images, with input dimension scaling up to one million and output dimension also in the million range. It works by ``imprinting'' the input data $\vxi \in \R^N$ onto a coherent light beam using a digital micro-mirror device (DMD) and shining the modulated light through a multiple scattering medium such as titanium dioxide white paint. The scattered lightfield in the sensor plane can then be written as
\begin{align*}
    \vy = \mA \vxi
\end{align*}
where $\mA \in \C^{M \times N}$ is the transmission matrix of the random medium with desirable properties.

One of the major challenges associated with this approach is that $\mA$ is in general \textit{unknown}. Though it could in principle be learned via calibration \cite{dremeau2015reference}, such a procedure is slow and inconvenient, especially at high resolution. On the other hand, the system can be designed so that the distribution of $\mA$ is approximately iid standard complex Gaussian. Luckily, this fact alone is sufficient for many algorithms and the actual values of $\mA$ are not required.

Another challenge is that common light sensors are only sensitive to intensity, so we can only measure the intensity of scattered light, $\abs{\vy}^2$, where $\abs{\cdot}$ is the elementwise absolute value. The phase information is thus lost. While the use of interferometric measurements with a reference could enable estimating the phase, the practical setup is more complex,  sensitive, and it does not share the convenience and simplicity of the one illustrated in Figure~\ref{fig:opu_setup} (left).

This motivates us to consider the \textit{measurement phase retrieval (MPR) problem}. The MPR sensor data is modeled as
\begin{align}
    \vb = \abs{\vy}^2 + \veta = \abs{\mA \vxi}^2 + \veta, \label{eq:phase_problem}
\end{align}
where $\vb \in \R^M$, $\vxi \in \R^N$, $\mA \in \C^{M \times N}$, $\vy \in \C^M$, and $\veta \in \R^M$ is noise. The goal is to recover the phase of each complex-valued element of $\vy$, $y_i$ for $1 \leq i \leq M$, from its magnitude measurements $\vb$ when $\vxi$ is \textit{known} and the entries of $\mA$ are \textit{unknown}. The classical phase retrieval problem which has received much attention over the last decade \cite{netrapalli2013phase,candes2015phase} has the same quadratic form as \eqref{eq:phase_problem} but with a known $\mA$ and the task being to recover $\vxi$ instead of $\vy$. While at a glance it might seem that not knowing $\mA$ precludes computing the phase of $\mA \vxi$, we show in this paper that it is in fact possible via an exercise in distance geometry.

\begin{figure}
\centering
\def4.2cm{7.3cm}
\fontsize{8}{8}\selectfont
\import{figures/}{explainer_with_DMD.pdf_tex}
\hspace{5em}
\def4.2cm{4.2cm}
\import{figures/}{explainer_plane.pdf_tex}
\caption{Left: The optical processing unit (OPU) is an example application of where the MPR problem appears. A coherent laser beam spatially encodes a signal $\left(\vx_q - \vx_r \right)$ via a digital micro-mirror device (DMD) which is then shined through \rev{a random medium}. A camera measures the squared magnitude of the scattered light which is equivalent to the Euclidean distance between complex numbers $y_{q} \in \C$ and $y_{r} \in \C$. Furthermore the camera takes quantized measurements;
Right: $y_q$ and $y_r$ are points on the two-dimensional complex plane. We can measure the squared Euclidean distance between points and use these distances to localize points on the complex plane and obtain their phase. Note that transformations such as rotations and reflections do not change the distances.}
\label{fig:opu_setup}
\end{figure} 

The noise $\veta$ is primarily due to quantization because standard camera sensors measure low precision values, 8-bit in our case (integers between 0 and 255 inclusive). Furthermore, cameras may perform poorly at low intensities. This is another data-dependent noise source which is modelled in \eqref{eq:phase_problem_with_mask} by a binary mask vector $\vw \in \R^M$ which is zero when the intensity is below some threshold and one otherwise; $\odot$ denotes the elementwise product.
\begin{align}
    \vb = \vw \odot \left( \abs{\vy}^2 + \veta \right) = \vw \odot \left( \abs{\mA \vxi}^2 + \veta \right) \label{eq:phase_problem_with_mask}
\end{align}

The distribution of $\mA$ follows from the properties of random scattering media \cite{liutkus2014imaging,dremeau2015reference}. It has iid standard complex Gaussian entries, $a_{mn} \sim \mathcal{N}(0,1) + j\mathcal{N}(0,1)$ for all $1 \leq m, n \leq M, N$. 

The usefulness of phase is obvious. While in some applications having only the magnitude of the random projection is enough (see \cite{saade2016random} for an example related to elliptic kernels), most applications require the phase. For example, with the phase one can implement a more diverse range of kernels as well as randomized linear algebra routines like randomized singular value decomposition (SVD). We report the results of the latter on real hardware in Section \ref{sec:experiments}.

\paragraph{Our contributions.}

We develop an algorithm based on distance geometry to solve the MPR problem \eqref{eq:phase_problem}. We exploit the fact that we control the input to the system, which allows us to mix $\vxi$ with arbitrary reference inputs. By interpreting each pixel value as a point in the complex plane, this leads to a formulation of the MPR problem as a pure distance geometry problem (see Section \ref{sub:mprdist} and Figure \ref{fig:opu_setup} (right)). With enough pairwise distances (corresponding to reference signals) we can localize the points on the complex plane via a variant of multidimensional scaling (MDS) \cite{torgerson1952multidimensional,dokmanic2015euclidean}, and thus compute the missing phase.

As we demonstrate, the proposed algorithm not only accurately recovers the phase, but also improves the number of useful bits of the magnitude information thanks to the multiple views. Established Euclidean distance geometry bounds imply that even with many distances below the sensitivity threshold and coarse quantization, the proposed algorithm allows for accurate recovery. This fact, which we verify experimentally, could have bearing on the design of future random projectors by navigating the tradeoff between physics and computation.

\subsection{Related work}

The classical phase retrieval problem looks at the case where $\mA$ is known and $\vxi$ has to be recovered from $\vb$ in \eqref{eq:phase_problem} \cite{fienup1982phase, shechtman2015phase, jaganathan2015phase}. A modified version of the classical problem known as holographic phase retrieval is related to our approach: a known reference signal is concatenated with $\vxi$ to facilitate the phase estimation \cite{barmherzig2019holographic}. Interference with known references for classical phase retrieval has also been studied for known (Fourier) operators \cite{beinert2017one, kim1990phase} . 

\rev{An optical random projection setup similar to the one we consider has been used for kernel-based classification \cite{saade2016random}, albeit using only magnitudes. A  phaseless approach to classification with the measured magnitudes fed into a convolutional neural network was reported by Satat et al. \cite{satat2017object}.

An alternative to obtaining the measurement phase is to measure, or calibrate, the unknown transmission matrix $\mA$. This has been attempted in compressive imaging applications but the process is impractical at even moderate pixel counts \cite{dremeau2015reference, liutkus2014imaging}. Estimating $\mA$ can take days and even the latest GPU-accelerated methods take hours for moderately sized $\mA$ \cite{sharma2019inverse}. Other approaches forego calibration and use the measured magnitudes to  \textit{learn} an inverse map of $\vx \mapsto |\mA \vx|^2$ for use with the magnitude measurements \cite{horisaki2016learning}. 

Leaving hardware approaches aside, there have been multiple algorithmic efforts to improve the speed of random projections \cite{le2013fastfood,yang2017randomized} for machine learning and signal processing tasks. Still, efficiently handling high-dimensional input  remains a formidable challenge.}






\section{The measurement phase retrieval problem}

We will denote the signal of interest by $\vxi \in \R^N$, and the $K$ reference \emph{anchor} signals by $\vr_k \in \R^N$ for $1 \leq k \leq K$. To present the full algorithm we will need to use multiple signals of interest which we will then denote $\vxi_1, \ldots, \vxi_S$; each $\vxi_s$ is called a frame. We set the last, $K$th anchor to be the origin, $\vr_K = \vzero$. We ascribe $\vxi$ and the anchors to the columns of the matrix $\mX \in \R^{N \times Q}$, so that $\mX = \left[\vxi, \vr_1, \vr_2, \cdots, \vr_K \right]$ and let $Q=K+1$. The $q$th column of $\mX$ is denoted $\vx_q$. For any $1 \leq q, r \leq Q$, we let $\vy_q = \mA \vx_q$ and $\vy_{qr} := \mA(\vx_q - \vx_r)$, with $y_{qr, m}$ being its $m$th entry. Finally, the $m$th row of $\mA$ will be denoted by $\va^{m}$ so that $y_{qr,m} = \inprod{\va^m, \vx_q-\vx_r}$. 


\subsection{Problem statement and recovery up to a reference phase and conjugation}

Since we do not know $\mA$, it is clear that recovering the absolute phase of $\mA \vxi$ is impossible. On the other hand, many algorithms do not require any knowledge of $\mA$ except that it is iid standard complex Gaussian, and that it does not change throughout the computations. 

Let $\mathsf{R}$ be an operator which adds a constant phase to each row of its argument (multiplies it by $\diag (e^{j \phi_1}, \ldots, e^{j \phi_m})$ for some $\phi_1, \ldots, \phi_m$) and conjugates a subset of its rows.
Since a standard complex Gaussian is circularly symmetric, $\mathsf{R} (\mA)$ has the same distribution as $\mA$. 
\rev{Therefore, since we do not know $\mA$, it does not matter whether \finalrev{we} work with $\mA$ itself or with $\mathsf{R}(\mA)$ for some possibly unknown $\mathsf{R}$. As long as the same \textit{effective} $\mathsf{R}$ is used for all inputs during algorithm operation, the relative phases between the frames will be the same whether we use $\mathsf{R}(\mA)$ or $\mA$.\footnote{Up to a sign.}}



\begin{problem} \label{prob:main_problem}
Given a collection of input frames $\vxi_1, \ldots, \vxi_S$ to be randomly projected and a device illustrated in Figure \ref{fig:opu_setup} (left) with an unknown transmission matrix $\mA \in \C^{M \times N}$ and a $b$-bit camera,  compute the estimates of projections $\hat{\vy}_1, \ldots, \hat{\vy}_S$ up to a global  row-wise phase and conjugation; that is, so that there exists some $\mathsf{R}$  such that $\hat{\vy}_s \approx \mathsf{R} (\vy_s)$ for all $1 \leq s \leq S$.
\end{problem}

\subsection{MPR as a distance geometry problem}
\label{sub:mprdist}

Since the rows of $\mA$ are statistically independent, we can explain our algorithm for a single row and then repeat the same steps for the remaining rows. We will therefore omit the row subscript/superscript $m$ except where explicitly necessary.

Instead of randomly projecting $\vxi$ and measuring the corresponding projection magnitude $|\mA \vxi|^2$, consider randomly projecting the difference between $\vxi$ and some reference vector, or more generally a difference between two columns in $\mX$, thus measuring $|\inprod{\va, \vx_q-\vx_r}|^2 = \abs{y_{q} - y_{r}}^2$. Interpreting $y_{q}$ and $y_{r}$ as points in the complex plane, we see that the camera sensor measures exactly the squared Euclidean distance between them. Since we control the input to the OPU, we can indeed set it to $\vx_q - \vx_r$ and measure $\abs{y_{q} - y_{r}}^2$ for all $1 \leq q,r \leq Q$.

This is the key point: as we can measure pairwise distances between a collection of two-dimensional vectors in the two-dimensional complex plane, we can use established distance geometry algorithms such as multidimensional scaling (MDS) to localize points and get their phase. This is illustrated in Figure \ref{fig:opu_setup} (right). The same figure also illustrates the well known fact that rigid transformations of a point set cannot be recovered from distance data. We need to worry about \rev{three} things: translations, \rev{reflections} and rotations.

The translation ambiguity can be easily dealt with if one notes that for any column $\vx_q$ of $\mX$, $|y_q| = |\inprod{\va, \vx_q}|$ gives us the distance of $y_q$ to the origin which is a fixed point, ultimately resolving the translation ambiguity. There is, however, no similar simple way to do away with the rotation \rev{and reflection} ambiguity, so it might seem that there is no way to uniquely determine the phase of $\inprod{\va, \vxi}$. This is where the discussion from the preceding subsection comes to the rescue. Since $\mathsf{R}$ is arbitrary, as long as \rev{it is kept} fixed for all the frames, we can arbitrarily set the orientation of any given frame and use \rev{it as a reference}, making sure that the relative phases are computed correctly.

\subsection{Proposed algorithm} \label{sec:algorithm}

As defined previously, the columns of $\mX \in \R^{N \times Q}$  list the signal of interest and the anchors. Recall that all the entries of $\mX$ are known. Using the OPU, we can compute a noisy (quantized) version of
\begin{align}
    \abs{y_{qr}}^2 = \abs{\inprod{\va, \vx_q - \vx_r}}^2 
    = \abs{y_{q} - y_{r}}^2,
\end{align}
for all $(q, r)$, which gives us $Q(Q-1)/2$ squared Euclidean distances between points $\set{y_q \in \C}_{q = 1}^Q$ on the complex plane. These distances can be  used to populate a Euclidean (squared) distance matrix $\mD \in \R^{Q \times Q}$ as
\(
\mD = (d_{qr}^2)_{q,r=1}^Q = (\abs{y_{qr}}^2)_{q,r=1}^Q,
\) 
which we will use to localize \rev{all complex points} $y_{q}$.

We start by defining the matrix of \rev{all the complex} points in $\R^2$ \rev{which we want to recover} as
\[
  \mat{\Upsilon} =
  \begin{bmatrix}
    \mathrm{Re}(y_1) & \mathrm{Re}(y_2) & \cdots & \mathrm{Re}(y_Q)\\
    \mathrm{Im}(y_1) & \mathrm{Im}(y_2) & \cdots & \mathrm{Im}(y_Q)
  \end{bmatrix} \in \R^{2 \times Q}.
\]
Denoting the $q$th column of $\mat{\Upsilon}$ by $\vec{\upsilon}_q$, we have $d_{qr}^2 = \norm{\vec{\upsilon}_q - \vec{\upsilon}_r}_2^2 = \vec{\upsilon}_q^\T \vec{\upsilon}_q - 2 \vec{\upsilon}_q^\T \vec{\upsilon}_r + \vec{\upsilon}_r^\T \vec{\upsilon}_r$ so that
\begin{align}
    \mD = \diag\left(\mG \right) \vec{1}_Q^\T  - 2\mG + \vec{1}_Q \diag\left(\mG \right)^\T 
    =: \mathcal{K}\left(\mG \right), \label{eq:distance_to_points}
\end{align}
where $\diag(\mG) \in \R^Q$ is the \rev{column} vector of the diagonal entries in \rev{the Gram matrix} $\mG := \mat{\Upsilon}^\T \mat{\Upsilon} \in \R^{Q \times Q}$ and $\vone_Q \in \R^Q$ is the \rev{column} vector of $Q$ ones. \rev{This establishes a relationship between the measured distances in $\mD$ and the locations of the complex points in $\R^2$ which we seek.} We denote by $\mJ$ the geometric centering matrix, $\mJ := \mI - \tfrac{1}{Q} \vone_Q \vone_Q^\T$ so that 
\begin{align}
  \wh{\mG} = -\tfrac{1}{2} \mJ \mD \mJ = \mJ \mG \mJ = (\mat{\Upsilon} \mJ)^\T (\mat{\Upsilon} \mJ) \label{eq:centered_gram}
\end{align}
is the Gram matrix of the centered point set in \rev{terms of} $\mat{\Upsilon}$. \rev{$\wh{\mG}$ and $\mJ$ are know as the Gram matrix of the centered point set and the geometric centering matrix because $\mat{\Upsilon}\mJ$ is the points in $\mat{\Upsilon}$ with their mean subtracted.} An estimate $\wh{\mat{\Upsilon}}$ of the centered point set\rev{, $\mat{\Upsilon}\mJ$,} is then obtained by eigendecomposition as $\wh{\mG} = \mV \diag(\lambda_1, \ldots, \lambda_Q) \mV^\T$ and taking $\wh{\mat{\Upsilon}} = [\sqrt{\lambda_1} \vv_1, \sqrt{\lambda_2} \vv_2]^\T$ \rev{where $\vv_1$ and $\vv_2$ are the first and second \finalrev{columns} of $\mV$} and assuming that the eigenvalue sequence is nonincreasing. This process is the classical MDS algorithm \cite{torgerson1952multidimensional,dokmanic2015euclidean}. Finally, the phases can be calculated via a four-quadrant inverse tangent, $\phi(y_q) = \mathrm{arctan}(\upsilon_{q2}, \upsilon_{q1})$.

\paragraph{Procrustes analysis.} 

A\rev{s we recovered a centered point set via MDS with a geometric centering matrix $\mJ$,} the point set will have its centroid at the origin. This is a consequence of the used algorithm, and not the ``true'' origin. As described above, we know that $|y_q|^2$ define\rev{s} squared distances to the origin and $y_Q = \inprod{\va, \vx_Q} = 0 + 0j$ \rev{(as $\vx_Q$ was set to the origin)}, meaning that we can \rev{correctly center the recovered points} by translating the point set\rev{, $\wh{\mat{\Upsilon}}$,} by $-\vec{\upsilon}_Q$. 

The correct absolute rotation \rev{and reflection} cannot be recovered. However, since we only care about working with some effective $\mathsf{R}(\mA)$ with the correct distribution, we only need to ensure that the relative phases between the frames are correct. We can thus designate the first frame as the reference frame and \rev{set} the rotation (which directly corresponds to the phase) \rev{and reflection (corresponding to conjugation)} arbitrarily. Once \rev{these are} chosen, the anchors $\vr_1, \ldots, \vr_K$ are fixed, which in turn fixes the phasing--conjugation operator $\mathrm{R}$.

Since $\mA$ is unknown, $\mathsf{R}$ is also unknown, but fixed anchors allow us to compute the correct relative phase with respect to $\mathsf{R}(\mA)$ for the subsequent \rev{frames}. Namely, upon receiving a new input $\vxi_s$ to be randomly projected, we now localize it with respect to a fixed set of anchors. This is achieved by Procrustes analysis. Denoting by $\wt{\mat{\Upsilon}}_1$ our reference estimate of the anchor positions in frame 1 (columns $2, \ldots, Q$ of $\wh{\mat{\Upsilon}}$ \rev{above which was recovered from $\wh{\mG}$ in \eqref{eq:centered_gram}}), and by $\wt{\mat{\Upsilon}}_s$ the MDS estimate of anchor positions in frame $s$, adequately \rev{centered}. Let $\wt{\mat{\Upsilon}}_s \wt{\mat{\Upsilon}}_1^\T = \mU \mSigma \mV^\T$ be the singular value decomposition of $\wt{\mat{\Upsilon}}_s \wt{\mat{\Upsilon}}_1^\T$. The optimal \rev{transformation} matrix in the least squares sense is then $\mR = \mV \mU^\T$ so that $\mR \wt{\mat{\Upsilon}}_s \approx \wt{\mat{\Upsilon}}_1$ \cite{schoenemann1964solution}.

Finally, we note that with a good estimate of the anchors, one can imagine not relocalizing them in every frame. The localization problem for $\vxi$ then boils down to multilateration, cf. Section \ref{sec:sr_ls} in the supplementary material.

\subsection{Sensitivity threshold and missing measurements} \label{sub:theory}

As we further elaborate in Section \ref{sec:practical_considerations} of the supplementary material, in practice some measurements fall below the sensitivity threshold of the camera and produce spurious values.  A nice benefit of multiple ``views'' of $\vxi$ via its interaction with reference signals is that we can ignore those measurements. This introduces missing values in $\mD$ which can be modeled via a binary mask matrix $\mW$. The recovery problem can be modeled as estimating $\mat{\Upsilon}$ from 
\(
    \mW \odot (\mD + \mE) \label{eq:edm_model}
\)
where $\mW \in \R^{N \times N}$ contains zeros for the entries which fall below some prescribed threshold, and ones otherwise.

We can predict the performance of the proposed method when modeling the entries of $\mW$ as iid Bernoulli random variables with parameter $p$, where $1 - p$ is the probability that an entry falls below the sensitivity threshold and $\mE$ as uniform quantization noise distributed as $\mathcal{U}\left(-\frac{\kappa}{2(2^b-1)}, \frac{\kappa}{2(2^b-1)} \right)$, where $b$ is the number of bits, and $\kappa$ an upper bound on the entries of $\mD$ (in our case $2^8 - 1 = 255$).

Adapting existing results on the performance of multidimensional scaling \cite{zhang2019localization} (by noting that $\mE$ is sub-Gaussian), we can get the following scaling of the distance recovery error with the number of anchors $K$ (for $K$ sufficiently large),
\begin{equation}
    \frac{1}{K} \expect{\norm{\hat{\mD} - \mD}_F} \lesssim \frac{\kappa}{\sqrt{pK}},
\end{equation}
where $\lesssim$ denotes inequality up to a constant which depends on the number of bits $b$, the sub-Gaussian norm of the entries in $\mE$, and the dimension of the ambient space (here $\R^2$). An important implication is that even for coarse quantization (small $b$) and for a large fraction of entries below the sensitivity threshold (small $p$), we can achieve arbitrarily small amplitude and phase errors \textit{per point} by increasing the number of reference signals $K$.

\paragraph{Refinement with gradient descent.} 

The output of the classical MDS method described above can be further refined via a local search. A standard differentiable objective called the squared stress is defined as follows,
\begin{align}
    \min_{\mZ} ~ f\left(\mat{\Upsilon} \right) = \min_{\mZ} ~ \norm{\mW \odot \left( \mD - \mathcal{K}\left(\mZ^\T \mZ \right) \right)}_F^2,    \label{eq:edm_nonconvex}
\end{align}
where $\mathcal{K}(\cdot)$ is as defined in \eqref{eq:distance_to_points} and $\mZ \in \R^{2 \times Q}$ is the point matrix induced by row $m$ of $\mA$. In our experiments we report the result of refining the classical MDS results via gradient descent on \eqref{eq:edm_nonconvex}.

Note that the optimization \eqref{eq:edm_nonconvex} is nonconvex. The complete procedure is thus analogous to the usual approach to nonconvex phase retrieval by spectral initialization followed by gradient descent \cite{netrapalli2013phase, candes2015phase}. Algorithm \ref{algo:mpr_algorithm} summarizes our proposed
method. 

\begin{algorithm}[t]
\caption{MPR algorithm for $S$ frames.}
\begin{algorithmic}[1]
\Require Squared distances $\big[\abs{y_{jQ, m} - y_{lQ, m}}^2\big]_s$ for all $1 \leq j,l \leq Q$ for frames $1 \leq s \leq S$ and rows $1 \leq m \leq M$; $[\,\cdot\,]_s$ denotes frame $s$
\Ensure $\mY \in \C^{M \times S}$ containing all localized points such that $\vy_s = \mathsf{R}(\mA)\vxi_s$ for some fixed $\mathsf{R}$.
\State $\mY \gets \vec{0}_{M \times S}$ \Comment{Initialize $\mY$}
\State $m \gets 1$
\While{$m \leq M$} \Comment{Solve each row separately}
\State Populate all frame $s = 1$ distances into distance matrix $\mD$ \Comment{$\mD \in \R^{Q \times Q}$}
\State $[\mat{\Upsilon}]_1 \gets$ \texttt{MDS}$(\mD)$ \Comment{$[\mat{\Upsilon}]_1 \in \R^{2 \times Q}$}
\State $[\mat{\Upsilon}]_1 \gets$ \texttt{GradientDescent}$(\mD, [\mat{\Upsilon}]_1)$ \label{algo:gd1}
\State $[\mat{\Upsilon}]_1 \gets [\mat{\Upsilon}]_1 - [\vec{\upsilon}_Q]_1 \vec{1}^\T$ \Comment{Translate to align with origin}
\State $s \gets 2$
\While{$s \leq S$}
\State Populate all frame $s$ distances into distance matrix $\mD$
\State $[\mat{\Upsilon}]_s \gets$ \texttt{MDS}$(\mD)$
\State $[\mat{\Upsilon}]_s \gets$ \texttt{GradientDescent}$(\mD, [\mat{\Upsilon}]_s)$ \label{algo:gd2}
\State $[\mat{\Upsilon}]_s \gets [\mat{\Upsilon}]_s - [\vec{\upsilon}_Q]_s \vec{1}^\T$
\State $\mR \gets $ \texttt{Procrustes}$([\vec{\upsilon}_2, \ldots, \vec{\upsilon}_Q]_1, [\vec{\upsilon}_2, \ldots, \vec{\upsilon}_Q]_s)$ \Comment{$\mR$ aligns frames 1 and $s$ anchors}
\State $[\mat{\Upsilon}]_s \gets $ \texttt{Align}$([\mat{\Upsilon}]_s, \mR, [\vec{\upsilon}_2, \ldots, \vec{\upsilon}_Q]_1)$ \Comment{Align anchors}
\State $s \gets s+1$
\EndWhile
\State $\mU \gets \big[[\vec{\upsilon}_1]_1, [\vec{\upsilon}_1]_2, \ldots , [\vec{\upsilon}_1]_S \big]$ \Comment{$\mU \in \R^{2 \times S}$}
\State $\vy^m \gets \vu^1 + j \vu^2$ \Comment{Multiply second row of $\mU$ with $j$ and add to first row}
\State $m \gets m+1$
\EndWhile
\end{algorithmic}
\label{algo:mpr_algorithm}
\end{algorithm}

\section{Experimental verification and application}

We test the proposed MPR algorithm via simulations and experiments on a real OPU. For hardware experiments, we use a scikit-learn interface to a publicly available cloud-based  OPU
\footnote{\url{https://www.lighton.ai/lighton-cloud/}. \\
Reproducible code available at \url{https://github.com/swing-research/opu_phase} under the MIT License.}

\paragraph{Evaluation metrics.}

The main challenge is to evaluate the performance without knowing the transmission matrix $\mA$. To this end, we propose to use the \textit{linearity error}. The rationale behind this metric is that with the phase correctly recovered, the end-to-end system should be linear. That is, if we recover $\vy$ and $\vz$ from $\abs{\vy}^2 = \abs{\mA\vxi_1}^2$ and $\abs{\vz}^2 = \abs{\mA\vxi_2}^2$, then we should get $(\vy + \vz)$ when applying the method to $\abs{\vv}^2 = \abs{\mA(\vxi_1 + \vxi_2)}^2$. With this notation, the relative linearity error is defined as
\begin{equation}
    \label{eq:linerr}
    \text{linearity error} = \frac{1}{M} \sum_{m=1}^M \frac {\abs{(y_m + z_m) - v_m}}{\abs{v_m}}.
\end{equation}

The second metric we use is the number of ``good'' or correct bits. This metric can only be evaluated in simulation since it requires the knowledge of the ground truth measurements. Letting $\abs{y}^2 = \abs{\inprod{\va, \vxi}}^2$ and $\hat{y}$ be our estimate of $
y$, the number of good bits is defined as 
\[
\text{good bits} = - \tfrac{20}{6.02}  \log\left({|  \abs{y}^2 - \abs{\hat{y}}^2|} \, / \, {\abs{y}^2} \right).
\]
It is proportional to the  signal-to-quantization-noise ratio if the distances uniformly cover all quantization levels.\footnote{Note that the quantity registered by the camera is actually the squared magnitude, hence the factor 20.}

\subsection{Experiments} \label{sec:experiments}

In all simulations, intensity measurements are quantized to 8 bits and all signals and references are iid standard (complex) Gaussian random vectors. 

We first test the phase recovery performance by evaluating the linearity error. In simulation, we draw random frames $\vxi_1$, $\vxi_2$,  and $\mA \in \C^{100 \times 64^2}$. We apply Algorithm \ref{algo:mpr_algorithm} to  $|\mA \vxi_1|^2$, $|\mA \vxi_2|^2$ and $|\mA(\vxi_1 + \vxi_2)|^2$ and calculate the linearity error \eqref{eq:linerr}. We use classical MDS and MDS with gradient descent (MDS-GD). Figure \ref{fig:linearity_simulation} shows that the system is indeed approximately linear and that the linearity error becomes smaller as the number of reference signals grows. In Figure \ref{fig:linearity_simulation_floor6}, we set the sensitivity threshold to $\tau = 6$ and zero the distances below the threshold per \eqref{eq:phase_problem_with_mask}. 
Again, the linearity error quickly becomes small as the number of anchors increases showing that the overall system is robust and that it allows recovery of phase for small-intensity signals.

Next, we test the linearity error with a real hardware OPU. The OPU gives 8-bit unsigned integer measurements. A major challenge is that the \rev{DMD (see Figure \ref{fig:opu_setup}) only allows binary input signals.} This \rev{is a property of the particular OPU we use and while it} imposes restrictions on reference \rev{design, the method is unchanged as our algorithm \finalrev{does} not assume a particular type of signal}. Section \ref{sec:practical_considerations} in the supplementary material describes how we create binary references and addresses other hardware-related practicalities. 

Figure \ref{fig:linearity_OPU} reports the linearity error on the OPU with suitably designed references \rev{and the same size $\mA$}. The empirically determined sensitivity threshold of the camera is $\tau = 6$, and the measurements below the threshold were not used. We ignore rows of $\mA$ which give points with small norms (less than two) because they are prone to noise and disproportionately influence the relative error. Once again, we observe that the end-to-end system with Algorithm \ref{algo:mpr_algorithm} is approximately linear and that the linearity improves as we increase the number of anchors.

Finally, we demonstrate the magnitude denoising performance. We draw $\va \in \C^{100}$, a random signal $\vxi \in \R^{100}$ and a set of random reference anchor signals. We run our algorithm for number of anchors varying between 2 and 15.
For each number of anchors, we recover $\hat{y}$ for $\abs{y}^2 = \abs{\inprod{\va, \vxi}}^2$ using either classical MDS or MDS-GD. We then measure the number of good bits. The average results over 100 trials are shown in Figure \ref{fig:good_bits}. Figure \ref{fig:good_bits_floor6} reports the same experiment with the sensitivity threshold set to $\tau = 6$ (that is, the entries below $\tau$ are zeroed in the distance matrix per \eqref{eq:phase_problem_with_mask}). 
Both figures show that the proposed algorithm significantly improves the estimated magnitudes in addition to recovering the phases. \rev{The approximately 1 additional good bit with gradient descent in Figure \ref{fig:good_bits_floor6} corresponds to the relative value of ${2^{1}}/{2^8}\approx 0.8\%$ which is consistent with the gradient descent improvement in Figure \ref{fig:linearity_simulation_floor6}.}

\begin{figure*}[t]
    \centering
    \begin{subfigure}[t]{0.329\textwidth}
        \centering
        \includegraphics[width=\textwidth]{figures/linearity_simulation.pdf}
        \caption{}
        \label{fig:linearity_simulation}
    \end{subfigure}
    \begin{subfigure}[t]{0.329\textwidth} 
        \centering
        \includegraphics[width=\textwidth]{figures/linearity_simulation_floor6.pdf}
        \caption{}
        \label{fig:linearity_simulation_floor6}
    \end{subfigure}
   
     
      \begin{subfigure}[t]{0.329\textwidth}
        \centering
        \includegraphics[width=\textwidth]{figures/linearity_OPU_floor6_selected.pdf}
        \caption{}
        \label{fig:linearity_OPU}
    \end{subfigure}
    \caption{\rev{Experiments in simulation and on real hardware to evaluate the linearity error as defined in \eqref{eq:linerr}. The input signals are of dimension $64^2$, $M$ in \eqref{eq:linerr} is 100 and the number of anchors signals are increased. The classical MDS and MDS with gradient descent (MDS-GD) are used. In all cases the error decreases as the number of anchors increases.} (\subref{fig:linearity_simulation}) In simulation \rev{with Gaussian signals and Gaussian reference signals}; (\subref{fig:linearity_simulation_floor6}) In simulation \rev{with Gaussian signals and Gaussian reference signals with sensitivity threshold} $\tau = 6$; 
    (\subref{fig:linearity_OPU}) On \rev{a real} OPU \rev{with binary signals and binary references.}}
\end{figure*}

\begin{figure*}[t]
    \centering
    \begin{subfigure}[t]{0.329\textwidth}
        \centering
        \includegraphics[width=\textwidth]{figures/good_bits_x.pdf}
        \caption{}
        \label{fig:good_bits}
    \end{subfigure}
   
     
    \begin{subfigure}[t]{0.329\textwidth}
        \centering
        \includegraphics[width=\textwidth]{figures/good_bits_x_floor6.pdf}
        \caption{}
        \label{fig:good_bits_floor6}
    \end{subfigure}
   
   
    \begin{subfigure}[t]{0.329\textwidth}
        \centering
        \includegraphics[width=\textwidth]{figures/srls_vs_MDS.pdf}
        \caption{}
        \label{fig:srls_vs_MDS}
    \end{subfigure}
   
     
    \caption{(\subref{fig:good_bits}) Magnitude denoising performance of MDS and MDS-GD \rev{over 100 trials. Input signals are Gaussian and of dimension 100}; (\subref{fig:good_bits_floor6}) Magnitude denoising performance \rev{of MDS and MDS-GD over 100 trials with 100-dimensional Gaussian signals and sensitivity threshold $\tau = 6$}; (\subref{fig:srls_vs_MDS}) \rev{Comparison between recovering a single point and recovering the point and anchors at the same time.} SR-LS is used to locate a single point when anchors are known and MDS is used to locate all points when anchors are unknown.}
\end{figure*}

We also test a scenario where the anchor positions on the complex plane are known exactly \rev{and we only have to localize a single measurement. We compare this to localizing the anchors and the measurements jointly}. \rev{Localizing a single point via multilateration} is performed by minimizing \rev{the} SR-LS objective (see  \eqref{eq:SRLS} in the supplementary material). \rev{The input signal dimension is $64^2$ and we recover $\hat{y}$ for $\abs{y}^2 = |\inprod{\va, \vxi}|^2$. We perform 100 trials} and calculate the SNR of the recovered complex points. Figure \ref{fig:srls_vs_MDS} shows that although having perfect knowledge of anchor locations helps, classical MDS alone does not perform much worse.

\paragraph{Optical randomized singular value decomposition.} 

We use Algorithm \ref{algo:mpr_algorithm} to implement randomized singular value decomposition (RSVD) \rev{as} described in \rev{Halko et al.} \cite{halko2011finding} \rev{on the OPU}. We use 5 anchors in all RSVD experiments. The original \rev{RSVD} algorithm and \rev{a variant with  adaptations for the OPU} are described in Algorithms \ref{algo:rsvd} and \ref{algo:opu_rsvd} in the supplementary material.

\rev{One of the steps in the} RSVD algorithm for an input matrix $\mB \in \R^{M \times N}$ requires the computation of $\mB \mOmega$ where $\mOmega \in \R^{N \times 2K}$ is a standard \rev{real} Gaussian matrix, $K$ is the target number of singular vectors, and $2K$ may be interpreted as the number of \rev{random} projections for each row of $\mB$. \rev{We use the OPU to compute this random matrix multiplication.} An interesting observation is that since \rev{in Algorithm \ref{algo:mpr_algorithm}} we recover \rev{the result of} multiplications by a complex matrix with independent real and imaginary parts, we can halve the number of projections \rev{when using} the OPU with respect to the original algorithm. By treating each row of $\mB$ as an input frame, we \rev{can} obtain $\mY \in \C^{K \times M}$ \rev{via Algorithm \ref{algo:mpr_algorithm}} when $\abs{\mY}^2 = \abs{\mA \mB^\T}^2$ \rev{with} $\mA$ \rev{as defined in Problem \ref{prob:main_problem} with} $K$ rows. Then, \rev{we can construct} $\mP = \begin{bmatrix} \mathrm{Re}(\mY^*) & \mathrm{Im}(\mY^*) \end{bmatrix} \in \R^{M \times 2K}$ \rev{which would be} equivalent to computing $\mB \mOmega$ for real $\mOmega$. Section \ref{sec:rsvd_details} in the supplementary material describes this in more detail.

\rev{Figure \ref{fig:rsvd} shows the results when the OPU is used to perform the random matrix multiplication of the RSVD algorithm on a matrix $\mB$.} Figure \ref{fig:rsvd} (left) reports experiments with \rev{a random binary matrix} $\mB \in \R^{10 \times 10^4}$, different numbers of random projections \rev{(number of rows in $\mA$)}, and ten trials per number of projections. We plot the average error per entry when reconstructing $\mB$ from its RSVD \rev{matrices and singular values}. Next, we take 500 $28 \times 28$ samples from the MNIST dataset \cite{lecun1998gradient}, threshold them to \rev{be} binary, vectorize \rev{them}, and stack \rev{them into} a matrix $\mB \in \R^{500 \times 28^2}$. Figure \ref{fig:rsvd} (right) shows the seven leading right singular vectors reshaped to $28 \times 28$. \rev{The top row shows the singular vectors that are obtained when using the OPU with 500 projections and the bottom row shows the result when using Python.} The error is negligible.
\vspace{-1mm}


\begin{figure*}[ht]
    \centering
    \begin{subfigure}[t]{0.27\textwidth}
        \centering
        \includegraphics[width=\textwidth]{figures/svd_OPU_recon_error.pdf}
       
        \label{fig:svd_OPU_recon_error}
    \end{subfigure}
   
     
    \begin{subfigure}[t]{0.72\textwidth}
        \centering
        \includegraphics[width=\textwidth]{figures/svd_MNIST_OPU_top10_singvecs_half_proj.pdf}
       
        \label{fig:svd_MNIST_OPU_top10_singvecs}
    \end{subfigure}
   
   
    \caption{Left: Average RSVD error \rev{over 10 trials} with varying number of projections on hardware \rev{with an input matrix of size $10 \times 1000$}; Right: \rev{Reshaped leading right singular vectors of an MNIST matrix of size $500 \times 28^2$.} The top rows shows the leading right singular vectors after performing RSVD with the OPU and using our algorithm. The bottom row shows the leading right singular vectors from Python. The relative error is below each singular vector.}
    \label{fig:rsvd}
\end{figure*}

\section{Conclusion}

Traditional computation methods are often too slow for processing tasks which involve large data streams. This motivates alternatives which instead use fast physics to ``compute'' the desired functions. In this work, we looked at using optics and multiple scattering media to obtain linear random projections. A common difficulty with optical systems is that off-the-shelf camera sensors only register the intensity of the scattered light. Our results show that there is nevertheless no need to reach for more complicated and more expensive coherent setups. We showed that measurement phase retrieval can be cast as a problem in distance geometry, and that the unknown phase of random projections can be recovered even without knowing the transmission matrix of the medium. 

Simulations and experiments on real hardware show that the OPU setup combined with our algorithm indeed approximates an end-to-end linear system. What is more, we also improve intensity measurements. The fact that we get full complex measurements allows us to implement a whole new spectrum of randomized algorithms; we demonstrated the potential by the randomized singular value decomposition. These benefits come at the expense of a reduction in data throughput. Future work will have to precisely quantify the smallest achievable data rate reduction due to allocating a part of the duty cycle for reference measurements, though we note that the optical processing data rates are very high to begin with.

\section*{Acknowledgement}
Sidharth Gupta and Ivan Dokmani\'c would like to acknowledge support from the National Science Foundation under Grant CIF-1817577.

\newpage


\input{opu.bbl}

\newpage
\input{supplementary.tex}

\end{document}

\section*{Supplementary material}

