\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}


\section{\uppercase{Introduction}}
\label{sec:introduction}

\noindent
Heart arrhythmias are caused by irregular electrical conduction in cardiac
tissue. Atrial fibrillation, which affects 1--2\% of the population~\cite{Camm2010},
is the most common one. Furthermore, its prevalence increases
with age, from $<$0.5\% at 40--50 years to 5--15\% at 80 years. While not
directly life-threatening, it can lead to serious complications~\cite{January2014}.
In particular, atrial fibrillation is associated with a
3--5 fold increased risk of stroke and a 2-fold increased risk of mortality~\cite{Kannel1998}.
It is also associated with a 3-fold risk of heart failure~\cite{Wang2003}.
Typical symptoms include heart palpitations, shortness of
breath, and fainting. However, about one third of the cases are asymptomatic
which prevents early diagnosis. This, in turn, precludes early therapies which
might protect the patient from the consequences of atrial fibrillation but also
from its progression. Indeed, atrial fibrillation causes electrical and
structural remodeling of the atria which facilitates its further
development~\cite{Frick2001,Nattel2008}.

The gold standard for diagnosing atrial fibrillation and other heart
arrhythmias is the 12-lead ECG. A trained electrophysiologist can select the
most appropriate treatment after reviewing ECG signals and the patient history.
This is, however, a time-consuming task, especially for long recordings such as
the ones collected with Holter monitors. To alleviate this task, several
approaches have been proposed to detect arrhythmias from ECG
signals~\cite{Owis2002,DeChazal2004}. Even without perfect detection accuracy, these
approaches are useful as they facilitate reviewing ECG by selecting relevant
signal excerpts.

Recently, neural networks have shown impressive performance in various
classification and regression tasks. Image processing was the first field where
deep networks surpassed existing approaches by a large margin~\cite{Krizhevsky2012}.
Since then, they have been extensively applied to fields
previously dominated by signal processing. In particular, several architectures
have been proposed to detect and classify heart arrhythmias from ECG signals.

In the context of the challenge of Computing in Cardiology 2017~\cite{Clifford2017},
a few neural network architectures were proposed to
classify single lead ECG signals into one of the following classes: normal
sinus rhythm, atrial fibrillation, other rhythm, and noise. One of these
architectures uses logarithmic spectrograms computed over sliding windows of
ECG as input to two-dimensional convolutional layers~\cite{Zihlmann2017}.
Aggregation of successive windows is done with either temporal averaging or a
recurrent layer. However, the convolutional and recurrent layers were trained
separately due to convergence issues. A similar approach used a 16-layer
convolutional neural network with skip connections to classify arrhythmias from
ECG signals~\cite{Xiong2017}. Each layer is composed of batch normalization,
ReLU activation, dropout, one-dimensional convolution, and averaging pooling.

Recently, a convolutional neural network was shown to reach cardiologist-level
arrhythmia detection~\cite{Rajpurkar2017}. This 34-layer network takes
advantage of a very large dataset of 64,121 ECG signals, recorded from 29,163
patients, to recognize 12 different heart arrhythmias including atrial
fibrillation, atrial flutter, and ventricular tachycardia. Another approach
applied convolutional neural networks for detecting atrial fibrillation from
time-frequency representations of ECG signals~\cite{Xia2018}. Two methods for
computing these representations were compared: the short-time Fourier transform
or the stationary wavelet transform. In this case, the neural network using
coefficients from the second transform yielded better classification accuracy.

Neural networks have thus shown promising results for the detection of abnormal
cardiac rhythms. Furthermore, as mentioned previously, it is of the utmost
importance to detect arrhythmias as early as possible to improve treatment
outcome. To tackle this issue, we developed a system composed of two main
elements: a smart vest to record ECG signals and an algorithm to detect and
classify arrhythmias. The smart vest includes two cooperative sensors to record
a single lead ECG signal and stream by Bluetooth the collected data to a
gateway. This gateway then forwards the ECG signal to the cloud where it is
processed by a neural network in order to detect abnormal rhythms.

This article is structured as follows. First, the dataset of ECG signals, our
neural network architecture, and the monitoring system are described in
Section~\ref{sec:methods}. Then, the results are presented in
Section~\ref{sec:results} and discussed in Section~\ref{sec:discussion}.
Finally, this article ends with a short conclusion in
Section~\ref{sec:conclusion}.

\section{\uppercase{Materials and Methods}}
\label{sec:methods}

\noindent

\subsection{Dataset}

We used the dataset from the challenge of Computing in Cardiology 2017 to train
a neural network to classify cardiac arrhythmias. This dataset includes 8528
single lead ECG signals recorded with an AliveCor device. The signals are
sampled at 300~Hz and have durations ranging from 9 to 60 seconds. Each signal
was acquired when the subject held each one of the two electrodes in each hand
resulting in a lead I (left arm -- right arm) ECG. As the device has no
specific orientation, many signals are inverted (right arm -- left arm).

All ECG signals are labeled with one of the following four classes: normal
sinus rhythm, atrial fibrillation, other rhythm, and noise. The proportion of
each class in the dataset varies from 3.27\% for noise to 59.52\% for normal
rhythm. The full breakdown of all classes is reported in
Table~\ref{table:breakdown}.

\begin{table}
    \centering
    \scriptsize
    \caption{Number of signals and mean duration for each class.}
    \renewcommand{\arraystretch}{1.3}
    \begin{tabular}{@{}lrrr@{}}
    \toprule
    Class & Count & Proportion & Mean duration [s] \\
    \midrule
    Normal rhythm & 5076 & 59.52\% & 32.11 \\
    Atrial fibrillation & 758 & 8.89\% & 32.34 \\
    Other rhythm & 2415 & 28.32\% & 34.30 \\
    Noise & 279 & 3.27\% & 24.38 \\
    Total & 8528 & 100\% & 32.50 \\
    \bottomrule
    \end{tabular}
    \label{table:breakdown}
\end{table}

The entries of the challenge of Computing in Cardiology 2017 were ranked
according to the following score evaluated on a private test set:
\begin{equation}
    S_{\text{CinC}} = \frac{F_{1n} + F_{1a} + F_{1o}}{3} \label{eq:cinc2017}
\end{equation}
where $F_{1n}$, $F_{1a}$, and $F_{1o}$ denote the $F_1$ scores for normal
rhythm, atrial fibrillation, and other rhythm. The four
winners~\cite{Teijeiro2017,Datta2017,Zabihi2017,Hong2017} reached a score of 0.83.

Several aspects of this dataset are challenging for arrhythmia classification.
First, many signals are inverted as mentioned previously. Second, the classes
are not balanced. There are very few signals labeled atrial fibrillation or
noise compared to the ones labeled normal rhythm. Furthermore, the durations of
the recordings are also different. They vary from 9 to 60 seconds.
Figure~\ref{fig:durations} illustrates this issue. Most ECG signals last 30
seconds but a significant portion has shorter or longer durations. In addition,
labeling is relatively coarse as a single label is associated with each ECG
signal. In some cases, several labels could be used for the same signal.
Finally, the ECG quality of a non-negligible part of the records is rather
poor. Four examples of signals are shown in Figure~\ref{fig:examples}. The
first two signals are labeled normal rhythm and atrial fibrillation and have
good overall quality. The third example is a normal rhythm record with
acceptable quality except for a short segment of noise. The last ECG signal is
labeled as atrial fibrillation but has very poor quality due to large shifts of
the baseline. It also illustrates that all signals do not share the same
duration.

\begin{figure}
    \centering
    \includegraphics[width=75mm]{durations-eps-converted-to.pdf}
    \caption{Record durations sorted in increasing order.}
    \label{fig:durations}
\end{figure}

\begin{figure}
    \centering
    \includegraphics[width=75mm]{examples-eps-converted-to.pdf}
    \caption{Examples of ECG records: (A) normal rhythm from record A00001, (B)
        atrial fibrillation from record A00004, (C) normal rhythm from record
        A00002, (D) atrial fibrillation from record A00015.}
    \label{fig:examples}
\end{figure}

The test set used during the challenge of Computing in Cardiology was not yet
publicly released. Therefore, we split the dataset (which was originally
intended only for training) into a training set of 7000 records and a test set
with the remaining 1528 records while approximately preserving class
repartition.

We applied the following pre-processing steps to the dataset. First, we
filtered the ECG signals with a digital Butterworth band-pass filter between
0.5 and 40~Hz. The filter is applied in both forward and backward directions to
avoid distortions. The analog filter used by the recording device has similar
cutoff frequencies but still leaked some components outside the pass-band.
Then, we downsampled the signals to 200~Hz to reduce the number of samples.
Finally, the signals are scaled by the mean of the standard deviations of all
signals from the training set. Scaling is helpful to accelerate training~\cite{LeCun2012}.

\subsection{Network Architecture}

An approach to handle signals with different durations is to extract windows
with the same length. The label of a signal is then used for all included
windows. Fixed size inputs would make possible to apply a convolutional neural
network to learn features useful to classify arrhythmias. However, this
approach is sub-optimal as illustrated in Figure~\ref{fig:examples} where a
signal labeled as normal rhythm includes a segment of noise. Furthermore, all
signals would need to be truncated to the same length. This would lead to a
large data loss due to the considerable variations in signal durations.

A more appropriate approach is to use a recurrent neural network which is
well-suited to process sequences with varying lengths. While such a neural
network can, by design, remember past values over long time intervals, they are
not as efficient as convolutional neural networks for learning complex
features.

After reviewing the advantages and drawbacks of these two approaches, we chose
to combine them and build a neural network that includes convolutional and
recurrent layers. Specifically, each ECG signal was divided in sliding windows
of 512 samples with 50\% overlap. This corresponds to a window duration
slightly above 2.5 seconds. The number of windows extracted from each signal
depends on its duration. Seven convolutional layers were applied to all windows
of a signal. Each convolutional layer is composed of a 1D convolution and a
max pooling operation. The convolution uses a kernel of size 5, zero padding,
and a ReLU activation~\cite{Hahnloser2000}. The pool size for max pooling was
set to 2. The first convolutional layer has 8 output channels (from the single
channel windows). Then, each following layers double the number of channels
while max pooling halves the number of samples. After the convolutional layers,
a global averaging pooling layer was applied. This results in 512 features for
each input window. The features were then processed with a long short-term
memory (LSTM) layer~\cite{Hochreiter1997} with 128 units. Finally, a softmax
layer outputs the probability of each class for the input ECG windows.
This results in a neural network with 1,203,364 trainable parameters including
874,656 for the convolutional part, 328,192 for the recurrent part, and 516 for
the final softmax layer. The complete network architecture is summarized in
Table~\ref{table:architecture}.

\begin{table}
    \centering
    \caption{Neural network architecture. The output size is given as
        $N_w \times N_s \times N_c$ where $N_w$ denotes the variable number of
        windows, $N_s$ the number of samples, and $N_c$ the number of
        channels.}
    \renewcommand{\arraystretch}{1.3}
    \begin{tabular}{@{}lc@{}}
        \toprule
        Layer & Output size \\
        \midrule
        Input windows & $N_w \times 512 \times 1$ \\
        Convolutional layer 1 & $N_w \times 256 \times 8$ \\
        Convolutional layer 2 & $N_w \times 128 \times 16$ \\
        Convolutional layer 3 & $N_w \times 64 \times 32$ \\
        Convolutional layer 4 & $N_w \times 32 \times 64$ \\
        Convolutional layer 5 & $N_w \times 16 \times 128$ \\
        Convolutional layer 6 & $N_w \times 8 \times 256$ \\
        Convolutional layer 7 & $N_w \times 4 \times 512$ \\
        Global average pooling & $N_w \times 512$ \\
        LSTM layer & $128$ \\
        Softmax layer & $4$ \\
        \bottomrule
    \end{tabular}
    \label{table:architecture}
\end{table}

\subsection{Data Augmentation}

As the dataset is relatively small for fitting a neural network, we applied
different strategies to synthetically augment the number of ECG signals
available during training. The first strategy is to simply flip the sign of
each signal with probability 0.5. We found it easier to let the neural network
learn to take into account inverted ECG signals instead of applying a
rectifying step during pre-processing.

Furthermore, when extracting sliding windows, it is not possible to use all
samples for the large majority of ECG signals. Indeed, the maximum number of
sliding windows $N_w$ in a signal with $N$ samples is given by
\begin{displaymath}
    N_w = \left\lfloor \frac{N - 512}{256} \right\rfloor + 1
\end{displaymath}
assuming $N \ge 512$. In the previous expression, $\lfloor\cdot\rfloor$ denotes
the floor function. We took advantage of this fact to place the first window
at a random offset from the start of the signal. This random offset is drawn
uniformly from
\begin{displaymath}
    \{0, 1, \dots, N - (N_w - 1) \cdot 256 - 512\}
\end{displaymath}
for each signal at each epoch. The main idea behind this strategy is to prevent
the neural network from learning the exact positions of the QRS complexes in
the training set. However, we always used the maximum number of sliding windows
possible for each signal to avoid wasting ECG samples.

The third strategy we applied is to resample each signal at each epoch with
probability 0.8 in order to simulate slightly slower or faster heart rate and
thus help the neural network to reach better generalization performance.
Naturally, the resampling operation should not change the heart rate too much.
Otherwise, there is a risk to confuse a cardiac rhythm for another one.
Therefore, if a signal needs resampling, its length is changed by a proportion
sampled uniformly between $-5$\% and $+5$\%.

\subsection{Training}

We implemented our neural network and strategies for data augmentation in
Python with the Keras library~\cite{Chollet2015}. We trained the neural network
for 200 epochs by minimizing the cross-entropy with the Adam
algorithm~\cite{Kingma2014}. We set the initial learning rate to 0.001. The learning rate
was divided by two if the cross-entropy evaluated on the test set did not
decrease for 5 consecutive epochs with a lower limit at $10^{-5}$.

We used a batch size of 50 signals. As a batch must include the same number of
sliding windows for each signal, we applied zero-padding. Specifically, too
short signals were prepended with all-zero windows. To limit zero-padding as
much as possible we sorted the signals by duration and grouped them in batches
of similar lengths. This resulted in batches with varying numbers of windows.

The LSTM layer was regularized  by applying dropout with a rate of 0.5 for both
the input and recurrent parts~\cite{Srivastava2014,Gal2016}. We monitored the
accuracy on the test set and selected the weights at the best epoch as the
final parameters of the neural network.

\subsection{Monitoring System}

The smart vest used to monitor ECG includes two cooperative sensors
illustrated in Figure~\ref{fig:sense_sensors}. These sensors use dry stainless
steel bi-electrodes. This technology, which was validated in a previous study,
yields high quality measurements of ECG and bio-impedance signals, even in
motion. No wetting of the electrode is required. Moreover, the electrical
connection linking both sensors does not have to be shielded, nor insulated,
which makes its integration in garment easier and cheaper (a conductive fabric
is sufficient). The length of the connection is not limited to very short
distances and is therefore placed in the back so as to have space for a central
zipper in the vest, which makes donning and doffing easy, even for the elderly.
The two watertight sensors are clipped in the vest and can be removed for
washing and recharging. Operation is simplified to its maximum: the sensors
automatically switch to record mode as soon as they are applied on the skin and
return to standby mode when removed.

\begin{figure}
    \centering
    \ifthenelse{\boolean{blind}}{
        \includegraphics[width=75mm]{pixelized_sense_sensors.jpg}
    }{
        \includegraphics[width=75mm]{sense_sensors.jpg}
    }
    \caption{Cooperative sensors included in the smart vest for monitoring
        ECG.}
    \label{fig:sense_sensors}
\end{figure}

This smart vest which was originally developed for athletes also monitors the
following biomedical signals: heart rate, transthoracic impedance, respiration
rate, skin temperature, activity class (resting, walking, running), and posture
(lying, standing/sitting). A quality index is associated with each signal so
that the reliability of the measurements can be easily assessed. However, these
additional signals and the quality indices were not used in the present study.
Figure~\ref{fig:sense_vest} shows the smart vest during the validation
protocol.

\begin{figure}
    \centering
    \ifthenelse{\boolean{blind}}{
        \includegraphics[width=75mm]{pixelized_sense_vest.jpg}
    }{
        \includegraphics[width=75mm]{sense_vest.jpg}
    }
    \caption{Smart vest during the validation protocol.}
    \label{fig:sense_vest}
\end{figure}

Another important feature of the smart vest consists in its capability to store
all the recorded signals locally as well as synchronize them in the cloud. This
is achieved by means of an accompanying gateway that acts as a relay for the
streamed data. The gateway is robust to poor network connectivity and uploads
the newly available data to the cloud as it is being streamed while an active
connection is available. In our case, the gateway is composed of a simple
Raspberry Pi that connects to the smart vest via Bluetooth. The biomedical data
is then collected, stored locally in a compact format and relayed via telemetry
messages simultaneously to all the clients connected to this gateway as well as
a private cloud. The messaging protocol chosen for this application is
MQTT\footnote{https://mqtt.org} (Message Queuing Telemetry Transport). This
protocol, which uses the publish/subscribe paradigm, is one of the most widely
employed telemetry protocols in IoT and real-time streaming. Once the ECG
signal is in the cloud, we can leverage the powerful computing capabilities and
detect heart arrhythmias with our neural network.

The ECG signals recorded with the smart vest device were pre-processed
similarly to the data from the challenge of Computing in Cardiology 2017.
First, the same band-pass filter between 0.5 and 40~Hz was used to remove the
baseline and high-frequency noise. Then, the signals were resampled at 200~Hz.
Finally, we scaled each signal with its standard deviation. Indeed, we could
not use the scaling factor computed on the training set as the two types of
ECG signals did not have the same range of values.

After pre-processing the ECG signals, we extracted sliding windows of 512
samples with 50\% overlap. We used groups of 25 such windows as input to the
neural network. We selected this specific number of windows as it corresponds
to segments of approximately 33 seconds which is close to the median length of
the signals used for training and testing the model. Applying the neural
network resulted in a rhythm prediction for each group of 25 windows.

\section{\uppercase{Results}}
\label{sec:results}

\noindent
We evaluated three configurations. The first one used the architecture
described in Table~\ref{table:architecture} except that it included only 6
convolutional layers. In addition, all data augmentation strategies were
applied during training except resampling. The second configuration used the
full architecture with 7 convolutional layers but again without resampling.
Finally, the last configuration was identical to the second one but with
resampling for further data augmentation. Figure~\ref{fig:training} shows the
evolution of the cross-entropy loss and the accuracy during training for these
three network configurations. An additional convolutional layer helped to
increase the accuracy and reduce the loss. Resampling the signals also slightly
improved performance. Despite our efforts, we could not completely eliminate
over-fitting as shown by the performance gap between training and test sets.
Indeed, additional regularization only decreased the performance of the
network. However, it is worth mentioning that resampling helped to reduce
over-fitting. The neural network with the third configuration reached an
accuracy of 87.37\% on the test set after 47 epochs. At the same epoch, the
accuracy on the training set was 90.90\%.

\begin{figure}
    \centering
    \includegraphics[width=75mm]{training-eps-converted-to.pdf}
    \caption{Cross-entropy loss (top) and accuracy (bottom) evaluated during
        training for a network with 6 convolutional layers without resampling
        (blue), a network with 7 convolutional layers without resampling
        (orange), and a network with 7 convolutional layers with resampling
        (green). Dashed lines denote results obtained on the training set and
        solid lines results obtained on the test set.}
    \label{fig:training}
\end{figure}

After selecting the best neural network, we evaluated it without zero-padding
by selecting a batch size of 1. In this case, the accuracy was 92.64\% on the
training set and 87.50\% on the test set. Sensitivity, specificity, and
$F_1$ score values are reported in Table~\ref{table:metrics} for all classes.
Unsurprisingly, the best $F_1$ score was obtained for the class with the most
samples (normal rhythm) and the worst one for the class with the least samples
(noise). Interestingly, the specificity for atrial fibrillation was relatively
high at 0.9784 while the sensitivity was lower at 0.8382. Thus, the number of
false positive is more larger than the number of false negative which is an
important property of the model. Indeed, it missed only a few atrial
fibrillation cases while false detections can always be disproved with
additional analyses such as a full 12-lead ECG. Furthermore, when evaluated in
terms of the score used during the challenge of Computing in Cardiology 2017
\eqref{eq:cinc2017}, the neural network yielded 0.9156 and 0.8495 on the
training and test sets. This is similar to the best score obtained by the
winners of the challenge. However, we could only evaluate our neural network on
1528 ECG signals from the original training set since the official test set was
not made public yet. Therefore, it is difficult to compare the performance of
our approach with the winning entries.

\begin{table}
    \centering
    \small
    \caption{Performance metrics for the best neural network.}
    \renewcommand{\arraystretch}{1.3}
    \begin{tabular}{@{}llcc@{}}
        \toprule
        Class & Metric & Training set & Test set \\
        \midrule
        \multirow{3}{*}{Normal rhythm} & Sensitivity & 0.9707 & 0.9263 \\
        & Specificity & 0.8955 & 0.8853 \\
        & $F_1$ score & 0.9509 & 0.9243 \\
        \addlinespace
        \multirow{3}{*}{Atrial fibrillation} & Sensitivity & 0.9084 & 0.8382 \\
        & Specificity & 0.9942 & 0.9784 \\
        & $F_1$ score & 0.9232 & 0.8143 \\
        \addlinespace
        \multirow{3}{*}{Other rhythm} & Sensitivity & 0.8380 & 0.7921 \\
        & Specificity & 0.9675 & 0.9352 \\
        & $F_1$ score & 0.8728 & 0.8099\\
        \addlinespace
        \multirow{3}{*}{Noise} & Sensitivity & 0.9345 & 0.7600 \\
        & Specificity & 0.9972 & 0.9871 \\
        & $F_1$ score & 0.9264 & 0.7103 \\
        \bottomrule
    \end{tabular}
    \label{table:metrics}
\end{table}

After training and evaluation on the dataset of the challenge of Computing in
Cardiology 2017, we applied the neural network to a few ECG signals recorded
with the smart vest. We also extracted the times between consecutive R-waves
from the ECG signals (a.k.a. RR intervals) to facilitate visualization of the
heart rate. An example is shown in Figure~\ref{fig:sense_normal_other_normal}.
In this case, the subject had a normal sinus rhythm with two successive
premature ventricular contractions characterized a very short RR interval
followed by a long recovery RR interval. These premature contractions are
clearly identifiable especially compared to the stable RR intervals of a normal
rhythm. The neural network classified the this ECG segment as normal rhythm
except for the groups of windows including the premature contractions which are
classified as other rhythm.

\begin{figure*}
    \centering
    \includegraphics[width=155mm]{normal_other_normal-eps-converted-to.pdf}
    \caption{Cardiac rhythm classification from an ECG signal recorded with the
        smart vest. The ECG signal (top) includes two premature ventricular
        contractions which are visible in the RR intervals (middle). The neural
        network correctly identifies this segment of the signal as reflected by
        the class probabilities (bottom).}
    \label{fig:sense_normal_other_normal}
\end{figure*}

Another example of classification is shown in Figure~\ref{fig:sense_noise}. In
this case, no ECG signal was collected during a few minutes due to poor
skin-electrode contact caused by motion. The neural network identified this
segment without signal as noise. Furthermore, two segments preceding signal
loss were classified as normal rhythm and other rhythm. The RR intervals were
stable in the first segment and there were a few premature ventricular
contractions in the second one. After the ECG signal was recovered, the RR
intervals varied widely and the neural network classified the data as other
rhythm.

\begin{figure*}
    \centering
    \includegraphics[width=155mm]{noise-eps-converted-to.pdf}
    \caption{Cardiac rhythm classification from an ECG signal recorded with the
        smart vest. The ECG signal (top) includes a segment with no R waves due
        to poor skin-electrode contact. Consequently, the RR intervals (middle)
        could not be extracted. The class probabilities (bottom) computed by
        the neural network are valid as the segment is labeled as noise.}
    \label{fig:sense_noise}
\end{figure*}

\section{\uppercase{Discussion}}
\label{sec:discussion}

\noindent
The neural network architecture we developed achieved a classification
performance similar to the winners of the challenge of Computing in Cardiology
2017. Although we could not evaluate our approach on the same data since the
true test set was not released yet, the results are promising. In particular,
the specificity for atrial fibrillation shows that only a limited number cases
are not detected. We could take into account signals of different lengths by
combining convolutional and recurrent layers. Indeed, the convolutional layers
extract features relevant for arrhythmia classification from sliding windows of
raw ECG signal. With this approach, there is no need to pre-process the data to
extract spectrogram or wavelet transform for instance. Furthermore, no feature
engineering is required as the neural network learns during the training phase
to extract high-level features useful for classification. The only
pre-processing we used is band-pass filtering as well as scaling to accelerate
training. We also used strategies for data augmentation to reduce over-fitting
and improve the generalization performance of the network. In particular,
randomly flipping the sign of ECG signals forced the neural network to learn to
take into account both regular and inverted waveforms. Without this simple
strategy, we would have needed to develop a robust method to detect signals
that had to be rectified.

In addition, we have shown that it is possible to apply a neural network
trained on a generic dataset to ECG signals recorded by the smart vest with
minimal changes. Indeed, there was no need to adapt the network architecture.
The method to standardize the ECG signals was the only element that required a
modification. Of course, if the two measurement systems, the AliveCor device
and the smart vest, did not record the same ECG lead, additional changes would
be needed. However, it is difficult to determine the extend of these changes
without evaluating the neural network on a dataset of ECG signals from another
lead.

Taken together, these results demonstrate that our system could be used to
monitor the cardiac activity of a subject for detecting abnormal rhythms over
long periods of time. Indeed, the smart vest is more comfortable to wear than
a traditional Holter and it can still collect high quality ECG signals since it
was initially developed for athletes. After transmitting the data to the cloud,
our neural network can quickly classify abnormal rhythms. While the accuracy
of our system is not perfect, it can still help to reduce the time spent
reviewing ECG by selecting segments with potential abnormalities that require
additional attention. These segments can then be analyzed by a trained
specialist. If needed, a full 12-lead ECG can be performed to refine or confirm
the diagnosis.

\section{\uppercase{Conclusion}}
\label{sec:conclusion}

\noindent
We presented a system composed of a smart vest to record a single lead ECG
signal and a neural network for detection and classification of heart
arrhythmias. We plan to aggregate several databases of ECG signals with rhythm
annotations to extend the types of arrhythmias that our algorithm can detect
and improve its accuracy. We will also investigate whether adding skip
connections in our network architecture further improves classification
performance.

\section*{\uppercase{Acknowledgements}}

\noindent
We would like to thank the anonymous reviewers for their valuable suggestions
and comments.

\bibliographystyle{apalike}
{\small