\section
In the spontaneous parametric downconversion (SPDC) process, portion of pump photons split into photon pairs created with energies strictly satisfying the energy conservation condition and scattered in directions promoted by the momentum conservation (or phase-matching) condition.
The SPDC photons are thus entangled by nature in energy-time \cite{hong1987measurement,franson1989bell} as well as in momentum-position \cite{mair2001entanglement,lima2006propagation,barbieri2007complete, walborn2010spatial,leach2010quantum}. However for the polarization, the remaining degree of freedom (DoF), its entanglement appears always as a consequence of subtle designs using either a single \cite{kwiat1995new,Hugues2006non,kim2006phase,ueno2012entangled} or more nonlinear crystals \cite{kwiat1999ultrabright,kim2001interferometric,ljunggren2006theory,hegazy2017orthogonal,Villar2018Experimental}.

In either case, the wavefunctions of the two polarization possibilities (HH and VV for $|\phi^\pm\rangle$, or HV and VH for $|\psi^\pm\rangle$) can be made almost identical in magnitude, while a comprehensive picture about the indistinguishability in the space and time domains is provided by the directional-spectral relative-phase function. Generally, a flatter relative-phase function in some domain implies higher indistinguishability in its conjugate Fourier-transform domain
 \cite{hegazy2017orthogonal,Hegazy2017CLEO}.
Therefore, the spatial and temporal compensations actually correct (or flatten) the varying relative-phase function in direction and frequency. Several works previously addressed spatial and temporal compensations at specific spatial modes \cite{kwiat1995new,nambu2002generation,ljunggren2006theory,trojek2008collinear} or over relatively narrow aperture \cite{altepeter2005phase,rangarajan2009optimizing,cialdi2010programmable} by the use of additional birefringent elements.
Over such limited spatial extents, well-designed birefringent elements can effectively clean for the linear dependence of the relative phase function and restore the indistinguishability \cite{hegazy2015tunable}.
This is not the case when photon pairs are collected over relatively wide angular extents, where higher-order contributions of the relative-phase function drastically
reduce the effectiveness of linear methods.

In this paper, we present and demonstrate experimentally a method for eliminating the spatial and temporal distinguishability of a polarization-entangled state over wide emission angles.
We initially derive the two-photon state emitted by a noncollinear SPDC; including the directional-spectral relative-phase function all over the SPDC cone.
It is noticed that while the time-delay map (the spatial distribution of spectral relative-phase gradient) has a negligible change over wide angular extent, the directional relative-phase exhibits obvious quadratic profile. To restore the temporal and spatial indistinguishability, the compensation of the former can be effectively done in the conventional way; using a birefringent element. The latter is manipulated using a two-dimensional spatial light modulator (2D SLM) loaded with the inverted relative-phase profile. A fast yet accurate experimental technique is introduced to directly measure the relative-phase profile; dispensing with the quantum state tomography.

Consider a diagonally polarized pump beam illuminating a couple of crossed type-I crystals which are identically cut to promote production of degenerate photons in a noncollinear SPDC geometry.
The pump beam is classically treated as a superposition of frequency and plane-wave  modes:
$E_{{p}}(\mathbf{x},z,t)= \int d\omega _{{p}}d\mathbf{q}_{{p}}~A_{p}(\omega_{{p}};\mathbf{q} _{{p}})\exp i\left(\kappa _{{p}}z+\mathbf{q}_{{p}}. \mathbf{x}-\omega _{{p}}t\right) +\mathrm{c.c.}$,
with angular frequency $\omega _{p}$ and wavevector $\mathbf{k}_{p}=(%
\mathbf{q}_{{p}},\kappa _{{p}})$ where $\mathbf{q}_{{p}%
} $ are the transverse components along the coordinates $\mathbf{x}=(x,y)$, and $\kappa _{p}$ is the longitudinal component.
The pump polarization components in horizontal/vertical (H/V) basis have a relative phase $\phi_p$ that can be manipulated by tilting a birefringent element; e.g., a quarter wave plate (QWP).

Similarly, the produced SPDC two-photon emission can be expressed –based on the SPDC creation operators– as a joint spectral and spatial expansion in monochromatic planar waves with the angular frequencies $\omega_{1,2}$ and the transverse wave vectors  $\mathbf{q}_{\mathrm{1,2}}$ (subscripts 1, 2 denote signal and idler photons, respectively).
The longitudinal wavevectors of the interacting waves can be expressed as dictated by Maxwell's equations for ordinary (\textit{o}) and extraordinary (\textit{e}) polarization as
\cite{born2013principles,walborn2010spatial,hegazy2017orthogonal}
\begin{equation}
\begin{split}
&\kappa _{{j}}^{o }(\omega_{{j}};\mathbf{q}_{{j}})=\sqrt{( {\omega_{{j}} ~n_{{j}}^{o}}/{c}) ^{2}-\vert \mathbf{q}_{{j}}\vert
^{2}}\\
&\kappa _j ^{e(H,V)}(\omega_j;\mathbf{q}_j)
=   q_{j(x,y)} \tan\rho^\perp_j \\
&~~~~~~~~~~~~+ n_j^{e^\perp} \left[ \left(\tfrac{\omega _j}{c} \right) ^{2} -         \tfrac{1}{(n_j^e)^2}   q_{j(y,x)}^2-       (\tfrac{{n_j^{e^\perp}}}{n_j^e n_j^o} )^2         q_{j(x,y)}^2 \right]^{\frac{1}{2}},\\
\end{split}
\label{eq:Longwavevector}
\end{equation}
where the subscript $j=p, 1, 2$, the superscripts $(H)$ and $(V)$ label the SPDC crystals with optic axes lying in the $xz$ (horizontal) and $yz$ (vertical) plane respectively, $c$ is the speed of light in space, $n_{{j}}^{o}$ and $n_{{j}}^{e}$ are principal values of refractive index, and $n_j^{e^\perp}$ and $\rho^\perp_j$ are  the refractive indices and walk-off angles of the extraordinary-polarized ray propagating along $z$ axis (perpendicular to the crystals interfaces).
The SPDC in the two crystals creates a two-photon state entangled in every degree of freedom; frequency, momentum, and polarization, which is expressed at small emission angles  \cite{migdall1997polarization} by the superposition
\begin{equation}
\begin{split}
\left\vert \psi \right\rangle   \sim  \mathbf{\int}d\omega _{1}d\omega _{2}d\mathbf{q}_{1}d%
\mathbf{q}_{2} &{\large \{}\Phi _\mathrm{HH}\left(\omega _{1},\omega _{2};\mathbf{q}%
_{1},\mathbf{q}_{2}\right)\left\vert H_1H_2
\right\rangle \\
&+~\Phi _\mathrm{VV}\left( \omega _{1},\omega
_{2};\mathbf{q}_{1},\mathbf{q}_{2}\right) \left\vert V_1V_2
\right\rangle{\large \}},
\end{split}
\label{ent1}
\end{equation}
where $\Phi _\mathrm{HH,VV}\left(\omega _{1},\omega _{2};\mathbf{q}_{1},\mathbf{q}_{2}\right)$ are the biphoton wavefunctions corresponding to the HH and VV possibilities.
 Assuming that the crystals are of infinite transverse extent and pumped by a non-depleted beam with the reflected waves at all crystal interfaces being negligible, the two-photon wavefunctions can be written as
\begin{equation*}
\begin{split}
\Phi_\mathrm{HH} = ~&\tfrac{\chi^{(2)}}{\sqrt[]{2}} A_p(\omega_1+\omega_2;\mathbf{q}_1+\mathbf{q}_2) \int_{-2L}^{-L} dz~ e^{i \int_0^z dz^{\prime } \Delta \kappa(z^{\prime })}, \\
\propto ~&
e^{-i L\left(\tfrac{1}{2}{\Delta \kappa_{V}^{ooe}} + \Delta \kappa_H^{eeo}  \right)}
L
A_p(\omega_1+\omega_2;\mathbf{q}_1+\mathbf{q}_2) ~ \mathrm{sinc} \left( \tfrac{ \Delta \kappa_{V}^{ooe} L}{2}\right), \\
\Phi_\mathrm{VV} = ~&  e^{-i\phi_p}~ \tfrac{\chi^{(2)}}{\sqrt[]{2}}
A_p(\omega_1+\omega_2;\mathbf{q}_1+\mathbf{q}_2) \int_{-L}^{0} dz~ e^{i \int_0^z dz^{\prime } \Delta \kappa(z^{\prime })}, \\
\propto ~&
 e^{-i\left(\phi_p+\tfrac{1}{2}{\Delta \kappa_{H}^{ooe}L} \right)}
L
A_p(\omega_1+\omega_2;\mathbf{q}_1+\mathbf{q}_2)  ~\mathrm{sinc} \left(\tfrac{\Delta \kappa_{H}^{ooe} L}{2}\right), \\
\end{split}
\label{Eq: amplitudes}
\end{equation*}
where
$L$ is the crystal length, $\chi ^{(2)}$ is the bilinear susceptibility,
$~\Delta \kappa^{ooe}_{H,V} =\kappa _{p}^{e(H,V)}-\kappa _{1}^{o}-\kappa _{2}^{o}$ are  wavevector mismatches  within the interacting $(H)$ and $(V)$ crystals,
$\Delta \kappa_H^{eeo}=\kappa _{p}^{o}-\kappa _{1}^{e(H)}-\kappa_{2}^{e(H)}$ is the mismatch of waves interacting in the first crystal when passing through the second one ($\Delta \kappa_H^{eeo} \gg \Delta \kappa^{ooe}_{H,V}$, therefore no downconversion is considered in the later case).

The relative phase of the produced state $|\psi\rangle$ is then
\begin{equation}
\begin{split}
\vartheta
=\mathsf{arg}\left\{  \Phi_\mathrm{VV} / \Phi_\mathrm{HH} \right\}=
\left[\tfrac{1}{2}({\Delta \kappa_{V}^{ooe} }-\Delta \kappa_{H}^{ooe}) + \Delta \kappa_H^{eeo}  \right] L -\phi_p,
\end{split}
\label{eq:relative_phase}
\end{equation}
which includes phase-matching terms and an initial phase term. We assume that the spatiotemporal amplitude of the pump beam is factorizable:
$A_p(\omega_p;\mathbf{q}_p)= ({\sqrt[]{\tau_p}}/{\sqrt[4]{\pi}})\,\exp [-\tfrac{1}{2} \tau_p^2 (\omega_{{p}}-\omega_p^0)^2 ] \,\,\, \delta(\mathbf{q}_p),$ which expresses a polychromatic plane wave with a coherence time $\tau_p$. Therefore, the biphoton spatial characteristics can be fully determined with reference to the signal photon alone ($\mathbf{q}_2 = -\mathbf{q}_1$). This assumption implies also that $\Delta \kappa_{V}^{ooe} \approx\Delta \kappa_{H}^{ooe}\equiv\Delta \kappa^{ooe}$ as dictated by Eq. (\ref{eq:Longwavevector}); and consequently $|\Phi_\mathrm{HH}|\approx|\Phi_\mathrm{VV}|$. By substituting Eq. (\ref{eq:Longwavevector}) into Eq. (\ref{eq:relative_phase}) and ignoring the insignificant higher-order terms, the relative phase writes
\begin{equation}
\begin{split}
    \vartheta  (\omega _{1},\omega _{2};  \mathbf{q}_{1},  & -\mathbf{q}_{1}) \approx\\
\Big \{ & \frac{(n_p^o - n_1^{e^\perp}) \omega_1 + (n_p^o -  n_2^{e^\perp}) \omega_2}{c}\\
&- q_{1,x} (\tan \rho_1^\perp   -   \tan \rho_2^\perp ) \\
&+  \frac{c q_{1,x}^2}{2}    \Big[ \frac{(n_1^{e^\perp})^3}{\omega_1  (n_1^e )^2 (n_1^o)^2}  + \frac{(n_2^{e^\perp} )^3}{\omega_2  (n_2^e )^2 (n_2^o )^2}\Big] \\
&+ \frac{c q_{1,y}^2}{2} \Big[ \frac{n_1^{e^\perp}}{\omega_1 (n_1^e )^2} + \frac{n_2^{e^\perp}}{\omega_2 (n_2^e )^2 }\Big] \Big\} L - \phi_p,
\end{split}
\label{eq:relative_phase_approx}
\end{equation}
which is a quadratic function in the transverse wavevector components of the SPDC photons.
 Moreover, the first-order derivatives of the relative phase $\vartheta$ with respect to $\omega_1$ and $\omega_2$ determine the time delay \textit{map} between HH and VV possibilities for the signal and idler photons, respectively [that is, for a biphoton emission in the directions specified by $(\mathbf{q}_{1},-\mathbf{q}_{1})$]. Equation (\ref{eq:relative_phase_approx}) thus presents closed forms for the relative-phase and time-delay maps which give theoretical predictions equivalent to those of the iterative approach given in Ref. \cite{hegazy2017relative}. For convenience, the components of $\mathbf{q}_1$ are substituted by the free-space emission angles in $xz$ (horizontal) and $yz$ (vertical) planes as $q_{1x,y}= \tan \theta_{x1,y1}~\sqrt[]{(\omega_1/c)^2-|\mathbf{q}_1|^2}$ (similarly, $\mathbf{q}_2$).

\begin{figure*} [t]
   \begin{center}
   \begin{tabular}{c}
   \includegraphics[width=12cm]{setup_station6}
   \end{tabular}
   \end{center}
   \caption[example]
   { \label{fig:setup}
Extended entangled-photons source. The photon pairs are created in two crossed type-I 0.5-mm BBO crystals cut at 29.3$^\mathsf{o}$ and pumped by 405-nm beam to produce noncollinear SPDC photons centered at ~3$^\mathsf{o}$ half-opening angle. A half-wave plate (HWP$_p$) and a tiltable quarter-wave plate (QWP$_p$) manipulate the pump polarization state and are also components of polarization interferometer used for relative-phase measurements.
The relative-phase function is manipulated in frequency using a temporal compensator (TC), and in direction via a two-dimenstional spatial light modulator (2D SLM) charged by the inverted grayscale relative-phase map. PM: Partial mirror; D: Classical detector; APD: Avalanche photodiode, IF: interference filter, L: lens.
}
\end{figure*}

The coincidence map, on the other hand, is determined by the two-photon wavefunction at the detection plane \cite{klyshko1988photons}
\begin{equation}
\begin{split}
\Psi
=\langle 0 | E_1^{(+)}(\theta_{x1},\theta_{y1};t_1) E_2^{(+)}(\theta_{x2},\theta_{y2};t_2) |\psi \rangle, \\
\end{split}
\end{equation}
where $|0\rangle$ is the vacuum state and
\begin{equation*}
\begin{split}
& E_{1}^{(+)
=\int d\omega~e^{-i\omega t_{1}} \sum_{\mathbf{\sigma}}
(\mathbf{e}_1 . \mathbf{\sigma}) G(\omega)~ \mathbf{a}_\mathbf{\sigma}(\omega;\theta_{x1},\theta_{y1}), \\
& E_{2}^{(+)
=\int d\omega~ e^{-i\omega t_{2}}  \sum_{\mathbf{\sigma}}
(\mathbf{e}_2 . \mathbf{\sigma}) G(\omega)~ \mathbf{a}_\mathbf{\sigma}(\omega;\theta_{x2},\theta_{y2}),
\end{split}
\end{equation*}
are the positive-frequency parts of signal and idler field operators at the space-time coordinates of biphoton detection; $(\theta_{x1},\theta_{y1};t_{1})$ and $(\theta_{x2},\theta_{y2};t_{2})$ \cite{saleh2000duality}. Here the unit vector $\mathbf{e}_i=(e_{ix},e_{iy})$ specifies the orientation of the polarization analyzer in the $i$th SPDC arm, $G(.)$ is the transmissivity profile of the spectral filters, and $\mathbf{a}_\sigma(\omega; \theta_{x},\theta_{y})$ is the photon annihilation operator for a mode of frequency $\omega$, polarization $\mathbf{\sigma}$, and direction defined by the emission angles $(\theta_{x},\theta_{y})$ in free space. The spatiotemporal two-photon wavefunction thus writes \cite{atature2002multiparameter}
\begin{equation}
\begin{split}
\Psi
= \frac{1}{2} \mathbf{\int}&d\omega _{1}d\omega _{2}
  e^{-i(\omega_1 t_1 + \omega_2 t_2)}
\\
\times & G(\omega_1) G(\omega_2) \left( e_{1x}e_{2x}\Phi _\mathrm{HH}+ e_{1y}e_{2y} \Phi _\mathrm{VV} \right),
\end{split}
\label{Eq:DetectedWavefn}
\end{equation}
and the probability that the signal photon is captured by the typically used slow-response detector in the direction $(\theta_{x1},\theta_{y1})$ is then given by
\begin{equation}
P(\theta_{x1},\theta_{y1})=\int dt_1 dt_2 d\theta_{x2} d\theta_{y2}|\Psi (\theta_{x1},\theta_{y1},\theta_{x2},\theta_{y2}; t_1, t_2)|^2.
\label{Eq:angdistribution}
\end{equation}

 The schematic of the experimental setup is shown in Fig. \ref{fig:setup}. A collimated 405-nm beam emitted by a continuous-wave diode laser with $\sim$300-fs coherence time is used to pump two abutted 0.5-mm type-I BBO crystals with the optic axis of the first (second) crystal lying in the vertical (horizontal) plane. The two crystals are anti-reflection-coated at all wavelengths of concern and cut at 29.3$^\mathsf{o}$ to produce degenerate SPDC cone centered at ~3$^\mathsf{o}$ half-opening angle.
 The polarizing beam splitter (PBS) improves the extinction ratio of the pump beam (from 20:1 to about 10000:1). The half-wave plate (HWP$_p$) and tiltable quarter-wave plate (QWP$_p$) manipulate the pump polarization state –by rotating the linear polarization and adding phase difference $\phi_p$, respectively– and are also used to conduct the phase measurements as follows.
 To determine the phase difference $\phi_p$, the partial mirror (PM) back-reflects a tiny portion ($\sim 2\%$) of the pump beam to be analyzed by the dash-outlined polarization interferometer.
 At its output, the normalized optical power measured by the classical detector (D) is proportionate to $\sin^2 \phi_p$ (see supplementary material). We used 10-nm interference filters centered at 810 nm. This narrowband detection at the degenerate frequency along with the collimated-beam pumping render the approximations: $\theta_{x1}\approx - \theta_{x2}$ and $\theta_{y1}\approx - \theta_{y2}$ reasonable. We will henceforth refer to $\theta_{x1}$ and $\theta_{y1}$ as  $\theta_{x}$ and $\theta_{y}$, respectively, for short.

\begin{figure*} [t]

   \hspace{- .5cm}
   \centering
   \includegraphics[width=17cm]{Exp_ang_coincidence_two0-5BBO_2_polari_11-12}
   \caption[example]
   { \label{fig:coincidence}

Coincidence and time-delay maps all over the SPDC cone.
(a) Theoretical prediction (left) and experimental measurements (right) for the angular distribution of biphotons coincidence counts after polarization projection into the diagonal basis state $|++\rangle$. The right inset shows the normalized coincidence rates (red points), experimentally measured behind a 2-mm iris positioned, 50 cm from the SPDC crystals, at the directions: $\theta_x=\{2.4^\mathsf{o}, 2.7^\mathsf{o}, 3^\mathsf{o}, 3.3^\mathsf{o}, 3.45^\mathsf{o}, 3.6^\mathsf{o}\};$ with $\theta_y=0$.  The yellow line shows the normalized rates as predicted by Eq. (\ref{Eq:angdistribution}). (b) The same as in (a) after polarization projection into the rectilinear basis state $|HH\rangle$.
(c) The directional time-delay map between the HH and VV possibilities for the degenerate biphoton emission as predicted by the derivative of the relative phase $\vartheta$ with respect to $\omega_1$.}
   \end{figure*}

In the experiment, to figure out the spatial distribution of the scattered biphotons, the coincidence counts are measured behind a 2-mm circular iris translated laterally at $\sim$50 cm from the SPDC crystals. The measurements are made at a number of transverse positions which represent one radial group.
 Figure \ref{fig:coincidence}(a) shows the coincidence counts recorded after polarization projection into the basis state $|++\rangle$. We noticed that this spatial distribution varied when polarization projection was performed with different orientation. Figure \ref{fig:coincidence}(b) shows the coincidence counts after polarization projection into the basis state $|HH\rangle$ which exhibits more angular spread than that of the diagonal projection. This observation is in good agreement with Eq. (\ref{Eq:DetectedWavefn}) and Eq. (\ref{Eq:angdistribution}), where for the projector $|HH\rangle\langle HH|$, $e_{1x}e_{2x}=1$ and $e_{1y}e_{2y}=0$, which yields angular distribution $P(\theta_{x},\theta_{y})\propto \mathrm{sinc}^2 \left(\Delta \kappa^{ooe} L/2\right)$. Alternatively, for the diagonal projector $|++\rangle\langle ++|$, the factors $e_{1x}e_{2x}=e_{1y}e_{2y}=1/2$, therefore $P(\theta_{x},\theta_{y})\propto \mathrm{sinc}^2 \left(\Delta \kappa^{ooe} L/2\right) \cos^2 (\vartheta/2)$. The more confined angular distribution in case of diagonal projection originates from the interference of SPDC biphotons created along the two crystals. This observation shows one side of the coupling effects between polarization and spatial DoFs which can be eliminated via a suitable spatial phase compensation.

Figure \ref{fig:coincidence}(c) and Figure \ref{fig:phasemap}(a) show the time-delay and relative-phase maps all over the SPDC scattering cone as obtained by Eq. (\ref{eq:relative_phase_approx}) for the parameters of the experiment.
It can be noticed that while the time-delay map is essentially flat (centered at $\sim$200 fs), the relative-phase map has strong quadratic dependence in the radial direction.
This time delay is about the same as the coherence time of the SPDC light (determined by the interference filters). Therefore, over that wide directional extent, the temporal distinguishability between the leading (HH) and the lagging (VV) possibilities can be effectively washed out following the conventional approach; adding an offset temporal delay between the parent pump components $|V_p\rangle$ and  $|H_p\rangle$ \cite{nambu2002generation,rangarajan2009optimizing,hegazy2015tunable}.
In the experiment, this is accomplished using a temporal compensator (TC): Two \textit{available} 0.8-mm 30$^\mathsf{o}$-cut BBO crystals set abutted to each other with their optic axes rotated up-side-down  with respect to each other in the vertical plane; thereby canceling the transverse walk-off. To verify the degree of temporal compensation, we measured the visibility of coincidence fringes after 2-mm iris while HWP$_1$ is rotated and HWP$_2$ remains fixed at $22.5^\mathsf{o}$.
High polarization visibility was observed at all the transverse positions in Fig. \ref{fig:coincidence}; implying the flatness of the time-delay map and the effectiveness of the unified time compensation.

 At each observation point, the high polarization visibility was achievable at a tilt angle of QWP$_p$ different from its neighboring points (thus, different $\phi_p$). This indicates that the phase-matching part of the relative phase in Eq. (\ref{eq:relative_phase}) varies from a point to another.

 \begin{figure} [t]
   \begin{center}
   \begin{tabular}{c}
   \includegraphics[width=8.5cm]{phase_theta_pe_0_phi_pe_0_shortened_6}
   \end{tabular}
   \end{center}
   \caption[example]
   { \label{fig:phasemap}
(a) The two-dimensional quadratic relative-phase map (in radians) as predicted by Eq. (\ref{eq:relative_phase_approx}) all over the SPDC cone.  (b) Purifying SLM pattern to produce the states $|\phi^\pm\rangle$ with high fidelity over wide emission angles. The pattern constitutes the inverted \textit{modulo-2$\pi$} grayscale of the map in (a). The white rectangle circumscribes the 2D SLM's active area on which the compensation image is loaded. This area covers azimuthal-angle range $\sim$57$^\mathsf{o}$ (of the 180$^\mathsf{o}$ representing half SPDC cone); therefore about one  third  of  the  noncollinear SPDC emission cone could be manipulated.  (c) Experimentally measured phase retardation introduced by the 8-bit SLM at 810 nm plotted versus the gray display level (black:0 and white:255). The phase measurement error ($\pm 0.005 \pi$) is smaller than the readout markers.
}
   \end{figure}

Based on this observation, we present here an experimental technique to directly measure the relative phase of the entangled state:
First, polarization measurements are prepared projecting the two-photon state to diagonal basis states $|++\rangle, $ $|--\rangle$. We then scan the tilt angle of QWP$_p$. Along the scan, reaching the maximum coincidence counts indicates constructing a state nearest to the Bell state: $|\phi^+\rangle= (|H_1H_2\rangle+|V_1V_2\rangle)/\sqrt{2}$ [Here, the zero relative phase indicates that the phase-matching part of the relative phase is complementary to $\phi_p$ as dictated by Eq. (\ref{eq:relative_phase})] \cite{note:auto}. This determines the phase-matching part of the relative phase, since $\phi_p$ value is obtained from the normalized power detected by the photodetector D as aforementioned (see supplementary material).

To compensate for the relative-phase variation, a phase-only two-dimensional spatial light modulator (2D SLM) in one of the SPDC arms is used to introduce a pixel-based phase retardation between the H and V polarization components.
This reflective 2D SLM (display area 15.36 mm$\times$8.64 mm and resolution 1920$\times$1080 pixels) is placed at the transverse plane, $\sim$24.4  cm  from  the  SPDC  crystals, and programmed by a personal computer (PC). After loading by a grayscale image, the SLM develops a phase retardation for each pixel proportional to the pixel's gray level as commonly used in the monochromatic digital display.
The SLM can introduce up to $\sim 2\pi$ phase shift at 810 nm with its  8-bit display offering 256 intermediate gray/phase levels (black:0 and white:255).
 To treat the relative-phase map, the SLM was charged by the image within the white rectangle in Fig. \ref{fig:phasemap}(b) which is the inverted grayscale image of the modulo-2$\pi$ relative-phase map in Fig. \ref{fig:phasemap}(a). Figure \ref{fig:phasemap}(c) shows the experimentally measured phase of the SLM for the display gray levels \{0 (black), 63, 127, 192, 255 (white)\}.
 A special gamma correction is applied to enhance the linearity of the phase retardation to the pixel gray level at 810 nm.
Figure \ref{fig:measuredPhase_gradient} depicts the experimentally measured relative phase at the transverse positions in Fig. \ref{fig:coincidence} before and after charging the SLM with the compensating image.
The improved flatness of the directional relative-phase function implies less coupling between polarization and spatial DoFs. This yields polarization-entangled state $|\phi^+\rangle$ with high fidelity verified by $\sim 97\%$ polarization visibility measured at all transverse positions in Fig. \ref{fig:coincidence} across the SPDC emission (see supplementary material). This high visibility of coincidence fringes was observed at each  position behind the 2-mm iris while HWP$_1$ was rotated, with HWP$_2$ fixed at $22.5^\mathsf{o}$, and QWP$_p$ set at fixed tilting. By tilting QWP$_p$, the extended high-fidelity emission of the state $|\phi^-\rangle$ could also be formed.
   \begin{figure} [t]
   \begin{center}
   \begin{tabular}{c}
   \includegraphics[width=8.5 cm]{Experimental_Phase_both_5}
   \end{tabular}
   \end{center}
   \caption[example]
   { \label{fig:measuredPhase_gradient}
Experimentally measured relative phase for the produced polarization-entangled state across the SPDC cone before (brown rhombuses) and after (blue squares) loading the SLM. The measurements are done behind a 2-mm iris positioned, 50 cm from the SPDC crystals, at the directions: $\theta_x=\{2.7^\mathsf{o}, 3^\mathsf{o}, 3.3^\mathsf{o}, 3.45^\mathsf{o}, 3.6^\mathsf{o}\};$ with $\theta_y=0$. The error ($\pm 0.005 \pi$) is smaller than the symbols. Curves are fits to the data.
}
   \end{figure}

Beside equalizing the quadratic relative phase, the PC-programmed 2D SLM can simultaneously superimpose phase patterns serving other DoFs (e.g., orbital angular momentum modulation/demodulation). However, the following factors have to be taken into account:
(i) The SLM introduces unavoidable loss factor that can be classified as bulk losses due to the non-ideal reflectivity/transmissivity and diffraction losses due to the pixelation splitting of the incident light into several diffraction orders. In our experiment, the SLM losses reached $\sim 40 \%$. (ii) Effective anti-reflection coating is essential to diminish the possibility that the photon is not manipulated (by reflecting off the first interface) or manipulated multiple times (by rebouncing across the SLM).

To sum up, we have presented an experimental method to compensate for the spatial-spectral relative-phase function of a polarization entangled state, and produce the maximally entangled states $|\phi^\pm\rangle$ with high fidelity over wide emission angle –involving one third of the SPDC cone.
 In frequency domain, the relative-phase gradient has been treated conventionally by the use of unified temporal compensator.
In momentum domain, the relative phase over extended emission angles exhibits obvious two-dimensional quadratic dependence which has been equalized by using a high resolution PC-programmed 2D SLM loaded by the inverted relative-phase function. We have also developed a simple, yet accurate, experimental approach to measuring the relative phase of the entangled photon pairs by witnessing a normalized portion of the emerging power of the pump beam after polarization interferometry.
The high visibility and flat relative phase of the constructed state has been verified in a set of lateral positions across the SPDC cone. The presented approach can be applied to equalize the entire SPDC emission by using 2D SLM of sufficiently wide active area covering half the SPDC emission cone. Furthermore, the Bell states $|\psi^\pm\rangle$ –extended over wide angles– can be straightforwardly produced by inserting a half-wave plate \cite{note:waveplate} into one SPDC arm.
\\

See supplementary material for further information
about the relative-phase measurement technique and visibility measurements.
\\

This work is supported by ITAC program (Grant No. CFP134), ITIDA, Ministry of Communications and Information Technology (MCIT), Egypt.\\
\\
\noindent\textbf{Data Availability}

The data that supports the findings of this study are available within the article and its supplementary material.

\section
In the spontaneous parametric downconversion (SPDC) process, portion of pump photons split into photon pairs created with energies strictly satisfying the energy conservation condition and scattered in directions promoted by the momentum conservation (or phase-matching) condition.
The SPDC photons are thus entangled by nature in energy-time \cite{hong1987measurement,franson1989bell} as well as in momentum-position \cite{mair2001entanglement,lima2006propagation,barbieri2007complete, walborn2010spatial,leach2010quantum}. However for the polarization, the remaining degree of freedom (DoF), its entanglement appears always as a consequence of subtle designs using either a single \cite{kwiat1995new,Hugues2006non,kim2006phase,ueno2012entangled} or more nonlinear crystals \cite{kwiat1999ultrabright,kim2001interferometric,ljunggren2006theory,hegazy2017orthogonal,Villar2018Experimental}.

In either case, the wavefunctions of the two polarization possibilities (HH and VV for $|\phi^\pm\rangle$, or HV and VH for $|\psi^\pm\rangle$) can be made almost identical in magnitude, while a comprehensive picture about the indistinguishability in the space and time domains is provided by the directional-spectral relative-phase function. Generally, a flatter relative-phase function in some domain implies higher indistinguishability in its conjugate Fourier-transform domain
 \cite{hegazy2017orthogonal,Hegazy2017CLEO}.
Therefore, the spatial and temporal compensations actually correct (or flatten) the varying relative-phase function in direction and frequency. Several works previously addressed spatial and temporal compensations at specific spatial modes \cite{kwiat1995new,nambu2002generation,ljunggren2006theory,trojek2008collinear} or over relatively narrow aperture \cite{altepeter2005phase,rangarajan2009optimizing,cialdi2010programmable} by the use of additional birefringent elements.
Over such limited spatial extents, well-designed birefringent elements can effectively clean for the linear dependence of the relative phase function and restore the indistinguishability \cite{hegazy2015tunable}.
This is not the case when photon pairs are collected over relatively wide angular extents, where higher-order contributions of the relative-phase function drastically
reduce the effectiveness of linear methods.

In this paper, we present and demonstrate experimentally a method for eliminating the spatial and temporal distinguishability of a polarization-entangled state over wide emission angles.
We initially derive the two-photon state emitted by a noncollinear SPDC; including the directional-spectral relative-phase function all over the SPDC cone.
It is noticed that while the time-delay map (the spatial distribution of spectral relative-phase gradient) has a negligible change over wide angular extent, the directional relative-phase exhibits obvious quadratic profile. To restore the temporal and spatial indistinguishability, the compensation of the former can be effectively done in the conventional way; using a birefringent element. The latter is manipulated using a two-dimensional spatial light modulator (2D SLM) loaded with the inverted relative-phase profile. A fast yet accurate experimental technique is introduced to directly measure the relative-phase profile; dispensing with the quantum state tomography.

Consider a diagonally polarized pump beam illuminating a couple of crossed type-I crystals which are identically cut to promote production of degenerate photons in a noncollinear SPDC geometry.
The pump beam is classically treated as a superposition of frequency and plane-wave  modes:
$E_{{p}}(\mathbf{x},z,t)= \int d\omega _{{p}}d\mathbf{q}_{{p}}~A_{p}(\omega_{{p}};\mathbf{q} _{{p}})\exp i\left(\kappa _{{p}}z+\mathbf{q}_{{p}}. \mathbf{x}-\omega _{{p}}t\right) +\mathrm{c.c.}$,
with angular frequency $\omega _{p}$ and wavevector $\mathbf{k}_{p}=(%
\mathbf{q}_{{p}},\kappa _{{p}})$ where $\mathbf{q}_{{p}%
} $ are the transverse components along the coordinates $\mathbf{x}=(x,y)$, and $\kappa _{p}$ is the longitudinal component.
The pump polarization components in horizontal/vertical (H/V) basis have a relative phase $\phi_p$ that can be manipulated by tilting a birefringent element; e.g., a quarter wave plate (QWP).

Similarly, the produced SPDC two-photon emission can be expressed –based on the SPDC creation operators– as a joint spectral and spatial expansion in monochromatic planar waves with the angular frequencies $\omega_{1,2}$ and the transverse wave vectors  $\mathbf{q}_{\mathrm{1,2}}$ (subscripts 1, 2 denote signal and idler photons, respectively).
The longitudinal wavevectors of the interacting waves can be expressed as dictated by Maxwell's equations for ordinary (\textit{o}) and extraordinary (\textit{e}) polarization as
\cite{born2013principles,walborn2010spatial,hegazy2017orthogonal}
\begin{equation}
\begin{split}
&\kappa _{{j}}^{o }(\omega_{{j}};\mathbf{q}_{{j}})=\sqrt{( {\omega_{{j}} ~n_{{j}}^{o}}/{c}) ^{2}-\vert \mathbf{q}_{{j}}\vert
^{2}}\\
&\kappa _j ^{e(H,V)}(\omega_j;\mathbf{q}_j)
=   q_{j(x,y)} \tan\rho^\perp_j \\
&~~~~~~~~~~~~+ n_j^{e^\perp} \left[ \left(\tfrac{\omega _j}{c} \right) ^{2} -         \tfrac{1}{(n_j^e)^2}   q_{j(y,x)}^2-       (\tfrac{{n_j^{e^\perp}}}{n_j^e n_j^o} )^2         q_{j(x,y)}^2 \right]^{\frac{1}{2}},\\
\end{split}
\label{eq:Longwavevector}
\end{equation}
where the subscript $j=p, 1, 2$, the superscripts $(H)$ and $(V)$ label the SPDC crystals with optic axes lying in the $xz$ (horizontal) and $yz$ (vertical) plane respectively, $c$ is the speed of light in space, $n_{{j}}^{o}$ and $n_{{j}}^{e}$ are principal values of refractive index, and $n_j^{e^\perp}$ and $\rho^\perp_j$ are  the refractive indices and walk-off angles of the extraordinary-polarized ray propagating along $z$ axis (perpendicular to the crystals interfaces).
The SPDC in the two crystals creates a two-photon state entangled in every degree of freedom; frequency, momentum, and polarization, which is expressed at small emission angles  \cite{migdall1997polarization} by the superposition
\begin{equation}
\begin{split}
\left\vert \psi \right\rangle   \sim  \mathbf{\int}d\omega _{1}d\omega _{2}d\mathbf{q}_{1}d%
\mathbf{q}_{2} &{\large \{}\Phi _\mathrm{HH}\left(\omega _{1},\omega _{2};\mathbf{q}%
_{1},\mathbf{q}_{2}\right)\left\vert H_1H_2
\right\rangle \\
&+~\Phi _\mathrm{VV}\left( \omega _{1},\omega
_{2};\mathbf{q}_{1},\mathbf{q}_{2}\right) \left\vert V_1V_2
\right\rangle{\large \}},
\end{split}
\label{ent1}
\end{equation}
where $\Phi _\mathrm{HH,VV}\left(\omega _{1},\omega _{2};\mathbf{q}_{1},\mathbf{q}_{2}\right)$ are the biphoton wavefunctions corresponding to the HH and VV possibilities.
 Assuming that the crystals are of infinite transverse extent and pumped by a non-depleted beam with the reflected waves at all crystal interfaces being negligible, the two-photon wavefunctions can be written as
\begin{equation*}
\begin{split}
\Phi_\mathrm{HH} = ~&\tfrac{\chi^{(2)}}{\sqrt[]{2}} A_p(\omega_1+\omega_2;\mathbf{q}_1+\mathbf{q}_2) \int_{-2L}^{-L} dz~ e^{i \int_0^z dz^{\prime } \Delta \kappa(z^{\prime })}, \\
\propto ~&
e^{-i L\left(\tfrac{1}{2}{\Delta \kappa_{V}^{ooe}} + \Delta \kappa_H^{eeo}  \right)}
L
A_p(\omega_1+\omega_2;\mathbf{q}_1+\mathbf{q}_2) ~ \mathrm{sinc} \left( \tfrac{ \Delta \kappa_{V}^{ooe} L}{2}\right), \\
\Phi_\mathrm{VV} = ~&  e^{-i\phi_p}~ \tfrac{\chi^{(2)}}{\sqrt[]{2}}
A_p(\omega_1+\omega_2;\mathbf{q}_1+\mathbf{q}_2) \int_{-L}^{0} dz~ e^{i \int_0^z dz^{\prime } \Delta \kappa(z^{\prime })}, \\
\propto ~&
 e^{-i\left(\phi_p+\tfrac{1}{2}{\Delta \kappa_{H}^{ooe}L} \right)}
L
A_p(\omega_1+\omega_2;\mathbf{q}_1+\mathbf{q}_2)  ~\mathrm{sinc} \left(\tfrac{\Delta \kappa_{H}^{ooe} L}{2}\right), \\
\end{split}
\label{Eq: amplitudes}
\end{equation*}
where
$L$ is the crystal length, $\chi ^{(2)}$ is the bilinear susceptibility,
$~\Delta \kappa^{ooe}_{H,V} =\kappa _{p}^{e(H,V)}-\kappa _{1}^{o}-\kappa _{2}^{o}$ are  wavevector mismatches  within the interacting $(H)$ and $(V)$ crystals,
$\Delta \kappa_H^{eeo}=\kappa _{p}^{o}-\kappa _{1}^{e(H)}-\kappa_{2}^{e(H)}$ is the mismatch of waves interacting in the first crystal when passing through the second one ($\Delta \kappa_H^{eeo} \gg \Delta \kappa^{ooe}_{H,V}$, therefore no downconversion is considered in the later case).

The relative phase of the produced state $|\psi\rangle$ is then
\begin{equation}
\begin{split}
\vartheta
=\mathsf{arg}\left\{  \Phi_\mathrm{VV} / \Phi_\mathrm{HH} \right\}=
\left[\tfrac{1}{2}({\Delta \kappa_{V}^{ooe} }-\Delta \kappa_{H}^{ooe}) + \Delta \kappa_H^{eeo}  \right] L -\phi_p,
\end{split}
\label{eq:relative_phase}
\end{equation}
which includes phase-matching terms and an initial phase term. We assume that the spatiotemporal amplitude of the pump beam is factorizable:
$A_p(\omega_p;\mathbf{q}_p)= ({\sqrt[]{\tau_p}}/{\sqrt[4]{\pi}})\,\exp [-\tfrac{1}{2} \tau_p^2 (\omega_{{p}}-\omega_p^0)^2 ] \,\,\, \delta(\mathbf{q}_p),$ which expresses a polychromatic plane wave with a coherence time $\tau_p$. Therefore, the biphoton spatial characteristics can be fully determined with reference to the signal photon alone ($\mathbf{q}_2 = -\mathbf{q}_1$). This assumption implies also that $\Delta \kappa_{V}^{ooe} \approx\Delta \kappa_{H}^{ooe}\equiv\Delta \kappa^{ooe}$ as dictated by Eq. (\ref{eq:Longwavevector}); and consequently $|\Phi_\mathrm{HH}|\approx|\Phi_\mathrm{VV}|$. By substituting Eq. (\ref{eq:Longwavevector}) into Eq. (\ref{eq:relative_phase}) and ignoring the insignificant higher-order terms, the relative phase writes
\begin{equation}
\begin{split}
    \vartheta  (\omega _{1},\omega _{2};  \mathbf{q}_{1},  & -\mathbf{q}_{1}) \approx\\
\Big \{ & \frac{(n_p^o - n_1^{e^\perp}) \omega_1 + (n_p^o -  n_2^{e^\perp}) \omega_2}{c}\\
&- q_{1,x} (\tan \rho_1^\perp   -   \tan \rho_2^\perp ) \\
&+  \frac{c q_{1,x}^2}{2}    \Big[ \frac{(n_1^{e^\perp})^3}{\omega_1  (n_1^e )^2 (n_1^o)^2}  + \frac{(n_2^{e^\perp} )^3}{\omega_2  (n_2^e )^2 (n_2^o )^2}\Big] \\
&+ \frac{c q_{1,y}^2}{2} \Big[ \frac{n_1^{e^\perp}}{\omega_1 (n_1^e )^2} + \frac{n_2^{e^\perp}}{\omega_2 (n_2^e )^2 }\Big] \Big\} L - \phi_p,
\end{split}
\label{eq:relative_phase_approx}
\end{equation}
which is a quadratic function in the transverse wavevector components of the SPDC photons.
 Moreover, the first-order derivatives of the relative phase $\vartheta$ with respect to $\omega_1$ and $\omega_2$ determine the time delay \textit{map} between HH and VV possibilities for the signal and idler photons, respectively [that is, for a biphoton emission in the directions specified by $(\mathbf{q}_{1},-\mathbf{q}_{1})$]. Equation (\ref{eq:relative_phase_approx}) thus presents closed forms for the relative-phase and time-delay maps which give theoretical predictions equivalent to those of the iterative approach given in Ref. \cite{hegazy2017relative}. For convenience, the components of $\mathbf{q}_1$ are substituted by the free-space emission angles in $xz$ (horizontal) and $yz$ (vertical) planes as $q_{1x,y}= \tan \theta_{x1,y1}~\sqrt[]{(\omega_1/c)^2-|\mathbf{q}_1|^2}$ (similarly, $\mathbf{q}_2$).

\begin{figure*} [t]
   \begin{center}
   \begin{tabular}{c}
   \includegraphics[width=12cm]{setup_station6}
   \end{tabular}
   \end{center}
   \caption[example]
   { \label{fig:setup}
Extended entangled-photons source. The photon pairs are created in two crossed type-I 0.5-mm BBO crystals cut at 29.3$^\mathsf{o}$ and pumped by 405-nm beam to produce noncollinear SPDC photons centered at ~3$^\mathsf{o}$ half-opening angle. A half-wave plate (HWP$_p$) and a tiltable quarter-wave plate (QWP$_p$) manipulate the pump polarization state and are also components of polarization interferometer used for relative-phase measurements.
The relative-phase function is manipulated in frequency using a temporal compensator (TC), and in direction via a two-dimenstional spatial light modulator (2D SLM) charged by the inverted grayscale relative-phase map. PM: Partial mirror; D: Classical detector; APD: Avalanche photodiode, IF: interference filter, L: lens.
}
\end{figure*}

The coincidence map, on the other hand, is determined by the two-photon wavefunction at the detection plane \cite{klyshko1988photons}
\begin{equation}
\begin{split}
\Psi
=\langle 0 | E_1^{(+)}(\theta_{x1},\theta_{y1};t_1) E_2^{(+)}(\theta_{x2},\theta_{y2};t_2) |\psi \rangle, \\
\end{split}
\end{equation}
where $|0\rangle$ is the vacuum state and
\begin{equation*}
\begin{split}
& E_{1}^{(+)
=\int d\omega~e^{-i\omega t_{1}} \sum_{\mathbf{\sigma}}
(\mathbf{e}_1 . \mathbf{\sigma}) G(\omega)~ \mathbf{a}_\mathbf{\sigma}(\omega;\theta_{x1},\theta_{y1}), \\
& E_{2}^{(+)
=\int d\omega~ e^{-i\omega t_{2}}  \sum_{\mathbf{\sigma}}
(\mathbf{e}_2 . \mathbf{\sigma}) G(\omega)~ \mathbf{a}_\mathbf{\sigma}(\omega;\theta_{x2},\theta_{y2}),
\end{split}
\end{equation*}
are the positive-frequency parts of signal and idler field operators at the space-time coordinates of biphoton detection; $(\theta_{x1},\theta_{y1};t_{1})$ and $(\theta_{x2},\theta_{y2};t_{2})$ \cite{saleh2000duality}. Here the unit vector $\mathbf{e}_i=(e_{ix},e_{iy})$ specifies the orientation of the polarization analyzer in the $i$th SPDC arm, $G(.)$ is the transmissivity profile of the spectral filters, and $\mathbf{a}_\sigma(\omega; \theta_{x},\theta_{y})$ is the photon annihilation operator for a mode of frequency $\omega$, polarization $\mathbf{\sigma}$, and direction defined by the emission angles $(\theta_{x},\theta_{y})$ in free space. The spatiotemporal two-photon wavefunction thus writes \cite{atature2002multiparameter}
\begin{equation}
\begin{split}
\Psi
= \frac{1}{2} \mathbf{\int}&d\omega _{1}d\omega _{2}
  e^{-i(\omega_1 t_1 + \omega_2 t_2)}
\\
\times & G(\omega_1) G(\omega_2) \left( e_{1x}e_{2x}\Phi _\mathrm{HH}+ e_{1y}e_{2y} \Phi _\mathrm{VV} \right),
\end{split}
\label{Eq:DetectedWavefn}
\end{equation}
and the probability that the signal photon is captured by the typically used slow-response detector in the direction $(\theta_{x1},\theta_{y1})$ is then given by
\begin{equation}
P(\theta_{x1},\theta_{y1})=\int dt_1 dt_2 d\theta_{x2} d\theta_{y2}|\Psi (\theta_{x1},\theta_{y1},\theta_{x2},\theta_{y2}; t_1, t_2)|^2.
\label{Eq:angdistribution}
\end{equation}

 The schematic of the experimental setup is shown in Fig. \ref{fig:setup}. A collimated 405-nm beam emitted by a continuous-wave diode laser with $\sim$300-fs coherence time is used to pump two abutted 0.5-mm type-I BBO crystals with the optic axis of the first (second) crystal lying in the vertical (horizontal) plane. The two crystals are anti-reflection-coated at all wavelengths of concern and cut at 29.3$^\mathsf{o}$ to produce degenerate SPDC cone centered at ~3$^\mathsf{o}$ half-opening angle.
 The polarizing beam splitter (PBS) improves the extinction ratio of the pump beam (from 20:1 to about 10000:1). The half-wave plate (HWP$_p$) and tiltable quarter-wave plate (QWP$_p$) manipulate the pump polarization state –by rotating the linear polarization and adding phase difference $\phi_p$, respectively– and are also used to conduct the phase measurements as follows.
 To determine the phase difference $\phi_p$, the partial mirror (PM) back-reflects a tiny portion ($\sim 2\%$) of the pump beam to be analyzed by the dash-outlined polarization interferometer.
 At its output, the normalized optical power measured by the classical detector (D) is proportionate to $\sin^2 \phi_p$ (see supplementary material). We used 10-nm interference filters centered at 810 nm. This narrowband detection at the degenerate frequency along with the collimated-beam pumping render the approximations: $\theta_{x1}\approx - \theta_{x2}$ and $\theta_{y1}\approx - \theta_{y2}$ reasonable. We will henceforth refer to $\theta_{x1}$ and $\theta_{y1}$ as  $\theta_{x}$ and $\theta_{y}$, respectively, for short.

\begin{figure*} [t]

   \hspace{- .5cm}
   \centering
   \includegraphics[width=17cm]{Exp_ang_coincidence_two0-5BBO_2_polari_11-12}
   \caption[example]
   { \label{fig:coincidence}

Coincidence and time-delay maps all over the SPDC cone.
(a) Theoretical prediction (left) and experimental measurements (right) for the angular distribution of biphotons coincidence counts after polarization projection into the diagonal basis state $|++\rangle$. The right inset shows the normalized coincidence rates (red points), experimentally measured behind a 2-mm iris positioned, 50 cm from the SPDC crystals, at the directions: $\theta_x=\{2.4^\mathsf{o}, 2.7^\mathsf{o}, 3^\mathsf{o}, 3.3^\mathsf{o}, 3.45^\mathsf{o}, 3.6^\mathsf{o}\};$ with $\theta_y=0$.  The yellow line shows the normalized rates as predicted by Eq. (\ref{Eq:angdistribution}). (b) The same as in (a) after polarization projection into the rectilinear basis state $|HH\rangle$.
(c) The directional time-delay map between the HH and VV possibilities for the degenerate biphoton emission as predicted by the derivative of the relative phase $\vartheta$ with respect to $\omega_1$.}
   \end{figure*}

In the experiment, to figure out the spatial distribution of the scattered biphotons, the coincidence counts are measured behind a 2-mm circular iris translated laterally at $\sim$50 cm from the SPDC crystals. The measurements are made at a number of transverse positions which represent one radial group.
 Figure \ref{fig:coincidence}(a) shows the coincidence counts recorded after polarization projection into the basis state $|++\rangle$. We noticed that this spatial distribution varied when polarization projection was performed with different orientation. Figure \ref{fig:coincidence}(b) shows the coincidence counts after polarization projection into the basis state $|HH\rangle$ which exhibits more angular spread than that of the diagonal projection. This observation is in good agreement with Eq. (\ref{Eq:DetectedWavefn}) and Eq. (\ref{Eq:angdistribution}), where for the projector $|HH\rangle\langle HH|$, $e_{1x}e_{2x}=1$ and $e_{1y}e_{2y}=0$, which yields angular distribution $P(\theta_{x},\theta_{y})\propto \mathrm{sinc}^2 \left(\Delta \kappa^{ooe} L/2\right)$. Alternatively, for the diagonal projector $|++\rangle\langle ++|$, the factors $e_{1x}e_{2x}=e_{1y}e_{2y}=1/2$, therefore $P(\theta_{x},\theta_{y})\propto \mathrm{sinc}^2 \left(\Delta \kappa^{ooe} L/2\right) \cos^2 (\vartheta/2)$. The more confined angular distribution in case of diagonal projection originates from the interference of SPDC biphotons created along the two crystals. This observation shows one side of the coupling effects between polarization and spatial DoFs which can be eliminated via a suitable spatial phase compensation.

Figure \ref{fig:coincidence}(c) and Figure \ref{fig:phasemap}(a) show the time-delay and relative-phase maps all over the SPDC scattering cone as obtained by Eq. (\ref{eq:relative_phase_approx}) for the parameters of the experiment.
It can be noticed that while the time-delay map is essentially flat (centered at $\sim$200 fs), the relative-phase map has strong quadratic dependence in the radial direction.
This time delay is about the same as the coherence time of the SPDC light (determined by the interference filters). Therefore, over that wide directional extent, the temporal distinguishability between the leading (HH) and the lagging (VV) possibilities can be effectively washed out following the conventional approach; adding an offset temporal delay between the parent pump components $|V_p\rangle$ and  $|H_p\rangle$ \cite{nambu2002generation,rangarajan2009optimizing,hegazy2015tunable}.
In the experiment, this is accomplished using a temporal compensator (TC): Two \textit{available} 0.8-mm 30$^\mathsf{o}$-cut BBO crystals set abutted to each other with their optic axes rotated up-side-down  with respect to each other in the vertical plane; thereby canceling the transverse walk-off. To verify the degree of temporal compensation, we measured the visibility of coincidence fringes after 2-mm iris while HWP$_1$ is rotated and HWP$_2$ remains fixed at $22.5^\mathsf{o}$.
High polarization visibility was observed at all the transverse positions in Fig. \ref{fig:coincidence}; implying the flatness of the time-delay map and the effectiveness of the unified time compensation.

 At each observation point, the high polarization visibility was achievable at a tilt angle of QWP$_p$ different from its neighboring points (thus, different $\phi_p$). This indicates that the phase-matching part of the relative phase in Eq. (\ref{eq:relative_phase}) varies from a point to another.

 \begin{figure} [t]
   \begin{center}
   \begin{tabular}{c}
   \includegraphics[width=8.5cm]{phase_theta_pe_0_phi_pe_0_shortened_6}
   \end{tabular}
   \end{center}
   \caption[example]
   { \label{fig:phasemap}
(a) The two-dimensional quadratic relative-phase map (in radians) as predicted by Eq. (\ref{eq:relative_phase_approx}) all over the SPDC cone.  (b) Purifying SLM pattern to produce the states $|\phi^\pm\rangle$ with high fidelity over wide emission angles. The pattern constitutes the inverted \textit{modulo-2$\pi$} grayscale of the map in (a). The white rectangle circumscribes the 2D SLM's active area on which the compensation image is loaded. This area covers azimuthal-angle range $\sim$57$^\mathsf{o}$ (of the 180$^\mathsf{o}$ representing half SPDC cone); therefore about one  third  of  the  noncollinear SPDC emission cone could be manipulated.  (c) Experimentally measured phase retardation introduced by the 8-bit SLM at 810 nm plotted versus the gray display level (black:0 and white:255). The phase measurement error ($\pm 0.005 \pi$) is smaller than the readout markers.
}
   \end{figure}

Based on this observation, we present here an experimental technique to directly measure the relative phase of the entangled state:
First, polarization measurements are prepared projecting the two-photon state to diagonal basis states $|++\rangle, $ $|--\rangle$. We then scan the tilt angle of QWP$_p$. Along the scan, reaching the maximum coincidence counts indicates constructing a state nearest to the Bell state: $|\phi^+\rangle= (|H_1H_2\rangle+|V_1V_2\rangle)/\sqrt{2}$ [Here, the zero relative phase indicates that the phase-matching part of the relative phase is complementary to $\phi_p$ as dictated by Eq. (\ref{eq:relative_phase})] \cite{note:auto}. This determines the phase-matching part of the relative phase, since $\phi_p$ value is obtained from the normalized power detected by the photodetector D as aforementioned (see supplementary material).

To compensate for the relative-phase variation, a phase-only two-dimensional spatial light modulator (2D SLM) in one of the SPDC arms is used to introduce a pixel-based phase retardation between the H and V polarization components.
This reflective 2D SLM (display area 15.36 mm$\times$8.64 mm and resolution 1920$\times$1080 pixels) is placed at the transverse plane, $\sim$24.4  cm  from  the  SPDC  crystals, and programmed by a personal computer (PC). After loading by a grayscale image, the SLM develops a phase retardation for each pixel proportional to the pixel's gray level as commonly used in the monochromatic digital display.
The SLM can introduce up to $\sim 2\pi$ phase shift at 810 nm with its  8-bit display offering 256 intermediate gray/phase levels (black:0 and white:255).
 To treat the relative-phase map, the SLM was charged by the image within the white rectangle in Fig. \ref{fig:phasemap}(b) which is the inverted grayscale image of the modulo-2$\pi$ relative-phase map in Fig. \ref{fig:phasemap}(a). Figure \ref{fig:phasemap}(c) shows the experimentally measured phase of the SLM for the display gray levels \{0 (black), 63, 127, 192, 255 (white)\}.
 A special gamma correction is applied to enhance the linearity of the phase retardation to the pixel gray level at 810 nm.
Figure \ref{fig:measuredPhase_gradient} depicts the experimentally measured relative phase at the transverse positions in Fig. \ref{fig:coincidence} before and after charging the SLM with the compensating image.
The improved flatness of the directional relative-phase function implies less coupling between polarization and spatial DoFs. This yields polarization-entangled state $|\phi^+\rangle$ with high fidelity verified by $\sim 97\%$ polarization visibility measured at all transverse positions in Fig. \ref{fig:coincidence} across the SPDC emission (see supplementary material). This high visibility of coincidence fringes was observed at each  position behind the 2-mm iris while HWP$_1$ was rotated, with HWP$_2$ fixed at $22.5^\mathsf{o}$, and QWP$_p$ set at fixed tilting. By tilting QWP$_p$, the extended high-fidelity emission of the state $|\phi^-\rangle$ could also be formed.
   \begin{figure} [t]
   \begin{center}
   \begin{tabular}{c}
   \includegraphics[width=8.5 cm]{Experimental_Phase_both_5}
   \end{tabular}
   \end{center}
   \caption[example]
   { \label{fig:measuredPhase_gradient}
Experimentally measured relative phase for the produced polarization-entangled state across the SPDC cone before (brown rhombuses) and after (blue squares) loading the SLM. The measurements are done behind a 2-mm iris positioned, 50 cm from the SPDC crystals, at the directions: $\theta_x=\{2.7^\mathsf{o}, 3^\mathsf{o}, 3.3^\mathsf{o}, 3.45^\mathsf{o}, 3.6^\mathsf{o}\};$ with $\theta_y=0$. The error ($\pm 0.005 \pi$) is smaller than the symbols. Curves are fits to the data.
}
   \end{figure}

Beside equalizing the quadratic relative phase, the PC-programmed 2D SLM can simultaneously superimpose phase patterns serving other DoFs (e.g., orbital angular momentum modulation/demodulation). However, the following factors have to be taken into account:
(i) The SLM introduces unavoidable loss factor that can be classified as bulk losses due to the non-ideal reflectivity/transmissivity and diffraction losses due to the pixelation splitting of the incident light into several diffraction orders. In our experiment, the SLM losses reached $\sim 40 \%$. (ii) Effective anti-reflection coating is essential to diminish the possibility that the photon is not manipulated (by reflecting off the first interface) or manipulated multiple times (by rebouncing across the SLM).

To sum up, we have presented an experimental method to compensate for the spatial-spectral relative-phase function of a polarization entangled state, and produce the maximally entangled states $|\phi^\pm\rangle$ with high fidelity over wide emission angle –involving one third of the SPDC cone.
 In frequency domain, the relative-phase gradient has been treated conventionally by the use of unified temporal compensator.
In momentum domain, the relative phase over extended emission angles exhibits obvious two-dimensional quadratic dependence which has been equalized by using a high resolution PC-programmed 2D SLM loaded by the inverted relative-phase function. We have also developed a simple, yet accurate, experimental approach to measuring the relative phase of the entangled photon pairs by witnessing a normalized portion of the emerging power of the pump beam after polarization interferometry.
The high visibility and flat relative phase of the constructed state has been verified in a set of lateral positions across the SPDC cone. The presented approach can be applied to equalize the entire SPDC emission by using 2D SLM of sufficiently wide active area covering half the SPDC emission cone. Furthermore, the Bell states $|\psi^\pm\rangle$ –extended over wide angles– can be straightforwardly produced by inserting a half-wave plate \cite{note:waveplate} into one SPDC arm.
\\

See supplementary material for further information
about the relative-phase measurement technique and visibility measurements.
\\

This work is supported by ITAC program (Grant No. CFP134), ITIDA, Ministry of Communications and Information Technology (MCIT), Egypt.\\
\\
\noindent\textbf{Data Availability}

The data that supports the findings of this study are available within the article and its supplementary material.


\section{Introduction}
\vspace{-1mm}

 Matrix factorization provides a concise representation of data. Despite desirable uniqueness conditions and computational simplicity of the well-known singular value decomposition (SVD), it comes with some fundamental shortcomings. The intrinsic structure of data is not inherited to the singular components.  Moreover, SVD implies orthogonality  on the components which is irrelevant to the underlying structure of the original data. This enforced  structure makes the bases, a.k.a. singular vectors, hard to interpret \cite{hamm2020perspectives}. On the other hand, it is shown that borrowing bases from the actual samples of a dataset provides a robust representation, which can be employed in interesting applications where sampling is  their key factor \cite{wang2013improving}. This problem is studied under the literature of column subset selection problem (CSSP) \cite{deshpande2010efficient,boutsidis2009improved} and CUR decomposition \cite{mahoney2009cur,boutsidis2017optimal}. A general problem for CSSP and CUR can be written in the following form:
\begin{equation}
\label{eq:main}
\small
    (\mathbb{S}^c,\mathbb{S}^r)=\underset{\mathbb{S}^c,\mathbb{S}^r}{\text{argmin}}\|{\boldsymbol{X}}-\pi_{\mathbb{S}^r}^{r}(\pi_{\mathbb{S}^c}^{c}({\boldsymbol{X}}))\|_F^2,
\end{equation}
where, $\!\small{\boldsymbol{X}}\!\!\in\!\! \mathbb{R}^{N\times M}\!$ is the data matrix containing $M$ data points in an $N$-dimensional space. Here, $\pi_{\mathbb{S}^c}^{c}(.)$ and $\pi_{\mathbb{S}^r}^{r}(.)$ indicate column space projection and row space projection, respectively. These operators project all columns (rows) to a low-dimensional subspace spanned by selected columns (rows) of matrix $\boldsymbol{X}$ indexed by the set $\mathbb{S}^c$ ($\mathbb{S}^r$). The chronological order in applying $\pi_{\mathbb{S}^c}^{c}(.)$ and $\pi_{\mathbb{S}^r}^{r}(.)$ does not affect the problem since these operators are linear. Moreover, substituting $\pi_{\mathbb{S}^r}^{r}$ with the identity projection simplifies the problem to CSSP. Fig. \ref{fig:unfold} illustrates the structure of CUR matrix decomposition as a self-representative approach.

A versatile metric for evaluating the performance of a data subset selection algorithm can be defined by the approximation error resulted from the projection of the entire data to the span of selected rows/columns.
How close to the optimal selection an algorithm can reach, is determined by comparing its approximation error to the best low-rank approximation error specified by the spectral decomposition.  Recently, we proposed a fast and accurate algorithm for solving CSSP which is called spectrum pursuit (SP) \cite{joneidi2020select}.
\begin{figure*}[h!]
    \centering
    \includegraphics[width=0.83\textwidth]{figs/decomposition.pdf}
    \vspace{-3mm}
    \caption{\small{Two columns and three rows from matrix $\boldsymbol{X}$ are selected and organized in matrix $\boldsymbol{C}$ and $\boldsymbol{R}$. The outer product of each pair of  a selected column and a selected row constructs a rank-1 matrix, i.e., $\boldsymbol{c}_i\boldsymbol{r}_j^T$. The  contribution amount of each pair is reflected in variable $u_{ij}$. The core matrix $\boldsymbol{U}$ is the collection of all $u_{ij}$'s. The goal is to minimize $\|\boldsymbol{X}-\hat{\boldsymbol{X}}\|_F$ where $\hat{\boldsymbol{X}}=\boldsymbol{CUR}$. } }
    \vspace{-4mm}
    \label{fig:unfold}
\end{figure*}
Inspired by SP, we propose a new algorithm to address the more general case of the CUR matrix decomposition.
The main contributions of our paper are summarized as follows:
\vspace{-2mm}
\begin{itemize}
\small
    \item A novel algorithm for CUR  decomposition, referred to as two-way spectrum pursuit (TWSP), is proposed. TWSP provides an accurate solution for CUR decomposition.
    \item TWSP  enjoys a linear complexity w.r.t. the number of columns and the number of rows of a matrix.
    \item TWSP is a parameter-free algorithm that only requires the number of desired columns and rows for selection. Thus, TWSP does not require any parameter fine-tuning.
    \item The  TWSP algorithm is put to the test and investigated in a set of synthetic and real experiments.
    \item The role of the core matrix $\boldsymbol{U}$ in $\boldsymbol{CUR}$ decomposition is illustrated which shows the connection between selected columns and rows. Based on analysis of $\boldsymbol{U}$, an interesting application for joint sensor/channel selection is presented.
\end{itemize}

\vspace{-2mm}
\textbf{Notations:} Throughout this paper, vectors, and matrices are denoted by bold lowercase, and bold uppercase letters, respectively. Moreover sets are indicated by blackboard bold characters. The $m^{\text{th}}$ column of matrix $\boldsymbol{X}$ is denoted by $\boldsymbol{X}(:,m)$ and the $n^{\text{th}}$ row of matrix $\boldsymbol{X}$ is denoted by $\boldsymbol{X}(n,:)$.

The rest of paper is organized as follows. Sec. \ref{sec:background} studies the problem of column subset selection as the conventional one-way selection scheme. In Sec. \ref{sec:twsp}, our main work is presented for joint column and row subset selection. Experimental results are exhibited in Sec. \ref{sec:exp}.
\section{Column Subset Selection Problem}

\label{sec:background}
Selecting the most diverse subset of data in an optimal sense is studied vastly \cite{deshpande2010efficient,li2017polynomial,joneidi2018optimal}. However, these methods do not guarantee that the unselected columns are  \emph{well represented} by the selected ones. Further, outliers are selected with a high probability using such algorithms due to their diversity \cite{joneidi2020select}.
 A more effective approach is selecting some \emph{representatives} which are able to approximate the rest of data accurately \cite{elhamifar2012see} as defined as a special case of \eqref{eq:main}.
This is an NP-hard problem \cite{shitov2020column} and there are several efforts for solving this  problem \cite{boutsidis2014near,deshpande2010efficient,paul2015column,zaeemzadeh2019iterative}. There are computationally expensive approaches based on convex relaxation \cite{elhamifar2012see,elhamifar2015dissimilarity} that are not computationally feasible for large datasets since their complexity is of order $O(M^3)$, where $M$ is the number of original columns. Recently, we proposed a new algorithm for solving CSSP with a linear complexity which is called Spectrum pursuit (SP)~\cite{joneidi2020select}. The SP algorithm finds $K$ columns of $\boldsymbol{X}$ such that their span is close to that of the best rank-$K$ approximation of $\boldsymbol{X}$. SP is an iterative approach where at each iteration one sample selection is optimized such that the ensemble of selected samples describes the whole dataset more accurately.  SP finds representatives such that the column space is spanned accurately via consecutive rank-$1$ approximations as theoretically analyzed in \cite{joneidi2020optimality}. In the present paper, we extend the SP algorithm for selecting columns and rows jointly such that their outer product can represent the whole matrix accurately. A naive approach is applying SP algorithm on the matrix of interest to select a subset of columns and applying SP on its transpose in order to find a subset of rows. However, this approach is not efficient and we will compare it with our proposed approach which is optimized  through a joint representation of selected columns and rows.
 \iffalse
\begin{algorithm}
\caption{Spectrum Pursuit Algorithm \cite{joneidi2020select}} \label{alg:SP}
\algsetup{
linenosize=\small,
linenodelimiter=:
}

\begin{algorithmic}[1]
\footnotesize
\REQUIRE $\boldsymbol{X}$ and $K$\\
\hspace{-3mm} \textbf{Output}:$\;\boldsymbol{X}_{\mathbb{S}}$\\
\STATE \textbf{Initialization:} \\
 $\mathbb{S}\leftarrow$A random subset of $\{1,\ldots,M\}$ with $|\mathbb{S}|=K$ \\
 $\{\mathbb{S}_k\}_{k=1}^K\!\!\leftarrow\!$ Partition $\mathbb{S}$ into $K$ subsets, each containing one element.\\
 iter$=0$\\ while a stopping criterion is not met
  \STATE $\quad$ $k=$mod(iter,K)+1

  \STATE $\quad$
$\boldsymbol{U}_{\overline{k}}=\;$normalize column($\boldsymbol{X}_{\mathbb{S}\backslash\mathbb{S}_k}$) \vspace{+0.2mm}
\STATE $\quad \boldsymbol{V}_{\overline{k}}=\boldsymbol{X}^T\boldsymbol{U}_{\overline{k}}
(\boldsymbol{U}^T_{\overline{k}} \boldsymbol{U}_{\overline{k}})^{-1}$
\STATE $\quad$  $\boldsymbol{E}_{\overline{k}}=\boldsymbol{X}-\boldsymbol{U}_{\overline{k}}{\boldsymbol{V}_{\overline{k}}}^T\;\;\;$ (null-space projection)
\STATE $\quad$ $\boldsymbol{u}_{k}=\;$find the first left singular vector of $\boldsymbol{E}_{\overline{k}}$
\STATE $\quad$ $\mathbb{S}_k\xleftarrow[]{}$ index of the most correlated column of $\boldsymbol{E}_{\overline{k}}$ with $\boldsymbol{u}_k$
\STATE $\quad$ $\mathbb{S}\xleftarrow[]{} \bigcup_{k'=1}^K \mathbb{S}_{k'}$
\STATE $\quad$ iter=iter+1
\\ end while
\end{algorithmic}
\label{alg1}
\vspace{-0.2mm}
\end{algorithm}
\vspace{-2mm}

\fi

\vspace{-3mm}
\section{Two-way Spectrum Pursuit }
\label{sec:twsp}
 The introduced joint column/row subset selection in \eqref{eq:main} can be written as a CUR decomposition in the following form in which factor matrices must be drawn from actual columns/rows of the original matrix as
\begin{align}
\small
\label{eq:problem_equ}
  (\boldsymbol{C},\boldsymbol{U},\boldsymbol{R})=\underset{\boldsymbol{C},\boldsymbol{U},\boldsymbol{R}}{\text{argmin}}\|{\boldsymbol{X}}-\boldsymbol{C}\boldsymbol{U}\boldsymbol{R}\|_F^2 , \\
  \text{s.t.}\;\;\boldsymbol{c}_k \in \mathbb{X}_c \;\text{and}\;\boldsymbol{r}_k\in\mathbb{X}_r.\nonumber
\end{align}
In this problem, $\mathbb{X}_c$, and $\mathbb{X}_r$ indicate the set of normalized columns and rows of matrix $\boldsymbol{X}$, respectively. Here, $\boldsymbol{c}_k$ and $\boldsymbol{r}^T_k$ denote the $k^{th}$ column and the $k^{th}$ row of $\boldsymbol{C}$ and $\boldsymbol{R}$, respectively. In other words, $\mathbb{X}_c\!=\!\{\boldsymbol{X}(:,m)/\|\boldsymbol{X}(:,m)\|\}$ for all columns and  $\mathbb{X}_r\!=\!\{\boldsymbol{X}(n,:)/\|\boldsymbol{X}(n,:)\|\}$ for all rows.
Please note that replacing constraints in \eqref{eq:problem_equ} with orthogonality constraint on $\boldsymbol{c}_k$'s and on $\boldsymbol{r}_k$'s results in the truncated singular value decomposition (SVD) with $K$ most significant components. In this case, $\boldsymbol{C}$  and $\boldsymbol{R}$ contain the first $K$ left singular vectors and the first $K$ right singular vectors of $\boldsymbol{X}$, respectively. Moreover, the core matrix will be diagonal and the entries will be singular values with diagonal entries as singular values. However, the underlying constraints in (\ref{eq:problem_equ}) turn the problem into a joint subset of row and column selection problem instead of matrix low-rank approximation problem.

To solve this complicated problem, we split it into two consecutive problems for optimization on the $k^{\text{th}}$ selected column/row. Our optimization approach is alternative, i.e., a random subset of columns and rows are picked. Then, one column or row is considered to be replaced with a more efficient one at each iteration. Since scaling  a vector does not change its span, without loss of generality we assume that the column or the row subject of the optimization lie on the unit sphere. At each iteration, a rank-$1$ component is optimized characterized by $\boldsymbol{c}\boldsymbol{g}^{T}$ or $\boldsymbol{h}\boldsymbol{r}^{T}$ given by

\begin{subequations}
\small
\vspace{-1mm}
\label{eq:seperated}
\begin{equation}
\label{eq:problem_sep1}
  \underset{\boldsymbol{c},\;\boldsymbol{W},\;\boldsymbol{g}}{\text{argmin}}\;\;\|\underbrace{\boldsymbol{X}-\boldsymbol{C}_k\boldsymbol{W}\boldsymbol{R}}_{\boldsymbol{E}^c}- \boldsymbol{c}\boldsymbol{g}^T\|_F^2 \;\text{s.t.}\;\;\|\boldsymbol{c}\|_2=1
  \end{equation}
  \begin{equation}
  \label{eq:problem_sep2}
   \underset{\boldsymbol{h},\;\boldsymbol{Y},\;\boldsymbol{r}}{\text{argmin}}\;\;\|\underbrace{\boldsymbol{X}-\boldsymbol{C}\boldsymbol{Y}\boldsymbol{R}_k}_{\boldsymbol{E}^r}- \boldsymbol{h}\boldsymbol{r}^T\|_F^2 \;\text{s.t.}\;\;\|\boldsymbol{r}\|_2=1
    \end{equation}
\end{subequations}

Matrix $\boldsymbol{C}_k$ is the set of selected columns except the $k^{\text{th}}$ one and $\boldsymbol{R}_k$ is the set of selected rows except the $k^{\text{th}}$ row.  The first subproblem can be solved easily w.r.t. $\boldsymbol{c}$ using singular value decomposition. In other words, $\boldsymbol{cg}^T$ and $\boldsymbol{hr}^T$ are the best rank-$1$ approximations of the residual $\boldsymbol{E}^c$ and $\boldsymbol{E}^r$, respectively. The obtained $\boldsymbol{c}$/$\boldsymbol{r}$ is the best column/row that can be added to the pool of selected columns/rows. However, the obtained vector is not available in the given dataset as a column or row since it is a singular vector which is a function of all columns/rows.  The following step re-imposes the underlying constraints at each iteration,

\begin{subequations}
\vspace{-4mm}
\small
\label{eq:matched}
\begin{equation}
\label{eq:problem_match1}
\mathbb{S}_k^c=\underset{m}{\text{argmax}}|\boldsymbol{x}_m^T\boldsymbol{c}|,\;\; \forall{\boldsymbol{x}_m}\in \mathbb{X}_c
  \end{equation}
  \begin{equation}
  \label{eq:problem_match2}
\mathbb{S}_k^r=\underset{n}{\text{argmax}}|\boldsymbol{x}_n^T\boldsymbol{r}|,\;\;\forall{\boldsymbol{x}_n}\in \mathbb{X}_r
    \end{equation}
\end{subequations}
Here, $\mathbb{S}_k^c$ indicates a singleton that contains $k^{\text{th}}$ selected column and $\mathbb{S}_k^r$ corresponds to the $k^{\text{th}}$ selected row. In each iteration, one column or one row is the subject of optimization. The impact of the latest estimation for that column/row on the representation is neglected. Then, an optimized replacement is found.
At each iteration of TWSP, sub-problems in \eqref{eq:seperated} are solved and their solutions are matched to the accessible column samples and row samples through matching equations in \eqref{eq:matched}. The new selected column or row is stored in $\mathbb{S}_k^c$ or $\mathbb{S}_k^r$, respectively. At each iteration we need to compute only the first singular vector and there are fast methods to do so \cite{comon1990tracking}. The pair of \eqref{eq:problem_sep1} and \eqref{eq:problem_match1} optimizes and matches a column. Similarly, performing \eqref{eq:problem_sep2} and \eqref{eq:problem_match2} provides us with an optimized row. However, we do not update both of them per each iteration. In fact, we need to perform a column update or a row update in each iteration. It should be determined that which update (column or row) is more efficient in each iteration of the algorithm.  To this aim, first, we choose a random previously selected column and a random previously selected row. Then, we find the best possible replacement column and the best possible replacement row.  Accordingly, we choose whichever that minimizes the cost function more. The best modified column-wise subset is denoted by $\tilde{\mathbb{S}}^c$ and $\tilde{\mathbb{S}}^r$ denotes the best row-wise modified subset. Alg. \ref{alg:jcrss} indicates the steps of TWSP algorithm.  Here, $\dagger$ refers to the Moore–Penrose pseudo-inverse operator. Iterations can be terminated either once CUR decomposition error is saturated or when a maximum number of iterations is reached.

\iffalse
\begin{algorithm}[h]
\caption{Two-way Spectrum Pursuit}
\label{alg:jcrss}
\begin{algorithmic}
\small{
\REQUIRE {\small{${\boldsymbol{X}}\!\in\! \mathbb{R}^{N\times M}\!$, $K_1$ and $K_2$.}}\\ \vspace{1.2mm}
\hspace{-3.3mm}\textbf{Output:} $\mathbb{S}^c$ and $\mathbb{S}^r$}.
\STATE \textbf{Initialization:} \\
 $\mathbb{S}^c\leftarrow$A random subset of $\{1,\ldots,M\}$ with $|\mathbb{S}^c|=K_1$ \\
  $\mathbb{S}^r\leftarrow$A random subset of $\{1,\ldots,N\}$ with $|\mathbb{S}^r|=K_2$ \\
 $\{\mathbb{S}^c_k\}_{k=1}^K\!\!\leftarrow\!$ Partition $\mathbb{S}^c$ into $K_1$ subsets.\\
  $\{\mathbb{S}^r_k\}_{k=1}^K\!\!\leftarrow\!$ Partition $\mathbb{S}^r$ into $K_2$ subsets\\
 iter$=0$\\
 while a stopping criterion is not met
  \STATE $\qquad$ $k=$mod(iter,$K_1$)$+1$\\
  $\qquad \mathbb{S}^c_{\overline{k}}=\mathbb{S}^c\backslash \mathbb{S}^c_k$\\
  $\qquad \boldsymbol{W}= \boldsymbol{C}_k^{\dagger}\boldsymbol{X}\boldsymbol{R}^{\dagger}$\\
  $\qquad \boldsymbol{E}^c=\boldsymbol{X}-\boldsymbol{C}_k\boldsymbol{WR}\;\;$ {\footnotesize (Null space projection for columns)}\\
  \STATE $\qquad$ $\boldsymbol{c}=\;$find the first left singular vector of $\boldsymbol{E}^c$ \eqref{eq:problem_sep1}\\
\STATE $\qquad$ $\mathbb{S}_k^c\xleftarrow[]{}$  the most correlated column of $\boldsymbol{E}$ with $\boldsymbol{c}$  \eqref{eq:problem_match1}
\STATE $\qquad$ $\mathbb{Z}^c\xleftarrow[]{} \bigcup_{k'=1}^{K_1} \mathbb{S}^c_{k'}$
  \vspace{1mm}
 \STATE $\qquad$ $k=$mod(iter,$K_2$)$+1$\\
  $\qquad \mathbb{S}^r_{\overline{k}}=\mathbb{S}^r\backslash \mathbb{S}^r_k$\\
    $\qquad \boldsymbol{Y}= \boldsymbol{C}^{\dagger}\boldsymbol{X}\boldsymbol{R}_k^{\dagger}$\\
  $\qquad \boldsymbol{E}^r=\boldsymbol{X}-\boldsymbol{CYR}_k \;\;${\footnotesize (Null space projection for rows)}\\
  \STATE $\qquad$ $\boldsymbol{r}=\;$find the first right singular vector of $\boldsymbol{E}^r$ \eqref{eq:problem_sep2} \\
\STATE $\qquad$ $\mathbb{S}_k^r\xleftarrow[]{}$ the most correlated row of $\boldsymbol{E}$ with $\boldsymbol{r}$ \eqref{eq:problem_match2}
\STATE $\qquad$ $\mathbb{S}^r\xleftarrow[]{} \bigcup_{k'=1}^{K_2} \mathbb{S}^r_{k'}$
\STATE $\qquad \text{iter}\xleftarrow[]{} \text{iter}+1$
\end{algorithmic}
\end{algorithm}

\fi

\begin{figure*}[t]
\label{fig:synthesis}
\begin{subfigure}{0.235\textwidth}
\centering
\includegraphics[width=1.45 in]{Figures/CUR.eps}
\footnotesize{\caption{Adaptive CUR \cite{mahoney2009cur}}}
\end{subfigure}
\begin{subfigure}{0.235\textwidth}
\centering
\includegraphics[width=1.45 in]{Figures/optimal.eps}
\footnotesize{\caption{Near-optimal \cite{wang2013improving}}}
\end{subfigure}
\begin{subfigure}{0.235\textwidth}
\centering
\includegraphics[width=1.45 in]{Figures/SP.eps}
\footnotesize{\caption{SP \cite{joneidi2020select}}}
\end{subfigure}
\begin{subfigure}{0.255\textwidth}
\centering
\includegraphics[width=1.63 in,height=1.03in]{Figures/main.eps}
\footnotesize{\caption{TWSP}}
\end{subfigure}

\footnotesize{\caption{(a)-(d) Performance comparison in terms of the normalized error of CUR decomposition.}}
\vspace{-1mm}
\end{figure*}

\begin{figure}[t]
\vspace{-3mm}
\centering
\includegraphics[width=2.9in, height=1.55in]{Figures/init_condition.eps}
  \caption{\footnotesize{The convergence behavior of the proposed algorithm w.r.t. the initial condition of selected subset. The initial cost function is corresponding to the initial set which is drawn randomly. The blue curves indicate the optimization path alongside iterations of the TWSP algorithm. Here, $100$ different realizations are studied.}}
  \vspace{-4mm}
\label{fig:init_convrgn}
\end{figure}

The proposed TWSP provides the selected columns and selected rows in order to form matrix $\boldsymbol{C}$ and $\boldsymbol{R}$ in CUR decomposition. It is straightforward to estimate the core matrix $\boldsymbol{U}$. Mathematically,
\vspace{-1mm}
\begin{equation}
\label{eq:core}
    \boldsymbol{U}=\boldsymbol{C}^{\dagger} \boldsymbol{X}\boldsymbol{R}^{\dagger.}
\vspace{-1mm}
\end{equation}
This matrix is a two-way compressed replica of the whole dataset and it contains valuable information in practice as will be discussed in Sec. \ref{subsec:joint}. In general, the number of selected columns may differ from the desired number of rows. Here, $K_1$ refers to the number of columns and $K_2$ points to the number of rows. It is worthwhile to mention that the  complexity order of TWSP is bottle-necked by computational burden for two pseudo-inverses and two singular vectors computation. Thus, the complexity can be expressed as  $\small{O(NK_1^2\!+\!MK_2^2\!+\!MN)}$ per iteration and the algorithm needs $\small{O(\text{max}(K_1,K_2))}$ iterations. Moreover, operator rnd($K$) refers to an integer random number less than or equal to $K$. In the next section, we evaluate the performance of our proposed algorithm.

The CSSP and CUR decomposition problems are NP-hard, i.e., a combinatorial search is required to find the best columns and rows. The proposed TWSP algorithm minimizes the main cost function \eqref{eq:main} in practice. However, there is no theoretical guarantee for convergence of the proposed TWSP algorithm Alg. \ref{alg:jcrss}.
In order to improve the convergence behavior of TWSP, we evade updating both columns and rows in each iteration. Rather, we prioritize updating  a row or a column depending on which exhibits a smaller projection error, i.e., which is a better minimizer for the cost function per each iteration. The implementation steps of TWSP are summarized in Alg. \ref{alg:jcrss}. As seen in Alg. \ref{alg:jcrss}, the method adds up a revision feature to remedy the greediness associated with consecutive selection by removing samples and optimizing to select a sample for replacement and reiterating this procedure.

\begin{algorithm}[h]
\caption{Two way spectrum pursuit (TWSP)}
\label{alg:jcrss}
\begin{algorithmic}
\small{
\REQUIRE {\small{${\boldsymbol{X}}\!\in\! \mathbb{R}^{N\times M}\!$, $K_1$ and $K_2$.}}\\ \vspace{1.2mm}
\hspace{-3.3mm}\textbf{Output:} $\mathbb{S}^c$ and $\mathbb{S}^r$}.
\STATE \textbf{Initialization:} \\
 $\mathbb{S}^c\leftarrow$A random subset of $\{1,\ldots,M\}$ with $|\mathbb{S}^c|=K_1$ \\
  $\mathbb{S}^r\leftarrow$A random subset of $\{1,\ldots,N\}$ with $|\mathbb{S}^r|=K_2$ \\
 $\{\mathbb{S}^c_k\}_{k=1}^K\!\!\leftarrow\!$ Partition $\mathbb{S}^c$ into $K_1$ subsets.\\
  $\{\mathbb{S}^r_k\}_{k=1}^K\!\!\leftarrow\!$ Partition $\mathbb{S}^r$ into $K_2$ subsets\\
 $i=\text{rnd}(K_1)$ and $j=\text{rnd}(K_2)$ and $\boldsymbol{X}=\boldsymbol{CUR}$ (CUR of $\boldsymbol{X}$)\\
 while a stopping criterion is not met \\%increase iter by $1$
 \vspace{+2mm}
  $\qquad \mathbb{S}^c_{\overline{i}}=\mathbb{S}^c\backslash \mathbb{S}^c_i$\\
  \vspace{+1mm}
  $\qquad \boldsymbol{C}_i\leftarrow\;\;\text{remove column}\;i\; \text{in matrix}\;\boldsymbol{C}$\\
  \vspace{+1mm}
  $\qquad \boldsymbol{W}= \boldsymbol{C}_i^{\dagger}\boldsymbol{X}\boldsymbol{R}^{\dagger}$\\
  $\qquad \boldsymbol{E}^c=\boldsymbol{X}-\boldsymbol{C}_i\boldsymbol{WR}\;\;$ {\footnotesize (Null space projection)}\\
  \STATE $\qquad$ $\boldsymbol{c}=\;$find the first left singular vector of $\boldsymbol{E}^c$ \eqref{eq:problem_sep1}\\
\STATE $\qquad$ $\mathbb{S}_i^c\xleftarrow[]{}$  the most correlated column of $\boldsymbol{E}$ with $\boldsymbol{c}$  \eqref{eq:problem_match1}
  \vspace{1mm}
\STATE $\qquad$ $\tilde{\mathbb{S}}^c\xleftarrow[]{} \bigcup_{i'=1}^{K_1} \mathbb{S}^c_{i'}$
  \vspace{1mm}
\STATE $\qquad \boldsymbol{C}=\boldsymbol{X}(:,\tilde{\mathbb{S}}^c)$
  \vspace{1mm}
\STATE $\qquad e^c=\underset{U}{min}\|\boldsymbol{X}-\boldsymbol{C}\boldsymbol{U}\boldsymbol{R}\|_F$\\
  \vspace{3mm}
  $\qquad \mathbb{S}^r_{\overline{j}}=\mathbb{S}^r\backslash \mathbb{S}^r_j$\\
  \vspace{+1mm}
        $\qquad \boldsymbol{R}_j\leftarrow\;\;\text{remove row}\;j\; \text{in matrix}\;\boldsymbol{R}$\\
    $\qquad \boldsymbol{Y}= \boldsymbol{C}^{\dagger}\boldsymbol{X}\boldsymbol{R}_j^{\dagger}$\\
    \vspace{+1mm}
  $\qquad \boldsymbol{E}^r=\boldsymbol{X}-\boldsymbol{CYR}_j \;\;${\footnotesize (Null space projection)}\\
  \STATE $\qquad$ $\boldsymbol{r}=\;$find the first right singular vector of $\boldsymbol{E}^r$ \eqref{eq:problem_sep2} \\
\STATE $\qquad$ $\mathbb{S}_j^r\xleftarrow[]{}$ the most correlated row of $\boldsymbol{E}$ with $\boldsymbol{r}$ \eqref{eq:problem_match2}
\STATE $\qquad$ $\tilde{\mathbb{S}}^r\xleftarrow[]{} \bigcup_{j'=1}^{K_2} \mathbb{S}^r_{j'}$ \vspace{+1mm}
\STATE $\qquad \boldsymbol{R}=\boldsymbol{X}(\tilde{\mathbb{S}}^r,:)$
  \vspace{1mm}
\STATE $\qquad e^r=\underset{U}{min}\|\boldsymbol{X}-\boldsymbol{C}\boldsymbol{U}\boldsymbol{R}\|_F$\\
\vspace{+2mm}
\STATE $\qquad \text{IF }e^c<e^r$
\STATE $\qquad \qquad \mathbb{S}^c\xleftarrow[]{}\tilde{\mathbb{S}}^c$\\
\STATE $\qquad \qquad i=\text{rnd}(K_1)$\\
\STATE $\qquad \text{else}$
\STATE $\qquad \qquad \mathbb{S}^r\xleftarrow[]{}\tilde{\mathbb{S}}^r$\\
\STATE $\qquad \qquad j=\text{rnd}(K_2)$\\
\end{algorithmic}
\end{algorithm}
\vspace{-2mm}
\section{Experimental Results}
\label{sec:exp}
In order to evaluate TWSP on machine-learning tasks in terms of CUR decomposition accuracy, we apply the proposed TWSP on synthetic data as well as three real applications.
\vspace{-2mm}
\subsection{CUR Decomposition on Synthetic Data}
In order to evaluate the general performance of TWSP, we compared it with the state-of-the-art methods for selecting columns and rows. In this regards, we created a $1000\times2000$ synthetic dataset. The dataset is generated by a rank-$30$ matrix contaminated with random noise. In Fig.~2, we have illustrated the CUR decomposition error for selecting a subset of rows and columns in the range of $2$ to $20$.  The reconstruction error of CUR is normalized by $\|\boldsymbol{X}\|_F^2$. We employ SP  as the state-of-the-art algorithm for column subset selection \cite{joneidi2020select}. We perform SP on the data matrix $\boldsymbol{X}$ and its transpose in order to, respectively, select a subset of columns and rows independently. Then, employing \eqref{eq:core} results in a CUR decomposition. We refer to the algorithm in \cite{mahoney2009cur} as adaptive CUR.
A more accurate algorithm for solving CUR decomposition results in a bigger blue region in Fig. 2. TWSP exhibits the best performance in this experiment. The convergence behavior of TWSP for this experiment is shown in Fig. 3  for selecting 20 columns and 20 rows. The final solution of the TWSP algorithm depends on the initial selected columns and selected rows. However, regardless of the initial condition, the TWSP algorithm minimizes the cost function of CUR decomposition. TWSP prioritizes in selection of columns or rows such that the largest decrease in the projection error is obtained. It is the main feature to obtain convergence in practice.

\vspace{-2mm}
\subsection{Joint Sensor Selection and Channel Assignment}

\label{subsec:joint}
The output products of CUR decomposition are not limited to a subset of columns and  rows. In some applications, interestingly, matrix $\boldsymbol{U}$ is the most important output of a CUR decomposition. Entry ($i,j$) in $\boldsymbol{U}$ indicates how important the contribution of the $i^{\text{th}}$ column and the $j^{\text{th}}$ row is to reconstruct the whole matrix $\boldsymbol{X}$. This interesting property is utilized for the problem of joint sensor selection and channel assignment in a cognitive radio network. To this aim the exact setup in \cite{zhang2020spectrum} is considered with $900$ grid points and $32$ frequency channels. The received power magnitudes are organized in a $900\times 32$ matrix.

The only difference here, is that the uniform sampling pattern of sensing  is replaced by the selection based on the CUR decomposition. Our proposed TWSP algorithm provides a fast and accurate solution for CUR decomposition. We select between $20$ and $80$ locations for spectrum sensing and all  $32$ channels. Each row of matrix $\boldsymbol{U}$  corresponds to a selected location and it has $32$ entries. We are to assign $F$ channels for each selected sensor. In other words, each location does not sense the whole spectrum and only $F$ frequency channels are assigned to each sensor.  The top-$F$  entries in each row with the highest absolute value show the most important channels for the corresponding locations to be sensed.  \\
Fig. \ref{fig:joint_sense} shows the cartography error of spectrum sensing for the conventional random selection as introduced in \cite{zhang2020spectrum} and our proposed optimized joint sensors and channels. For each sampled locations $F=8$ channels out of $32$ channels are sensed. The sampled spectrum map is interpolated using thin plane splines method \cite{ureten2012comparison} for both sampling methods. In addition to visual superiority in reconstruction of the spectrum map, channel assignment based on TWSP provides a better quantitative error as observed in Fig. \ref{fig:joint_sense}

\begin{figure}[t]
\vspace{-2mm}
\begin{subfigure}{0.48\textwidth}
\centering
\includegraphics[width=3.3 in,height=0.9in]{Figures/joint_chnl_snsr2.eps}
\end{subfigure}
\begin{subfigure}{0.48\textwidth}
\centering
\includegraphics[width=2.5 in]{Figures/joint_sns_ch_err.eps}
\end{subfigure}
\vspace{-1mm}
\label{fig:joint_channel}
\vspace{-3mm}
\caption{\footnotesize{The original spectrum map and its comparison with the interpolated map using random sampling and our proposed method. The interpolation error is depicted versus number of sensed locations.}\vspace{-3mm}}
\label{fig:joint_sense}
\end{figure}

\subsection{Supervised Sampling}
The proposed TWSP algorithm is an unsupervised data selection algorithm. In supervised settings, labels can infuse information to perform a more viable joint selection \cite{esmaeili2018transduction}. A naive approach is to select representatives from each class independently. However, considering classes jointly is more effective for data reduction. Assume we are given two data classes as $\boldsymbol{X}_1\in \mathbb{R}^{N\times M_1}$ and $\boldsymbol{X}_2\in \mathbb{R}^{N\times M_2}$. The goal is to select $K_1$ samples from class 1 and $K_2$ samples from class 2. To this aim, we propose to construct the cross correlation of two classes as a kernel representation for both classes jointly.  Matrix $\boldsymbol{X}={\boldsymbol{X}_2^T\boldsymbol{X}_1}$ which has $M_1$ columns and $M_2$ rows is fed to TWSP algorithm in order to select $K_1$ columns and $K_2$ rows jointly.

As an initial experiment, supervised sampling is performed on Kaggle cats and dogs dataset. The features are obtained by a trained Resnet-18 deep learning model as explained in \cite{he2016deep}. Three mutually exclusive data subsets for training, validation, and testing are partitioned randomly from $2000$ images of each class. The classification accuracy of $97.5$\% is achieved from a fine-tuned Resnet-18 using the whole training set containing $1000$ samples for each class. Afterwards, samples are selected by applying TWSP on the kernel feature matrix and the Resnet-18 is fine-tuned by using the sampled data. The model's accuracy is compared on the testing set with other sampling methods. Fig. \ref{fig:supervised} shows the performance of selection algorithms for different numbers of representatives per class. Using only two samples from each class, a classification accuracy of $82.3\%$ can be achieved which is more than $15\%$ improvement compared to random selection and more than $5\%$ improvement compared to other competitors.

\begin{figure}[!b]
\vspace{-4mm}
\centering
\includegraphics[width=2.8in, height=1.55in]{Figures/super_vised1.eps}
  \caption{\small{The classification accuracy of a fine-tuned Resnet-18 network using a few selected data per each class.}}
\label{fig:supervised}
\end{figure}

We have conducted further experiments to study the effectiveness of the proposed algorithm in the multi-class image classification problem. For this study, we use the Resnet-34 deep learning model pre-trained on CIFAR10 and trained on a subset of ImageNet Dataset comprising of $10$ classes. The classes used for this experiment are Tench, Goldfish, Great white shark, Tiger shark, Hammerhead, Electric ray, Stingray, Cock, Hen and Ostrich. The original training data consists of $1300$ images of each class, and the idea is to use TWSP to select the data samples such that $K$ samples of each class are used for training.
After training Resnet-34 for 10 classes of CIFAR10, the feature vectors of all the training images are extracted such that a 1300X512 matrix is obtained for each class. Since, TWPS is applicable for the two-class problem, we employ the one-versus-all approach. In other words, a separate cross correlation matrix is generated for each class such that $\boldsymbol{X}_1$ has the dimensions 512X1300 and $\boldsymbol{X}_2$ has the dimensions 512X11700 where $\boldsymbol{X}_1$ represents the feature vector of the class for which samples will be selected while $\boldsymbol{X}_2$ represents the feature vectors of the remaining $9$ classes. Hence, the number of cross correlation matrices is equal to the number of the classes. Then, $K_1$ rows and $K_2$ columns are selected separately by applying TWSP on each cross correlation matrix. This step of generating correlation matrices is different from binary classification where only one correlation matrix was formed. Pre-trained Resnet-34 on CIFAR10 has been refined again by using the sampled data from Imagenet dataset. The model’s accuracy is subsequently compared on the test set with other sampling methods in Fig. \ref{fig:multi_supervised}. As it can be seen, Our proposed TWSP algorithm shows a superior performance comparison to the state of the art methods in sampling.
\begin{figure}[!t]
\centering
\includegraphics[width=2.9in, height=1.55in]{Figures/super_vised1.eps}
  \caption{\small{Test accuracy in terms of improvement compared with the test accuracy obtained by random selection.}}
\label{fig:multi_supervised}
\vspace{-3mm}
\end{figure}
\vspace{-1mm}
\subsection{Informative Users/Contents Detection}
Another problem gaining a lot of interest by streaming  service providers is choosing a set of users to decisively reflect their feedbacks about different products. Therefore, it is crucial for such companies to find a subset of users and media products, reviews of which contains information about other users' unknown behavior.
Each user has a limited scope of interest.
For example, a user who only loves romance and action genres corresponds to a specific personality that is able to represent a cluster of users accurately. Moreover, such reviews for their areas of interest are more valuable, not reviews for all genres.
Moreover, there exist movies well representing their genres.

\begin{figure}[!t]
\centering
\includegraphics[width=2.9in, height=1.55in]{Figures/result_netflix.eps}
  \caption{\small{Comparison of the normalized prediction error  with state-of-the-art algorithms obtained by CUR decomposition for simultaneous movies and users subset selection from Netflix dataset.}}
  \vspace{-4mm}
\label{fig:Netflix}
\end{figure}

As a result, there exists a demand for a reliable algorithm to simultaneously choose the most informative subset of users and movies. Such a subset is desirable for streaming companies to the extent that they are willing to give the users incentives to leave comprehensive reviews for certain products. In this regards, we have evaluated our algorithm on Netflix Prize dataset containing $17,770$ movies and $480,189$ users. We have reduced the dataset to $990$ movies and $4,727$ users by considering only movies and users with most reviews. Then, we completed the dataset by Lin et al. method to have a ground truth \cite{lin2010augmented}.
We select a subset of rows and columns (users and media) from the completed dataset. Fig. \ref{fig:Netflix} reveals that TWSP shows the best performance in terms of predicting scores for all users/movies based on a few selected users/movies.

\vspace{-1mm}
\section{Conclusion}
Two-way spectrum pursuit is proposed as an accurate and efficient algorithm for solving CUR decomposition. TWSP can be employed for joint selection of columns and rows such that their outer products is able to reconstruct the whole matrix as accurate as possible. Some applications of the proposed algorithm are presented to show the efficacy of the proposed method. However, they are not limited to the mentioned applications. Moreover, the proposed algorithm can be extended to $n$-way spectrum pursuit for efficient tensor subset selection.

\bibliographystyle{plain}