\section{Introduction}

3D lane detection, comprising of an accurate estimation of the 3D position of the drivable lanes relative to the host vehicle, is a crucial enabler for autonomous driving. Two complementary technological solutions exist: loading pre-mapped lanes generated off-line~\cite{Urmson_boss} and perception-based real-time lane detection~\cite{aharon}. The off-line solution is geometrically accurate given precise host localization (in map coordinates) but complex to deploy and maintain. The most common perception-based solution uses a monocular camera as the primary sensor for solving the task. Existing camera-based methods detect lanes in the image domain and then project them to the 3D world by assuming a flat ground~\cite{aharon}, leading to inaccuracy not only in the elevation but also in the lane \textit{curvature} when the assumption is violated. Inspired by recent success of convolutional neural networks (CNNs) in monocular depth estimation~\cite{Liu_Depth}, we propose instead to directly detect lanes in 3D. More formally, given a single image taken from a front-facing camera, the task is to output a set of 3D curves in camera coordinates, each describing either a lane delimiter or a lane centerline.

\begin{figure}[ht]
	\begin{center}
		\includegraphics[width=\linewidth]{figures/example_result2.pdf}
	\end{center}
	\caption{\textbf{End-to-end approach illustrated}. \textbf{Left:} Output represented in top view. \textbf{Top-right: }result visualized in 3D. \textbf{Bottom-right:} result projected to original input image}
	\label{fig:tusimple}
\end{figure}

3D-LaneNet, our proposed solution, is a deep CNN that performs 3D lane detection. The network, trained end-to-end, outputs in each longitudinal road slice, the confidence that a lane passes through the slice and its 3D curve in camera coordinates. Our approach is schematically illustrated in Figure~\ref{fig:tusimple}. Our direct, single-shot approach avoids post-processing used in existing methods such as clustering and outlier rejection. The network's backbone is based on a novel dual pathway architecture that uses several in-network projections of feature maps to virtual bird-eye-view. This dual representation gives the network an enhanced ability to infer 3D in a road scene, and may possibly be used for other tasks requiring this ability (e.g. 3D car detection). The output is represented by a new column-based anchor encoding which makes the network horizontally invariant and enables the end-to-end approach. Each output is associated to an anchor in analogy to single-shot, anchor-based object detection methods such as SSD~\cite{SSD} and YOLO~\cite{YOLO}. Effectively, our approach casts the problem as an object detection one, in which each lane entity is an object, and its 3D curve model is estimated just like bounding box for an object.

\begin{figure*}[th]
	\begin{center}
		\begin{tabular}{c|c|c}
			\includegraphics[trim={2.2cm 7.4cm 2.15cm 2.8cm},clip=true,width=0.3\linewidth]{figures/3D_Examples/0196.png} &
			\includegraphics[trim={2.2cm 7.4cm 2.15cm 2.8cm},clip=true,width=0.3\linewidth]{figures/3D_Examples/0157.png} &
			\includegraphics[trim={2.2cm 9.4cm 5.15cm 2.8cm},clip=true,width=0.3\linewidth]{figures/examples_real/0720.png} \\

			\includegraphics[trim={2.2cm 7.4cm 2.15cm 2.8cm},clip=true,width=0.3\linewidth]{figures/3D_Examples/0090.png} &
			\includegraphics[trim={2.2cm 7.4cm 2.15cm 2.8cm},clip=true,width=0.3\linewidth]{figures/3D_Examples/0143.png} &
			\includegraphics[trim={2.2cm 9.4cm 5.15cm 2.8cm},clip=true,width=0.3\linewidth]{figures/examples_real/0162.png} \\

			\includegraphics[trim={2.2cm 7.4cm 2.15cm 2.8cm},clip=true,width=0.3\linewidth]{figures/3D_Examples/0575.png} &
			\includegraphics[trim={2.2cm 7.4cm 2.15cm 	2.8cm},clip=true,width=0.3\linewidth]{figures/3D_Examples/0841.png} &
			\includegraphics[trim={2.2cm 9.4cm 5.15cm 2.8cm},clip=true,width=0.3\linewidth]{figures/examples_real/0030.png} \\

		\end{tabular}
	\end{center}
	\caption{\textbf{Result visualization on test images. Centerline detection on } \textit{synthetic-3D-lanes} \textbf{examples (Left and middle columns) and Delimiter detection on } \textit{3D-lanes} \textbf{real image examples (Right column)}.  Detections with confidence $>0.5$ are shown. Ground truth \emph{\color{blue}(blue)} and method result  \emph{\color{red}(red)} shown in each image alongside a 3D visualization. Note that 3D axes are scene adaptive. Ignored lanes are marked in \emph{\color{cyan}cyan}. The leftmost bottom example shows a failure in correctly assigning a lane split, probably caused by occlusion.}
	\label{fig:results}

\end{figure*}

We validate our approach on three different datasets. The primary dataset used to develop the approach is a new computer-graphics dataset, \textit{synthetic-3D-lanes}\footnote{https://sites.google.com/view/danlevi/3dlanes} providing full access to the exact 3D position of each lane element. While several driving simulators exist~\cite{CARLA, AirSim}, they are not focused on the 3D lane detection task, and are limited in the variability of the relevant scene properties (e.g. lane curvature). Our main accomplishment in this domain, is the ability to randomly generate road segments with highly variable 3D shapes and lane topology. We therefore used it as the primary dataset for evaluation and ablation study. To validate our approach on real-world images we collected an additional dataset, \textit{3D-lanes}, from a vehicle-mounted front camera. 3D Lane annotation was accomplished using a Lidar scanner in a semi-manual annotation scheme. Finally, to compare with state-of-the-art lane detection methods, which only operate in the image domain, we adapted our method to this task, and demonstrated end-to-end image-only lane detection. This image-only version is tested on the tuSimple dataset~\cite{tu_simple}, reaching results competitive to state-of-the-art without the common post-processing techniques. To summarize, our main contributions are:

\begin{itemize}\addtolength{\itemsep}{-0.5\baselineskip}
	\item Introduction of a novel problem: single-frame 3D lane detection without geometric assumptions, along with new evaluation metrics
	\item A novel dual-pathway architecture deploying intra-network feature map IPM projections
	\item An new anchor-based output representation of lanes enabling a direct, end-to-end trained network, for both 3D and image-based lane detection.
	\item A methodology for generating randomly synthetic examples with variation in lane topology (i.e. number of lanes, merges, splits) and 3D shape.
\end{itemize}

\section{Related Work}

Traditional lane detection systems (e.g.~\cite{Gopalan}) combine low-level operations such as directional filters with high-level heuristics such as the Hough transform to detect continuous lanes in the image. A common pipeline includes 4 stages: local lane-feature extraction (1), lane model fitting (2), image-to-world correspondence (3) and temporal aggregation (4). Bar-Hillel \textit{et al.}~\cite{aharon} provide this modular decomposition alongside a detailed overview of traditional systems. In recent years, the \textit{local feature extraction stage} is performed by applying one or more CNNs to the image, but the overall pipeline remains very similar and the latter post processing stages remain.

Initially, CNNs were used to improve feature extraction by either enhancing the edge map (Kim and Lee~\cite{kim_lee}) or classifying candidate patches (He \textit{et al.}~\cite{He}). Huval \textit{et al.}~\cite{huval} detects local lane line segments with an object detection CNN. VPGNet (Lee \textit{et al.}~\cite{lee}), follows a similar concept and additionally detects other road markings and the vanishing point to improve lane detection. Kim and Park~\cite{kim_park} re-formulate the local-feature extraction stage as a semantic-segmentation problem, with two classes corresponding to the left and right lane delimiters, extending the reach of the network to perform clustering. However, a world-coordinate lane model must still be fitted to each cluster, and multiple lanes are not handled. Neven \textit{et al.}~\cite{neven} make an attempt towards end-to-end multi-lane detection, by training a CNN not only to create a binary lane pixel mask but also a feature embedding used for clustering lane points. Ghafoorian \textit{et al.}~\cite{Ghafoorian} propose applying a generative adversarial network to make the semantic segmentation network output more realistic in the context of lane detection. Several works (e.g. Meyer \textit{et al.}~\cite{Meyer}, Oliveira \textit{et al.}~\cite{Oliveira}) are built on a similar approach in which the host and possibly adjacent lanes are the semantic classes (lane interior rather than the lane delimiters).

As opposed to all presented methods, 3D-LaneNet unifies the first three stages of the common pipeline by providing a full multi-lane representation in 3D world coordinates directly from the image in a single feed-forward pass. In addition, previous methods use the flat ground assumption for the image-to-world correspondence while our method fully estimates the parametrized 3D curves defining the lanes. Only a few methods directly address 3D lane estimation: ~\cite{Nedevschi2004}, using stereo, and ~\cite{Xiong2018, Coulombeau_Laurgeau2002} which follow a multi-view geometry approach and assume a known constant road / lane width to solve the depth ambiguity. Instead, we use a data driven approach and make no geometric assumptions.

Inverse perspective mapping (IPM) generates a virtual top-view (sometimes called bird-eye-view) of the scene from a camera-view as in the example in Figure~\ref{fig:tusimple}. It was introduced in the context of obstacle detection by Mallot \textit{et al.}~\cite{Mallot} and first used for lane detection by Pomerleau~\cite{Pomerleau}. IPM has since been extensively used for lane detection (e.g.~\cite{Borkar,Aly}) since lanes are ordinarily parallel in this view and their curvature can be accurately fitted with low-order polynomials. In addition, removing the perspective effects causes lane markings to look similar (except for blurring effects) regardless of their distance from the camera. Most recently He \textit{et al.}~\cite{He} introduced a ``Dual-view CNN'' which is composed of two separate sub-networks, each producing a descriptor (one per view) which are then concatenated and applied to candidate image locations. Li \textit{et al.}~\cite{Li} use a CNN to detect lane markings along with geometrical attributes, such as local position and orientation, directly on a top-view image which preserves invariance to these properties. In addition they deploy a second, recurrent network, that traverses the image to detect consistent lanes. Neven \textit{et al.}~\cite{neven} use the horizon, predicted in each image by a sub-network (``H-net''), to project the lanes to top-view for improved curve fitting. In contrast to previous work, we exploit both views in a synergistic single network approach.

More generally, we propose the first method that uses an end-to-end trained CNN to directly detect multiple lanes and estimate the 3D curvature for each such lane. We also show that our method is applicable both to centerlines and delimiters with an ability to handle splits and merges as well, without any further post-processing.

\section{Method}

Our method gets as input a single image taken from a front facing camera mounted on a vehicle as illustrated in Figure~\ref{fig:geometry}. We assume known intrinsic camera parameters $\kappa$ (e.g. focal length, center of projection). We also assume that the camera is installed at zero degrees roll \textit{relative to the local ground plane}. We do not assume a known camera height and pitch since these may change due to vehicle dynamics. The lanes in a road scene can be described both by the set of centerlines $\left\lbrace C_i\right\rbrace_{i=1}^{N_C}$ of each lane and by the set of lane delimiters $\left\lbrace D_i\right\rbrace_{i=1}^{N_D}$ as illustrated in Fig.~\ref{fig:out}. Each such lane entity (centerline or delimiter) is a curve in 3D expressed in camera coordinates ($\mathcal{C}_{camera}$).  The task is to detect either the set of lane centerlines and/or lane delimiters given the image.

\begin{figure}[t]
	\begin{center}
		\includegraphics[width=0.79\linewidth]{figures/geometry2.pdf}
	\end{center}
	\caption{\textbf{Camera position and road projection plane}}
	\label{fig:geometry}

\end{figure}

\begin{figure*}[t]
	\begin{center}
		\includegraphics[width=\linewidth]{figures/network_architecture2.pdf}
	\end{center}
	\caption{\textbf{3D-LaneNet network architecture.}}
	\label{fig:net}
\end{figure*}

\begin{figure}[t]
	\begin{center}
		\includegraphics[width=0.7\linewidth]{figures/top_view2.pdf}
	\end{center}
	\caption{\textbf{Output representation.} Note that the number of anchors ($N$) equals the output layer width (marked by $\hat{w}/8$ in Fig.~\ref{fig:net}). The geometry representation for one lane \textit{centerline} is shown. Lane \textit{delimiters} (white dashed curves) are similarly represented.}
	\label{fig:out}
\end{figure}

\subsection{Top-view projection}\label{sec:IPM}
We briefly review Inverse Perspective Mapping (IPM). In short, IPM is a homography that warps a front view image to a virtual top view image as depicted in the top-left image of Figure~\ref{fig:tusimple}. It is equivalent to applying a camera rotation homography (view is rotated downwards) followed by an anisotropic scaling~\cite{Hartley2004}. In our implementation we want to ensure that each pixel in the top view image corresponds to a predefined position on the road, \textit{independent} of the camera intrinsics and its pose relative to the road.

See Figure~\ref{fig:geometry} for an illustration of the following definitions. The camera coordinates $\mathcal{C}_{camera}=\left(\acute{x},\acute{y},\acute{z}\right)$  are set such that $\acute{y}$ is the camera viewing direction. Let $P_{road}$ be the plane \textit{tangent} to the local road surface. We define the road coordinates $\mathcal{C}_{road}=\left(x,y,z\right)$ as follows: $z$ direction is the normal to $P_{road}$, $y$ is the projection of $\acute{y}$ onto $P_{road}$ and the origin is the projection of the camera center onto $P_{road}$. Let $T_{c2r}$ be the 6-D.O.F. transformation between $\mathcal{C}_{camera}$ and  $\mathcal{C}_{road}$ (3D translation and 3D rotation). Since we assume zero camera roll, $T_{c2r}$ is uniquely defined by the camera pitch angle $\theta$ and its height above the ground $h_{cam}$. The homography $H_{r2i}:\mathbb{P}^2\mapsto\mathbb{P}^2$, mapping each point on $P_{road}$ to image plane coordinates, is determined by $T_{c2r}$ and $\kappa$ (See ~\cite{Hartley2004}, Section 8.1.1). Finally, the IPM is obtained from $H_{r2i}$ using a fixed set of parameters $IPM_{Params}$  defining the top view region boundaries and an anisotropic scaling from meters to pixels. The top view image is generated using bilinear interpolation defined by a sampling grid $S_{IPM}$.

\label{sec:dual_block}

\subsection{Network structure}

An overview of the 3D-LaneNet is illustrated in Figure~\ref{fig:net}. Information is processed in two parallel streams or pathways: the image-view pathway and the top-view pathway. We call this the \emph{dual-pathway backbone}. The image-view pathway processes and preserves information from the image while the top-view pathway provides the features with translation invariance and is used to predict the 3D lane detection output. The architecture for the image-view pathway is based on VGG16~\cite{VGG} while the top-view pathway is similarly structured. Information flows to the top-view pathway through four \textit{projective transformation} layers as follows.

\subsection{The projective transformation layer}

A main building block in our architecture is the \emph{projective transformation layer} marked in blue in Fig.~\ref{fig:net}. This layer is a specific realization, with slight variations, of the \textit{spatial transformer  module}~\cite{spatial_transform}. It performs a differentiable sampling of input feature maps corresponding spatially to the \textit{image plane}, to output feature maps corresponding spatially to a \textit{virtual top view} of the scene while preserving the \# of channels. The differential sampling is achieved through a grid generated as described in Sec. ~\ref{sec:IPM}, using a IPM predicted by the \textit{Road projection prediction branch} as described in the next section. The resulting projected feature maps, except for the first set, are concatenated to downstream feature maps from the top-view pathway. A subsequent neuron, operating on the concatenated feature maps combines the following two desirable properties for lane detection. First, translational invariance in the top-view plane. This is valuable since in the top view lanes have similar appearance and geometry across the space. Second, preservation of a dual information context - in both image and top view. The additional image-view context encodes information which is not present in the top view such as fences, skyline and trees which are crucial for deciphering the 3D structure of the scene. Particularly, in the far range, the image-view context is much richer in visual information, and represents a much larger actual area compared to the top view.

\subsubsection{Road projection prediction branch}
\label{horizon_estimation}
The first intermediate output of the image-view pathway network is an estimation of  the ``road projection plane'' $P_{road}$. Essentially, this branch predicts $T_{c2r}$, the camera ($\mathcal{C}_{camera}$) to road ($\mathcal{C}_{road}$) transformation. It is trained in a supervised manner. $T_{c2r}$ determines the top-view homography $H_{r2i}$ and sampling grid $S_{IPM}$ as explained in Section~\ref{sec:IPM}, and is therefore needed for the feed-forward step of the top-view pathway. At inference time, it is used also to translate the network output which is expressed in $\mathcal{C}_{road}$, back to $\mathcal{C}_{camera}$. As described in Section~\ref{sec:IPM}, $T_{c2r}$ is defined in our case by the camera height $h_{cam}$ and pitch $\theta$, and therefore these are the two outputs of this branch.

\subsubsection{Lane prediction head}

At the heart of our end-to-end approach lies the \textbf{anchor-based lane representation}. Inspired by object detection, we use anchors to define lane candidates and a refined geometric representation to describe the precise 3D lane shape for each anchor. The output coordinate system is the estimation of $\mathcal{C}_{road}$ determined by $h_{cam}$, $\theta$. Our anchors correspond to longitudinal lines in this coordinate system and the refined lane geometry to 3D points relative to the respective anchor. As illustrated in Figure~\ref{fig:out}, we define the anchors by equally spaced vertical (longitudinal) lines in $x$-positions $\left\lbrace X_A^i \right\rbrace _{i=1}^{N}$. Per anchor $X_A^i$, a $3D$ lane is represented by $2\cdot  K$ output neurons activation $\left(\mathbf{x}^i,\mathbf{z}^i\right)=\left\lbrace \left(x_j^i,z_j^i\right)\right\rbrace_{j=1}^K$, which together with a fixed vector of $K$ predefined $y$ positions ($\mathbf{y}=\left\lbrace y_j \right\rbrace_{j=1}^K$) define a set of $3D$ lane points. The values $x_j^i$ are horizontal offsets \emph{relative to the anchor position $X_A^i$}. Meaning, the output $\left(x_j^i,z_j^i\right)$ represents the point $\left(x_j^i+X_A^i,y_j,z_j^i\right)\in\mathbb{R}^3$, in $C_{road}$ coordinates. In addition, for each anchor $i$, we output the confidence $p^i$ that there is a lane associated with the anchor. We use a predefined longitudinal coordinate $Y_{ref}$ for the association. The anchor $X_A^i$ associated to a lane is the one closest to the x-coordinate of the lane at $y=Y_{ref}$.

Per anchor, the network outputs up to \emph{three} types ($t$) of lane descriptors (confidence and geometry), the first two ($c_1,c_2$) represent lane centerlines and the third type ($d$) a lane delimiter. Assigning two possible centerlines per anchor yields the network support for merges and splits which may often result in having the centerlines of two lanes coincide at $Y_{ref}$ and separating at different road positions as in the rightmost example in Figure~\ref{fig:out}. The topology of lane \textit{delimiters} is generally more complicated compared to centerlines and our representation cannot capture all situations (for example the lane delimiters not crossing $y=Y_{ref}$ in Fig.~\ref{fig:out}). The prediction head of the 3D-LaneNet is designed to produce the described output. Through a series of convolutions with no padding in the $y$ dimension, the feature maps are reduced, and finally the prediction layer size is $3\cdot \left(2\cdot K +1\right)\times 1 \times N$ s.t. each column $i\in\left\lbrace 1\ldots N\right\rbrace$ corresponds to a single anchor $X_A^i$. Per anchor, $X_A^i$ and type $t\in \left\lbrace c_1, c_2, d\right\rbrace$ the network output is denoted by $\left(\mathbf{x}_t^i,\mathbf{z}_t^i,p_t^i\right)$.
The final prediction performs a 1D non-maximal suppression as common in object detection: only lanes which are locally maximal in confidence (compared to the left and right neighbor anchors) are kept. Each remaining lane, represented by a small number ($K$) of 3D points, is translated to a smooth curve using spline interpolation.

\subsection{Training and ground truth association}\label{sec:training}

Given an image example and its corresponding 3D lane curves, $\left\lbrace C_i\right\rbrace_{i=1}^{N_C}$ (centerlines) and $\left\lbrace D_i\right\rbrace_{i=1}^{N_D}$ (delimiters), the training proceeds as follows. First, the ground truth (GT) coordinate system $\mathcal{C}_{road}$ is defined for the local road tangent plane as described in Sec.~\ref{sec:IPM} using the known pitch ($\hat{\theta}$) and camera height ($\hat{h}_{cam}$). Next, each lane curve, projected to the $x-y$ plane of $\mathcal{C}_{road}$, is associated with the closest anchor at $Y_{ref}$. The leftmost lane delimiter and leftmost centerline associated with an anchor are assigned to the $c_1$ and $d$ output types for that anchor. If an additional centerline is associated to the same anchor it is assigned to output type $c_2$. This assignment defines the GT per example in the same format as the output: per anchor $X_A^i$ and type $t$ the associated GT is denoted by $\left(\mathbf{\hat{x}}_t^i,\mathbf{\hat{z}}_t^i,\hat{p}_t^i\right)$, where $\hat{p}_t^i$ is an anchor/type assignment indicator, and the coords in $\mathcal{C}_{road}$.

Both in training time and in the evaluation, entire lanes are ignored if they do not cross $Y_{ref}$ inside valid top-view image boundaries, and lane points are ignored if occluded by the terrain (i.e. beyond a hill top). The overall loss function of the network is given in Eq.~\ref{eq:loss}. It combines three equally weighed loss terms: lane detection (Cross-entropy-loss), lane geometry and road plane estimation ($L_1$-loss).

\begin{align}
	\label{eq:loss}
	\begin{split}
		\mathcal{L}=&-\sum_{t\in\left\lbrace c_1,c_2, d\right\rbrace}\sum_{i=1}^N \left(\hat{p}_t^i\log p_t^i + \left(1-\hat{p}_t^i\right)\log\left(1-p_t^i\right)\right) \\
		& + \sum_{t\in\left\lbrace c_1,c_2, d\right\rbrace}\sum_{i=1}^N  \hat{p}_t^i\cdot\left(\left\|\mathbf{x}_t^i-\mathbf{\hat{x}}_t^i \right\|_1   +\left\|\mathbf{z}_t^i-\mathbf{\hat{z}}_t^i \right\|_1\right)  \\
		& + \left| \theta -\hat{\theta}\right| + \left| h_{cam} -\hat{h}_{cam}\right|
	\end{split}
\end{align}

\section{Experiments}
Our experimental work is presented as follows. We first present the methodology used for generating a new synthetic dataset \emph{synthetic-3D-lanes}, which is used to derive most of this study conclusions. Next, we introduce the \emph{3D-lanes} dataset generated for validation on real-world imagery. Using a newly proposed evaluation method for 3D lane detection, we then present results on both datasets, including an ablation study carefully examining the contribution of each concept in our overall approach. Finally, we compare an image-only version of 3D-LaneNet to existing state-of-the-art methods on the tuSimple benchmark~\cite{tu_simple}.

\begin{figure}[h]
	\begin{center}
		\begin{tabular}{c}
			\begin{tabular}{cc}
			(a)\includegraphics[width=0.47\linewidth]{figures/surface_486.png} & (b)\includegraphics[width=0.35\linewidth]{figures/topology_486.png}
			\end{tabular}\\
			(c)\includegraphics[width=0.7\linewidth]{figures/surface_with_road_486.png}
		\end{tabular}
	\end{center}
	\caption{\textbf{Synthetic scene generation example.} (a) Surface  (b) Road topology and curvature (c) Road on surface}
	\label{fig:synthetic}
\end{figure}

\subsection{Synthetic 3D lane dataset}\label{sec:synthetic}
We generated the \emph{synthetic-3D-lanes} dataset using the open source graphics engine \textit{blender}~\cite{blender}. Our programmatic approach allows us to randomize each of the modeled elements, from the 3D geometry of the scene to object types as illustrated in Figure~\ref{fig:synthetic}. The process of generating \textit{each scene} is composed of the following steps:

\textbf{Terrain 3D.} The terrain is modeled by a Mixture of Gaussians distribution with the number of Gaussians and their parameters randomized. Figure~\ref{fig:synthetic}(a) shows an example of such a terrain.

\textbf{Lane topology.} Number of lanes on the main road is selected. Then we choose if there is a secondary road and the number of lanes in it. Depending on the later direction of the camera in the scene the junction of the secondary road is viewed as a merge or a split.

\textbf{Lane top view geometry.} The geometry of the main road in top view is modeled by a 4\textsuperscript{th} degree polynomial producing mild to extreme curvatures. The junction point for merges / splits as well as the lane width are chosen. This results in a top view lane-level map as shown in Fig.~\ref{fig:synthetic}(b).

\textbf{Lane 3D.} The top-view lane map is placed on the terrain, and secondary roads are lifted to simulate common road topography. Fig.~\ref{fig:synthetic}(c) shows the result of this stage.

\textbf{Terrain and road appearance.} The texture of the road and the terrain are chosen from a set of textures. The type and color of the lane markings is also randomized.

\textbf{Objects.} Cars and trees selected from a set of models are placed in the scene, on and off the road respectively.

\textbf{Scene rendering.} The host vehicle camera is positioned on the main road by choosing its lane and a lateral offset around lane center. The camera height is set randomly between 140cm and 190cm and a downward pitch between 0 and 5 degrees is selected. Finally, illumination is set and the scene is rendered from the camera view. The 3D points of each lane centerline and delimiter are translated to camera coordinates to generate the ground truth.

Each generated example consists of an image ($360\times480$ pixels) and its associated ground truth: 3D lanes, camera height and pitch. Figure~\ref{fig:results}(Left and middle columns) presents several examples showing the resulting diversity and complexity. The exact parameters used in the random generation process are listed in Appendix I. The generated dataset contains 300K train and 5K test examples. An additional 1K validation set was used for learning rate scheduling and choosing the best performing snapshot.

\subsection{Real-world 3D lane dataset}

Acquiring ground truth labeled data with 3D for the task is an endeavor that requires a complex multiple-sensor setup and possibly also expensive HD maps. To this end we introduce a new such dataset, \textit{3D-lanes}, created using a multi-sensor setup including a forward-looking camera, a Velodine HDL32 lidar scanner and a high-precision IMU all synchronized and accurately aligned. The data was collected in 6 drives each on a different road segment totaling nearly 2 hours of driving. Using the Lidar and IMU we generate aggregated lidar top view images as in~\cite{Urmson_boss}, which are then used together with a semi-manual annotation tool for generating ground truth. In total, $85K$ images were annotated, out of which $1K$, consisting of a separate drive, were used as the test set and the remaining as the train set. The lidar information is additionally used to provide the full 3D curve of each lane. A disadvantage of this approach is that lanes not sufficiently visible to the lidar, due to occlusions or limited resolution at distance, are missing from the ground truth. Therefore, the labeling is somewhat noisy as can be observed in Fig.~\ref{fig:results}(Right column). In addition, the dataset variability in terms of geometry and topology is modest compared to the \emph{synthetic-3D-lanes} dataset. We therefore used the synthetic data which has perfect ground truth to develop the method and conduct ablation studies while the real-world dataset is used for validating transferability of the approach to real data and qualitative analysis.

\subsubsection{Evaluation results}

\textbf{Evaluation metrics.} We propose an evaluation of 3D lane detection that separates the detection accuracy from the geometric estimation accuracy. \textit{Detection accuracy} is computed via the standard average precision (AP) measure of the precision-recall curve. We first compute a curve to curve distance between a GT and a detected lane as a weighted sum of point-wise Euclidean distances. We measure distances on a set of predefined y-values along the curves, every 80cm in the range 0-80 meters. Weight is decreased for farther away points. We then perform a one-to-one (curve) matching by choosing pairs in decreasing similarity. A matching is considered correct if the weighted distance is below a certain, rather permissive, threshold (1.5 meters). Iterating over lane confidence thresholds we generate the precision-recall curve.

For matched detections we assess the \textit{geometric estimation accuracy} by measuring the distribution of the error (point-wise Euclidean distance) over the same points used to measure the curve to curve distance. We further divide the entire dataset to lane points in the \textbf{near range} (0-30m) and \textbf{far range} (30-80m) due to the differences in the magnitude of errors. We then compute per range the $1\sigma$ error, as the $68$ error percentile and the $2\sigma$ error as the $95$ percentile. Lane centerline and delimiter detection are separately evaluated using this methodology. Irrelevant lane points are ignored in the evaluation as in the training phase.

\textbf{Implementation details.} 3D-LaneNet and all the variants brought in this section were initialized and trained using an identical protocol. The image-view pathway is initialized from VGG16 trained on imagenet~\cite{imagenet_cvpr09}. We train with Adam optimization~\cite{adam} and with initial learning rate $5\cdot 10^{-4}$. We use a variation on the cyclic learning rate regime described in~\cite{cyclic_lr} with a minimal learning rate of $10^{-6}$. The y-range of the top view representation is 80 meters and the x-range is 20 meters. The IPM scale is different in x and y: in the first top-view feature map each pixel corresponds to 16cm laterally (x) and 38.4cm longitudinally (y). The last top-view feature map is $\times 8$ smaller and since there is an anchor per column the distance between anchors is $16\times8=128cm$. We set the $K(=6)$ vertical reference points to be $\mathbf{y}=\left\lbrace5,20,40,60,80,100 \right\rbrace$ and $Y_{ref}=20m$.

\textbf{Results on \textit{synthetic-3D-lanes} dataset.} Typical network results on the test set are shown in Figure~\ref{fig:results}, with ground truth marked. The first row in Table \ref{tab:ablation} shows the quantitative results of 3D-LaneNet for centerline detection. A valid concern with synthetic datasets is that its variability is too limited and that the learned network memorizes the entire example space instead of learning to generalize. A positive indication that this is not the case is that test AP (0.952) is well below train AP (0.966) as are the geometric evaluation measures. All networks trained in the ablation tests presented here were initialized from VGG16 just as the 3D-LaneNet was and were trained with the same training parameters and number of epochs.

We first examine the role of the \textit{\textbf{dual-pathway architecture}} by comparing it to alternative architectures. The \textbf{image-view} only version connects the image-view pathway directly to the lane detection head which outputs the representation in $\mathcal{C}_{road}$ exactly as 3D-LaneNet does. The anchor positions $X_A$ in this case are determined by the columns in the last feature map: for each column we pick a pixel at a predefined image y-coordinate and project it to top-view to determine the anchor corresponding to the column. The \textbf{top-view} only version first projects the \textit{image} itself to top view and continues the same computation as the top-view pathway. In addition, we tested two versions which include a limited version of the dual-pathway. The \textbf{early IPM} includes a single dual context module (the first amongst the four in the full network). The \textbf{late IPM} similarly contains only the last dual context module out of the four. The results,  summarized in table~\ref{tab:ablation}, show that the full dual-pathway architecture has superior performance compared to all other variants. In particular, the worst result is delivered by the \textbf{image-view} only version, stressing the importance of the top-view processing pathway. Note that the \textbf{late stage IPM}, consisting of a trimmed version of the dual pathway, delivers the second best accuracy, but with a reduced computational cost, making it a good candidate for real-time implementations.

We also tried alternative \textit{\textbf{definitions of the road projection plane}}. One approach takes into consideration the entire scene when fitting the road plane and not only the local road normal. To test it we devised a ground truth generation algorithm which takes the farthest visible road point and connects it to the local road position to determine the pitch. This method, is termed \textbf{horizon} in Table~\ref{tab:ablation} since it resembles horizon estimation methods. Evidently, it performed slightly worse in general, although we observed consistently cases in which the scene topography favors this definition. We also tried assuming a \textbf{fixed position} of the camera, in which the average pitch ($2.5^\circ$) and camera height (165cm) were used to define $T_{c2r}$. Finally, we note that learning to predict the best road projection plane per scene \textit{without explicit supervision}, as proposed in~\cite{neven}, failed to produce satisfying results for our task.

The last row in Table \ref{tab:ablation} (\textbf{flat ground}) is brought to stress the importance of full 3D lane estimation compared to the current existing approach: image-only detection and \textit{image-to-world} translation using the flat-ground assumption. Image-only detection is obtained by projecting the 3D-LaneNet results to the image plane. For the image-to-world stage we need to choose the plane to project the image result to. We tried two options, both computed using the ground truth: the road plane $P_{road}$ and the plane defined by \textbf{horizon} as described in the previous experiment. As one may expect, the \textbf{horizon} based method, which essentially uses the best planar fit for the entire scene, produced the better results, which are still inferior to those of 3D-LaneNet which performs full 3D estimation.

The \textit{delimiter detection} performance obtained by 3D-LaneNet is $0.971$ AP (positional errors: 12.9cm@$1\sigma$, 33cm@$2\sigma$ near range;  30cm@$1\sigma$, 106cm@$2\sigma$ far range). These metrics show a slightly better performance compared to centerline detection. A possible explanation is that delimiters are clearly marked on the road while centerlines are indirectly inferred. Since output is transformed from road to camera coordinates using an estimated $T_{c2r}$ we also measured the quality of this estimation and its effect on the results. The median values of the absolute errors for pitch ($\theta$) and camera height ($h_{cam}$) are $0.09^\circ$ and $2.4cm$ respectively. To eliminate the contribution of this error we evaluated performance in road coordinates $C_{road}$ by taking the raw network output (before transforming to $C_{camera}$) and got a negligible difference in measured performance.

\textbf{Results on \textit{3D-lanes} dataset.} For operation on real-world data we trained the 3D-LaneNet on the train part of the \textit{3D-lanes} dataset. Result examples from the respective test set are shown in Fig.~\ref{fig:results} (Right column). Note that since the camera is mounted with a downward pitch the 3D lanes are detected as rising upward. Evaluation metrics are presented in Table~\ref{tab:real}. As in the synthetic data, using the \textbf{flat ground} assumption on the real data degrades performance, achieving a 4 times larger error in the far range.

\setlength{\heavyrulewidth}{1.5pt}
\setlength{\abovetopsep}{4pt}
\begin{table}[!htbp]
	\centering
	\caption{Centerline detection results on \emph{synthetic-3D-lanes} dataset}
	\begin{tabular}
		{*6c}
		\toprule
		& \textbf{AP} &   \multicolumn{2}{c}{Error \textbf{near} (cm)} & \multicolumn{2}{c}{Error \textbf{far} (cm)}   \\
		&    & $1\sigma$ & $2\sigma$ & $1\sigma$ & $2\sigma$\\
		\textbf{3D-LaneNet} & \textbf{0.952}	& \textbf{13.3}& 	\textbf{34.4}& 	\textbf{33.1}& 	\textbf{122}\\
		\textbf{image-view} & 0.819	& 20.3 & 50	& 74.7	& 241  \\
		\textbf{top-view} & 0.929 &	17.5 &	39.6 & 49.5 & 208  \\
		\textbf{early IPM } & 0.934 & 13.7	 & 35.5 & 43.5 & 189 \\
		\textbf{late IPM } & 0.948	&14.5	&37.2	& 37.4	&139  \\
		\midrule

		\textbf{horizon} &0.949	&14.8&	40.4&	36.7&	132	  \\
		\textbf{fixed position }	& 0.948	& 13.6	& 37.3	& 35.4	& 139  \\
		\midrule
		\textbf{flat ground}&0.566	& 46.9 & 114	& 99	& 289  \\
		\bottomrule
	\end{tabular} 	\label{tab:ablation}
\end{table}

\setlength{\heavyrulewidth}{1.5pt}
\setlength{\abovetopsep}{4pt}
\begin{table}[!htbp]
	\centering
	\caption{Delimiter detection results on \emph{3D-lanes} dataset}
	\begin{tabular}
		{*6c}
		\toprule
		& \textbf{AP} &   \multicolumn{2}{c}{Error \textbf{near} (cm)} & \multicolumn{2}{c}{Error \textbf{far} (cm)}   \\
		&    & $1\sigma$ & $2\sigma$ & $1\sigma$ & $2\sigma$\\
		\textbf{3D-LaneNet} & \textbf{0.918}	& \textbf{7.5}& 	\textbf{19.6}& 	\textbf{12.4}& 	\textbf{33}\\
		\midrule
		\textbf{flat ground} & 0.894 & 19.1 &	37.4 & 64.1	& 137  \\
		\bottomrule
	\end{tabular} 	\label{tab:real}
\end{table}

\subsection{Evaluation of image-only lane detection}

The purpose of this experiment is to compare our approach to the current state of the art, which exists for image-only lane detection. The tuSimple lane dataset~\cite{tu_simple} consists of 3626 training and 2782 test images. Unfortunately, today there is no access to the labels for the test images. We therefore divide the original training set to our own train/validation sets (90\% train and 10\% validation). While we are aware that there may be deviations between our evaluation (obtained on the validation set) and the one on the test set, we can expect a similar performance and reach the same conclusion qualitatively. Since this dataset does not contain 3D information, we train a variation of 3D-LaneNet, which detects the lanes in the image domain. Instead of a 3D representation, the network output was reduced to 2D points on the road projection plane by eliminating the elevation ($\mathbf{z}_t^i$) component. Only the delimiter output type is maintained ($t=d$) since the marked entities in the dataset are lane delimiters. A fixed homography, $H_{tuSimple}$, between image plane and road projection plane, was manually selected, such that straight lanes become parallel in top view. The lanes directly predicted by the network are transformed to lanes in the image view using $H_{tuSimple}$. Since $H_{tuSimple}$ is fixed, the road projection plane prediction branch is not used. Other than the aforementioned, the network is identical to the 3D-LaneNet as configured for the synthetic-3D-lanes dataset. The tuSimple main evaluation metric (\textit{acc})~\cite{tu_simple} is the average ratio of detected ground truth points per image. Using our end-to-end approach on our validation set we reached an accuracy of $0.951$ which is competitive with the one achieved by the tuSimple 2017 competition winning method~\cite{SCNN}, ($0.965$). This result is encouraging and somewhat surprising given that our entire approach was designed towards the 3D estimation task. In particular, our geometric loss (Eq.~\ref{eq:loss}) is computed in top view coordinates, giving in practice a much higher weight to distant lane points while in the tuSimple $acc$ metric all points equally contribute.

\section{Conclusions}
We presented a novel problem, 3D multiple lane detection, along with an end-to-end learning-based solution, 3D-LaneNet. The approach has been developed using a newly introduced synthetic dataset and validated on real data as well. The approach is applicable in principle to all driving scenarios except for complex urban intersections. Finally, we believe that the \textit{dual-pathway} architecture can facilitate additional on-road important 3D estimation tasks such as 3D vehicle detection.
{
	\small
	\bibliographystyle{ieee_fullname}


\section{Introduction}

The widespread use of machine learning algorithms and fully autonomous systems has greatly transformed the industrial landscape of the twenty-first century. However, with these great advances comes a responsibility for researchers, developers, and regulators to consider the impact of these systems on the broader society. In 2014, the US presidential administration published a report on big data collection and analysis, finding that ``big data technologies can cause societal harms beyond damages to privacy" \citep{united2014big}. The report raised the concern that algorithmic decisions inferred from big data may have harmful biases,  potentially leading to discrimination against disadvantaged groups.

This drive towards ethical practices in machine learning has led to many developments in \emph{algorithmic fairness}. One such advance has been towards developing algorithms which display group fairness, also referred to as statistical, conditional or demographical parity. From a regulatory viewpoint, group fairness is particularly interesting as affirmative action policies have already been passed to address discrimination against caste, race and gender \cite{weisskopf2004affirmative, dumont1980homo, deshpande2017affirmative}. However, it is worth noting there may be a considerable cost involved in achieving such fairness in some cases \cite{corbett2017algorithmic}.

In a machine learning context, there have primarily been two approaches towards developing systems which demonstrate group fairness; data alteration endeavors to modify the original dataset in order to prevent discrimination between groups \cite{luong2011k, kamiran2009classifying} in contrast to regularisation which penalizes models for unfair behavior \cite{kamishima2011fairness, berk2017convex, calders2013controlling, calders2010three, raff2017fair}.

More recently, there has been an effort towards constraining models such that they prohibit unfair behavior, a stricter assertion than regularisation. This work directly follows \cite{aistats_submission} in which group fairness in expectation for regression models is investigated, defined as:

\textbf{Group Fairness in Expectation (GFE):} \emph{A regressor $f(\cdot): X \rightarrow Y$ achieves group fairness in expectation with respect to groups $A, B \subseteq X$ and generative distributions $p_A(x)$ and $p_B(x)$ respectively iff,}
    \[ \mathbb{E}[f(x_i) | x_i \in A] - \mathbb{E}[f(x_i) | x_j \in B] = 0 \]
    \[ \int \left(p_A(x) - p_B(x) \right) f(x) dx =  0 \]

This work addresses an important issue raised in \cite{dwork2018group}, a model that satisfies conditional parity with respect to race and gender independently may fail to satisfy conditional parity with respect to the conjunction of race and gender. In the social science literature concerns about, potentially discriminated against, sub-demographics are referred to as intersectionality \cite{mccall2008complexity}. More formally, this work proposes a simple approach to ensure group fairness in expectation across an arbitrary set of subgroups. Applied to the popular decision trees, it is shown that provided the number of parity conditions is negligible compared to the number of training points, the order of computational and memory complexity is not increased.

\begin{figure}[b]
\label{fig:kernel}
\begin{center}
\begin{minipage}[b]{0.45\linewidth}\centering
\[K = \begin{tiny}\begin{bmatrix}
    \int \int q(x) K(x, x') q(x') \textnormal{d}x \textnormal{d}x' & \int q(x) K(x, x_0)  \textnormal{d}x & \int q(x) K(x, x_1)  \textnormal{d}x\\
     \int q(x) K(x, x_0)  \textnormal{d}x & K(x_0, x_0) & K(x_1, x_0) \\
     \int q(x) K(x, x_1)  \textnormal{d}x & K(x_1, x_0) & K(x_1, x_1)\\
\end{bmatrix}\end{tiny} \rightarrow \]
\vspace{8mm}
\end{minipage}
\hspace{20mm}
\begin{tikzpicture}[scale=.30]
  \begin{scope}
    \filldraw[fill=black!20!white] (1,8.) rectangle (10,10);
    \filldraw[fill=black!20!white] (0,0.) rectangle (2,9);

    \filldraw[fill=black!60!white] (0,8.) rectangle (2,10);

    \draw[pattern=north west lines, pattern color=black] (2,8.) rectangle (3,7);
    \draw[pattern=north west lines, pattern color=black] (3,7) rectangle (4,6);
    \draw[pattern=north west lines, pattern color=black] (4,6.) rectangle (5,5.);
    \draw[pattern=north west lines, pattern color=black] (5,5.) rectangle (6,4);
    \draw[pattern=north west lines, pattern color=black] (6,4) rectangle (7,3);
    \draw[pattern=north west lines, pattern color=black] (7,3.) rectangle (8,2.);
    \draw[pattern=north west lines, pattern color=black] (8,2.) rectangle (9,1);
    \draw[pattern=north west lines, pattern color=black] (9,1) rectangle (10,0);

    \draw (0, 0) grid (10, 10);
    \draw[very thick, scale=3] (3.333, 2.6666) -- (0, 2.6666);
    \draw[very thick, scale=3] (0.666, 0) -- (0.666, 3.333);
\end{scope}
\end{tikzpicture}

\end{center}
\caption{The above figure visualizes how the kernel matrix of a regressor can incorporate quadrature constraints. In the equations, $q(x)$ denotes the difference between two sub-populations $p_A(x)$ and $p_B(x)$. In the visual matrix representation (right) there are three components; the diagonal sub-matrix represents the leaves of the decision tree being independent of one another, the grey cells represent the relationship between the leaves and the constraints (denoted as $z$ in this work) and, finally, the dark upper left cells represents the interactions between constraints (denoted as $\bar{z}^T\bar{z}$ in this work).}
\end{figure}

\section{Constrained Kernel Methods}

As shown in \cite{aistats_submission}, kernel regressors may be constrained in terms of their expectation by adding auxiliary noiseless quadrature observations. Take for instance a Gaussian distribution with two dimensions. Given the distribution is zero mean, without loss in generality, correlation $\rho$ and variance  $\sigma_a^2$ and $\sigma_b^2$ respectively. With independent identically distributed noise $\sigma_n^2$, we can constrain the values of the expected outcomes by incorporating a noiseless observation on $\mu_a - \mu_b$ as follows,

\[
K = \begin{bmatrix}
    \sigma_a^2 + \sigma_b^2 - 2\rho\sigma_a\sigma_b & \sigma_a^2 - \rho\sigma_a\sigma_b & \sigma_b^2 - \rho\sigma_a\sigma_b\\
     \sigma_a^2 - \rho\sigma_a\sigma_b & \sigma_n^2 + \sigma_a^2 & \rho\sigma_a\sigma_b\\
     \sigma_b^2 - \rho\sigma_a\sigma_b& \rho\sigma_a\sigma_b& \sigma_n^2 +\sigma_b^2\\
\end{bmatrix}
\]

The above covariance matrix, $K$, has rank 2 with one perfectly observed value; namely the mean equality constraint. Thus inference on the two dimensions is constrained for $\sigma_n^2 >0$.

As shown visually in Figure 1, this principle can be extended to Gaussian processes and other kernel regression techniques by using the differences in quadrature observations in order to incorporate mean equality constraints. Multiple constraints can also be created by simply adding more columns and rows to the kernel matrix accordingly.

\section{Constrained Decision Trees}

While the above, constrained kernel inference is interesting, the widespread impact comes into play when we extend the result to decision tree regression, random forests, boosted trees and other ensemble techniques. This is said not to dismiss the importance or value in kernel methods more broadly, but rather due to their popularity amongst \emph{data scientists}, a profession more common than machine learning researchers \cite{wu2008top}.

While decision trees can be represented in either compressed or explicit kernel representation \cite{aistats_submission}, for the sake of conciseness this work will present results only for compressed representation. Thus, we will endeavor to minimize the perturbations induced on a per leaf bases, irrespective of the number of data points per leaf. The core difference between single and multiple constraints is that we can no longer use the arrowhead matrix lemma, instead, we must work out the update using the block matrix inversion lemma. Importantly, $p_A(x)$ and $p_B(x)$ for each constraint are defined as the empirical distributions of each subgroup considered. This is an important point as small subgroups may have empirical distributions which are not good approximations of the true generative distributions and hence our constrained space for inference may not constrain predictions to equate accordingly.

The kernel function in the compressed representation is simple the identity matrix, $K = \textnormal{I}$, that is to say each leaf is independent of one another. The kernel regression equation can be denoted as,

\[
f(\bar{\textnormal{x}}) = \mathbb{E}[\bar{\textnormal{y}}] = K_{\bar{x}, x} K_{x, x}^{-1} \textnormal{y}
\]

where the first $L$ values (number of leaves) of $K_{\bar{x}, x}$ indicate to which leaf $\bar{\textnormal{x}}$ belongs and the remaining values are set as zero, as the point to predict will not contribute to the empirical distributions of the subgroups under an inductive learning paradigm.

\begin{figure}[b!]
\centering
  \includegraphics[width=\textwidth]{ProPublica_multiple_train.png}
  \caption{The figure shows the effect of GFE constraints in the inferred scores of the ProPublica dataset between African American, Hispanic and all other defendants. Before perturbation is in blue and after in orange. Vertical lines indicate the mean of the distributions.}
  \label{fig:demo1}
\end{figure}

Using the block matrix inversion lemma we find,

\[K^{-1} = \begin{bmatrix}
    (1 + \sigma_n^2)\textnormal{I} & z\\
    z^T & \bar{z}^T\bar{z}\\
\end{bmatrix}^{-1} = {\tiny \begin{bmatrix}
    (1 + \sigma_n^2)^{-1} \textnormal{I} - (1 + \sigma_n^2)^{-2} z \left(\bar{z}^T\bar{z} - (1 + \sigma_n^2)^{-1} z^Tz \right)^{-1} z^T  & - (1 + \sigma_n^2)^{-1} z \left(\bar{z}^T\bar{z} - (1 + \sigma_n^2)^{-1} z^Tz \right)^{-1}\\
    - (1 + \sigma_n^2)^{-1} \left(\bar{z}^T\bar{z} - (1 + \sigma_n^2)^{-1} z^Tz \right)^{-1} z^T & \left(\bar{z}^T\bar{z} - (1 + \sigma_n^2)^{-1} z^Tz \right)^{-1} \\
\end{bmatrix}} \]

By simply inserting this into the kernel regression equation and noting that the elements of $\bar{z}$ are necessarily zero, the following update to the expected mean can be found as,

\[
f(\bar{\textnormal{x}}) = \frac{1}{1 + \sigma_n^2} \left( \textnormal{y}_j -  z_j \left(z^Tz \right)^{-1} z^T\textnormal{y} \right),
\]

with $z_j$ indicating the row of $z$ relating to the difference of subgroup distributions on leaf $j$. The effect of the noise can be removed by post multiplying by $1 + \sigma_n^2$.

\section{Experiments}

\subsection{ProPublica \& the COMPAS System}

The first experiment reproduces the experiment in \cite{aistats_submission} which uses a random forest to estimate the recidivism decile scores of the COMPAS algorithm applied to the ProPublica dataset while adding a GFE constraint between African Americans and Non-African Americans. However, it can also be noted that Hispanics also receive a similar discrimination. Figure 2, visualizes the effect of GFE constraints on the predicted distributions of the three demographics.

\subsection{Illinois State Employee Salaries}

The Illinois state employee salaries\footnote{\url{https://data.illinois.gov/datastore/dump/1a0cd05c-7d17-4e3d-938d-c2bfa2a4a0b1}} since 2011 can be seen to have a gender bias and bias between veterans and non-veterans. The motivation of this experiment was if one wished to predict a fair salary for future employees based on current staff. Gender labels were inferred using the employees' first names, parsed through the \emph{gender-geusser} python library. GFE constraints were applied between all intersections of gender and veteran / non-veterans, the marginals of gender and the marginals of veteran / non-veterans. Table 1 shows the expected outcome of each group before and after GFE constraints are applied and Figure 3 visualizes the perturbations to the marginals of each demographic intersection due to the GFE constraints. The train-test split was set as 80\%-20\% and the incorporation of the GFE constraints increase the root mean squared error from \$12,086 to \$12,772, the cost of fairness.

\begin{figure}[t!]
\centering
  \includegraphics[width=\textwidth]{illinois_intersectionality.png}
  \caption{The figure visualizes the distribution of salaries before and after perturbations due to GFE constraints. It is clear that female veterans benefit the most from such a constraint, while male non-veterans lose out. Colors and lines denoting the same meaning as Figure 2. The figures are cropped to the main mode of salaries to facilitate visual comparisons.}
  \label{fig:demo2}
\end{figure}

\begin{table}[t!]
\hspace{-13mm}{\small
\begin{tabular}{c| c c c c || c c|| c c }
\toprule
Group& Female Non-Vet. & Male Non-Vet. & Female Vet. & Male Vet. & Male & Female & Vet. &Non-Vet.\\
\hline
Original&47,334&52,777&41,890&51,063&46,962&52,215&49,555&49,805\\
Perturbed&\emph{48,695} & \emph{48,693}& \emph{48,694} & \emph{48,693}&\emph{48,695}&\emph{48,698}&\emph{48,775}&\emph{48,775}\\
\bottomrule
\end{tabular}
}
\caption{The above table shows the expected outcome of a random tree regressor with and without GFE constraints applied to four sub-demographics, between gender and between veterans and non-veterans.}
\end{table}

\section{Conclusion}

Regulatory bodies have shown precedent in developing affirmative action and other group fairness policy. This work extends previous efforts to develop group fairness constrained machine learning techniques. While relatively simple to understand and easy to incorporate into models used by practitioners, the methodology of this paper has a direct impact to four of the ten top data science algorithms according to \cite{wu2008top}. All source code used in the experiments are available at \url{https://github.com/OxfordML/Fair_Regression.git}.

\medskip

\small

\bibliographystyle{abbrv}


\section{Introduction}

Principal component analysis (PCA) is a fundamental technique for dimension reduction used widely in data analysis.
PCA projects data along a few directions that explain most of the variance of observed data. One can also view this as linearly transforming the original set of variables into a (smaller) set of uncorrelated variables called principal components.

Recent work in high-dimensional statistics has focused on sparse principal component analysis (SPCA), as ordinary PCA estimates become inconsistent in this regime \cite{johnstone2012consistency}. In SPCA, we restrict the principal components to be sparse, meaning they have only a few nonzero entries in the original basis. This has the advantage, among others, that the components are more interpretable \cite{jolliffe2003modified, zou2006sparse}, while components may no longer be uncorrelated. We study SPCA under the Gaussian $\emph{(single) spiked covariance model}$ introduced by \cite{johnstone2001distribution}: we observe $n$ samples of a random variable $X$ distributed according to a Gaussian distribution $\mathcal{N}(0, I_d+\theta uu^\top)$, where $||u||_2=1$ with at most $k$ nonzero entries,\footnote{Sometimes we will write this latter condition as $u \in B_0(k)$ where $B_0(k)$ is the $``\ell_0$-ball'' of radius $k$.} $I_d$ is the $d\times d$ identity matrix, and $\theta$ is the signal-to-noise parameter.
We study two settings of the problem, hypothesis testing and support recovery.

Sparsity assumptions have played an important role in a variety of other problems in high-dimensional statistics, in particular linear regression. Linear regression is also ill-posed in high dimensions, so via imposing sparsity on the regression vector we recover tractability.

Though the literature on two problems are largely disjoint, there is a striking similarity between the two problems, in particular when we consider statistical and computational trade-offs. The natural information-theoretically optimal algorithm for SPCA \cite{berthet2013optimal} involves searching over all possible
supports of the hidden spike. This bears resemblance to the minimax optimal algorithm for SLR \cite{raskutti2011minimax}, which optimizes over all sparse supports of the regression vector. Both problems appear to exhibit gaps between statistically optimal algorithms and known computationally efficient algorithms, and conditioned on relatively standard complexity assumptions, these gaps seem irremovable \cite{berthet2013complexity, wang2016statistical, zhang2014lower}.

\subsection{Our contributions}

In this paper we give algorithmic evidence that this similarity is likely not a coincidence. Specifically, we give a simple, general, and efficient procedure for transforming a black-box solver for sparse linear regression to an algorithm for SPCA. At a high level, our algorithm tries to predict each coordinate\footnote{From here on, we will use ``coordinate'' and ``variable'' interchangeably.} linearly from the rest of the coordinates using a black-box algorithm for SLR. The advantages of such a black-box framework are two fold: theoretically, it highlights a structural connection between the two problems;
practically, it allows one to simply plug in any of the vast majority of solvers available for SLR and directly get an algorithm for SPCA with provable guarantees. In particular,

\begin{itemize}
\item For hypothesis testing: we match state of the art provable guarantee for computationally efficient algorithms; our algorithm successfully distinguishes between isotropic and spiked Gaussian distributions at signal strength $\theta \gtrsim \sqrt{\frac{k^2 \log d}{n}}$. This matches the phase transition of diagonal thresholding \cite{johnstone2012consistency} and Minimal Dual Perturbation \cite{berthet2013optimal} up to constant factors.

\item For support recovery: for general $p$ and $n$, when each non-zero entry of $u$ is at least $\Omega(1/\sqrt{k})$ (a standard assumption in the literature), our algorithm succeeds with high probability for signal strength $\theta \gtrsim \sqrt{\frac{k^2 \log d}{n}}$, which is nearly optimal.\footnote{In the scaling limit $d/n \rightarrow \alpha$ as $d,n \rightarrow \infty$, the covariance thresholding algorithm \cite{deshpande2014sparse} theoretically succeeds at a signal strength that is an order of $\sqrt{\log d}$ smaller. However, our experimental results indicate that with an appropriate choice of black-box, our $\mathsf{Q}$ algorithm outperforms covariance thresholding}

\item In experiments, we demonstrate that using popular existing SLR algorithms as our black-box results in reasonable performance.

\item We theoretically and empirically illustrate that our SPCA algorithm is also robust to \emph{rescaling} of the data, for instance by using a Pearson correlation matrix instead of a covariance matrix.\footnote{Solving SPCA based on the correlation matrix was suggested in a few earlier works \cite{zou2006sparse, vu2013fantope}.} Many iterative methods rely on initialization via first running diagonal thresholding, which filters variables with higher variance; rescaling renders diganoal thresholding useless, so in some sense our framework is more robust.
\end{itemize}

\section{Preliminaries}
\label{sec:prelim}

\subsection{Problem formulation for SPCA}

\textbf{Hypothesis testing}
Here, we want to distinguish whether $X$ is distributed according to an isotropic Gaussian or a spiked covariance model. That is, our null and alternate hypotheses are:
\begin{align*}
H_0: X \sim \mathcal{N}(0, I_d) \text{ and } H_1: X \sim \mathcal{N}(0, I_d+\theta uu^\top),
\end{align*}

Our goal is to design a test $\psi\colon \mathbb{R}^{n \times d} \rightarrow \{0,1\}$ that discriminates  $H_0$ and $H_1$. More precisely,
we say that $\psi$ discriminates between $H_0$ and $H_1$ with probability $1-\delta$ if both type I and II errors have a probability smaller than $\delta$:
\[ \mathbf{P}_{H_0}(\psi(X) = 1)\leq \delta \text{ and } \mathbf{P}_{H_1}(\psi(X) = 0)\leq \delta.\]

We assume the following additional condition on the spike $u$:
\begin{enumerate}[label=(C\arabic*)]
    \item\label{cond1} $\nicefrac{c_{min}^2}{k} \le u_i^2 \le 1-\nicefrac{c_{min}^2}{k}$ for at least one $i \in [d]$ where $c_{min} > 0$ is some constant.
\end{enumerate}

The above condition says that at least one coordinate has enough mass, yet the mass is not entirely concentrated on just that singlecoordinate.
Trivially, we always have at least one $i\in[d]$ s.t. $u_i^2 \ge \nicefrac{1}{k}$, but this is not enough for our regression setup, since we want at least one other coordinate $j$ to have sufficient correlation with coordinate $i$.
We remark that the above condition is a very mild technical condition. If it were violated, almost all of the mass of $u$ is on a single coordinate, so a simple procedure for testing the variance (which is akin to diagonal thresholding) would suffice.

\textbf{Support recovery}
The goal of support recovery is to identify the support of $u$ from our samples $\mathbf{X}$. More precisely, we say that a support recovery algorithm succeeds if the recovered support $\widehat{S}$ is the same as $S$, the true support of $u$. As standard in the literature \cite{amini2009high,meinshausen2006high}, we need to assume a minimal bound on the size of the entries of $u$ in the support.

 For our support recovery algorithm, we will assume the following condition (note that it implies Condition \ref{cond1} and is much stronger):
\begin{enumerate}[label=(C\arabic*)]
\setcounter{enumi}{1}
    \item\label{cond2} $|u_i| \ge c_{min}/\sqrt{k}$ for some constant $0 < c_{min} < 1$ $\forall i \in [d]$
\end{enumerate}

Though the settings are a bit different, this minimal bound along with our results are  consistent with lower bounds known for sparse recovery. These lower bounds (\cite{fletcher2009necessary, wainwright2007information}; bound of \cite{fletcher2009necessary} is a factor of $k$ weaker) imply that the number of samples must grow roughly as $n \gtrsim (\nicefrac{1}{u_{min}^2}) k \log d$ where $u_{min}$ is the smallest entry of our signal $u$ normalized by $1/\sqrt{k}$, which is qualitativley the same threhold required by our theorems.

\subsection{Background on SLR}
In linear regression, we observe a response vector $y \in \mathbb{R}^n$ and a design matrix $\mathbb{X} \in \mathbb{R}^{n\times d}$ that are linked by the linear model $y = \mathbb{X}\beta^\ast + w$, where  $w \in \mathbb{R}^n$ is some form of observation noise, typically with i.i.d. $\mathcal{N}(0,\sigma^2)$ entries. Our goal is to recover $\beta^\ast$ given noisy observations $y$.
While the matrices $\mathbb{X}$ we consider arise from a (random) correlated design (as analyzed in \cite{wainwright2007information}, \cite{wainwright2009sharp}), it will make no difference to assume the matrices are deterministic by conditioning, as long as the distribution of the design matrix and noise are independent, which we will demonstrate in our case.  Most of the relevant results on sparse linear regression pertain to deterministic design.

In sparse linear regression, we additionally assume that $\beta^\ast$ has only $k$ non-zero entries, where $k \ll d$. This makes the problem well posed in the high-dimensional setting. Commonly used performance measures for SLR are tailored to prediction error ($\nicefrac{1}{n}\|\mathbb{X}\beta^\ast - \mathbb{X}\widehat{\beta}\|_2^2$ where $\widehat{\beta}$ is our guess), support recovery (recovering support of $\beta^\ast$), or parameter estimation (minimizing $\|\beta^\ast - \widehat{\beta}\|$ under some norm). We focus on prediction error, analyzed over random realizations of the noise. There is a large amount of work on SLR and we defer a more in-depth overview to Appendix \ref{app:slr}.

Most efficient methods for SLR impose certain conditions on $\mathbb{X}$. We focus on the \emph{restricted eigenvalue} condition, which roughly stated makes the prediction loss strongly convex near the optimum:
\begin{definition}[Restricted eigenvalue \cite{zhang2014lower}]
\label{def:re}
First define the cone $\mathbb{C}(S) = \{\beta \in \mathbb{R}^d \; | \; \|\beta_{S^c}\|_1 \le 3 \| \beta_S \|_1 \}$,
where $S^c$ denotes the complement, $\beta_T$ is $\beta$ restricted to the subset $T$.
The \emph{restricted eigenvalue} (RE) constant of $\mathbb{X}$, denoted $\gamma(\mathbb{X})$, is defined as the largest constant $\gamma > 0$ s.t.
  \[\nicefrac{1}{n}\|\mathbb{X}\beta\|_2^2 \ge \gamma \|\beta\|_2^2 \;\;\;\; \text{for all $\beta \in \bigcup_{|S|=k, S \subseteq [d]} \mathbb{C}(S)$} \]
\end{definition}
For more discussion on the restricted eigenvalue, see Appendix \ref{app:slr}.

\textbf{Black-box condition}
Given the known guarantees on SLR, we define a condition that is natural to require on the guarantee of our SLR black-box, which is invoked as $\mathsf{SLR}(y,\mathbb{X},k)$.
\begin{condition}[Black-box condition]
\label{blackcond}
Let $\gamma(\mathbb{X})$ denote the restricted eigenvalue of $\mathbb{X}$.
There are universal constants $c, c', c''$ such that
$\mathsf{SLR}(y,\mathbb{X},k)$ outputs $\widehat{\beta}$ that is $k$-sparse and satisfies:
\[\frac{1}{n} \|\mathbb{X}\widehat{\beta} - \mathbb{X}\beta^\ast\|_2^2 \le \frac{c}{\gamma(\mathbb{X})^2}  \frac{(\sigma^2 k \log d)}{n}  \;\;\;\;\; \text{$\forall\beta^\ast \in \Bz{k}$ } \text{w.p. $\ge 1 - c'\exp(-c''k\log d)$}\]
\end{condition}

\section{Algorithms and main results}
\label{sec:algo}

We first discuss how to view samples from the spiked covariance model in terms of a linear model. We then give some intuition motivating our statistic. Finally, we state our algorithms and main theorems, and give a high-level sketch of the proof.

\subsection{The linear model}
\label{sec:linearmodel}
Let $X^{(1)},X^{(2)},\ldots,X^{(n)}$ be $n$ i.i.d. samples from the spiked covariance model; denote as $\mathbf{X} \in \mathbb{R}^{n \times d}$ the matrix whose rows are $X^{(i)}$. Intuitively, if variable $i$ is contained in the support of the spike, then the rest of the support should allow to provide a nontrivial prediction for $\mathbf{X}_i$ since variables in the support are correlated. Conversely, for $i$ not in the support (or under the  isotropic null hypothesis), all of the variables are independent and other variables are useless for predicting $\mathbf{X}_i$.  So we regress $\mathbf{X}_i$ onto the rest of the variables.

Let $\mathbf{X}_{-i}$ denote the matrix of samples in the SPCA model with the $i$th column removed.  For each column $i$, we can view our data as coming from a linear model with design matrix $\mathbb{X} =\mathbf{X}_{-i}$ and the response variable $y= \mathbf{X}_i$.

The ``true'' regression vector depends on $i$. Under the alternate hypothesis $H_1$, if $i \in S$,
we can write $y = \mathbb{X} \beta^\ast + w$ where $\beta^\ast = \frac{\theta u_i}{1+(1-u_i^2)\theta} u_{-i}$
 and $w \sim \mathcal{N}(0,\sigma^2)$ with $\sigma^2=1 + \frac{\theta u_i^2}{1 + (1-u_i^2)\theta}$.\footnote{By the theory of linear minimum mean-square-error (LMMSE) confirms that this choice of $\beta^\ast$ minimizes the error $\sigma^2$. See Appendix \ref{app:lmmse}, \ref{app:calc} for details of this calculation.}
If $i \not\in S$, and for any $i \in [d]$ under the null hypothesis, $y = w$ where $w = \mathbf{X}_i \sim \mathcal{N}(0,1)$ (implicitly $\beta^\ast = \mathbf{0}$).

\subsection{Designing the test statistic}
Based on the linear model above, we want to compute a test statistic that will indicate when a coordinate $i$ is on support.
Intuitively, we predictive power of our linear model should be higher when $i$ is on support.
Indeed, a calculation shows that the variance in $X_i$ is reduced by approximately $\theta^2/k$. We want to measure this reduction in noise to detect when $i$ is on support or not.

Suppose for instance that we have access to $\beta^\ast$ rather than $\widehat{\beta}$ (note that this is not possible in practice since we do not know the support!).
Since we want to measure the reduction in noise when the variable is on support, as a first step we might try the following statistic:
  \[ Q_i = \nicefrac{1}{n}\|y - \mathbb{X} \beta^\ast \|_2^2 \]
Unfortunately, this statistic will not be able to distinguish the two hypotheses, as the reduction in the above error is too small (on the order of $\theta^2/k$ compared to overall order of $1+\theta$), so deviation due to random sampling will mask any reduction in noise. We can fix this by adding the variance term $\|y\|^2$:
\[ Q_i = \nicefrac{1}{n}\|y\|_2^2 - \nicefrac{1}{n}\|y - \mathbb{X} \beta^\ast\|_2^2 \]

On a more intuitive level, including $\|y\|_2^2$ allows us to measure the relative \emph{gain} in predictive power without being penalized by a possibly large variance in $y$. Fluctuations in $y$ due to noise will typically be canceled out in the difference of terms in $Q_i$, minimizing the variance of our statistic.

We have to add one final fix to the above estimator. We obviously do not have access to $\beta^\ast$, so we must use the estimate  $\widehat{\beta} = \mathsf{SLR}(y,\mathbb{X},k)$ ($y,\mathbb{X}$ are as defined in Section \ref{sec:linearmodel}) which we get from our black-box. As our analysis shows, this substitution does not affect much of the discriminative power of $Q_i$ as long as the SLR black-box satisfies prediction error guarantees stated in Condition \ref{blackcond}. This gives our final statistic:\footnote{As pionted out by a reviewer, Note that this statistic is actually equivalent to $R^2$ up to rescaling by sample variance. Note that our formula is slightly different though as we use the sample variance computed with population mean as opposed to sample mean, as the mean is known to be zero.}
\[ Q_i = \nicefrac{1}{n}\|y\|_2^2 - \nicefrac{1}{n}\|y - \mathbb{X} \widehat{\beta} \|_2^2. \]

\subsection{Algorithms}
Below we give algorithms for hypothesis testing and for support recovery, based on the $Q$ statistic:
\begin{figure}[!ht]
\begin{minipage}[t]{6cm}

  \begin{algorithm}[H]
    \caption{$Q$-hypothesis testing}
    \label{algohyp}
    \begin{algorithmic}
      \STATE Input: $\mathbf{X} \in \mathbb{R}^{d \times n} , k$
      \STATE Output: $\{0,1\}$
      \FOR{ $i = 1, \dots, d$}
      \STATE $\widehat{\beta}_i = \mathsf{SLR}(\mathbf{X}_i, \mathbf{X}_{-i}, k) $
      \item $Q_i =\frac{1}{n}\|\mathbf{X}_i\|_2^2-\frac{1}{n}\|\mathbf{X}_i-\mathbf{X}_{-i}\widehat{\beta}_i\|_2^2$
      \IF{ $Q_i>\frac{13k\log{\frac{d}{k}}}{n}$} \STATE return $1$
      \ENDIF
      \ENDFOR
      \STATE Return $0$
      \end{algorithmic}
      \end{algorithm}
    \end{minipage}%
\hfill
    \begin{minipage}[t]{6cm}
     \null
\begin{algorithm}[H]
    \caption{$Q$-support recovery}
    \label{algosupp}
    \begin{algorithmic}
      \STATE Input: $\mathbf{X} \in \mathbb{R}^{d \times n} , k$
      \STATE $\widehat{S}=\varnothing$
      \FOR{ $i = 1,\ldots, d$}
      \STATE $\widehat{\beta}_i = \mathsf{SLR}(\mathbf{X}_i, \mathbf{X}_{-i}, k) $
      \item $Q_i = \frac{1}{n} \|\mathbf{X}_i\|_2^2-\frac{1}{n}\|\mathbf{X}_i-\mathbf{X}_{-i}\widehat{\beta}_i\|_2^2$
      \IF{ $Q_i>\frac{13k \log {\frac{d}{k}}}{n}$} \STATE $\widehat{S}\leftarrow \widehat{S}\cup\{i\}$
      \ENDIF
      \ENDFOR
      \STATE Return $\widehat{S}$
      \end{algorithmic}
      \end{algorithm}
    \end{minipage}
\hfill
\centering
\end{figure}

\bigskip

Below we summarize our guarantees for the above algorithms. The proofs are simple, but we defer them to Appendix \ref{app:analysis}.

\begin{theorem}[Hypothesis test]
\label{thm:hyptest}
Assume we have access to $\mathsf{SLR}$ that satisfies Condition \ref{blackcond} and with runtime $T(d,n,k)$ per instance.
Under Condition \ref{cond1},
there exist universal constants $c_1,c_2,c_3,c_4$ s.t.
if $\theta^2 > \frac{c_1}{c_{min}^2} \frac{k^2\log d}{n}$ and $n > c_2 k \log d$,
Algorithm 1 outputs $\psi$ s.t.
\[\mathbf{P}_{H_0}(\psi(X) = 1) \vee \mathbf{P}_{H_1}(\psi(X) = 0) \le c_3\exp(-c_4 k \log d)\]
in time $O(dT + d^2n)$.
\end{theorem}

\begin{theorem}[Support recovery]
\label{thm:supp}
Under the same conditions as above plus Condition \ref{cond2}, if $\theta^2 > \frac{c_1}{c_{min}^2} \frac{k^2\log d}{n}$,
Algorithm 2 above finds
$\widehat{S} = S$ with probability at least $1-c_3\exp(-c_4 k \log d)$ in time $O(dT + d^2n)$.
\end{theorem}

\subsection{Comments}

\textbf{RE for sample design matrix}
Because population covariance $\Sigma = \mathbb{E}[\mathbb{X}\X^\top]$ has minimum eigenvalue 1, with high probability the sample design matrix $\mathbb{X}$ has constant restricted eigenvalue value given enough samples, i.e. $n$ is large enough (see Appendix \ref{app:prop} for more details), and the prediction error guarantee of Condition \ref{blackcond} will be good enough for our analysis.

\textbf{Running time}
The runtime of both Algorithm \ref{algohyp} and \ref{algosupp} is $\tilde{O}(nd^2)$. The discussion presented at the end of Appendix
\ref{app:analysis} details why this is competitive for (single spiked) SPCA, at least theoretically.

\textbf{Unknown sparsity}
Throughout the paper we assume that the sparsity level $k$ is known. However, if $k$ is unknown, standard techniques could be used to adaptively find approximate values of $k$ (\cite{do2008sparsity}). For instance, for hypothesis testing, we can start with an initial overestimate $k'$, and keep halving until we get enough coordinates $i$ with $Q_i$ that passes the threshold for the given $k'$.

\subsection*{Robustness of Q statistic to rescaling}
\label{sec:rescale}

Intuitively, our algorithms for detecting correlated structure in data should be invariant to rescaling of the data; the precise scale or units for which one variable is measured should not have an impact on our ability to find meaningful structure underlying the data.
Our algorithms based on $Q$ are robust to rescaling, perhaps unsurprisingly, since correlations between variables in the support remain under rescaling.

On the other hand, diagonal thresholding, an often-used preprocessing step for SPCA which filters variables strictly based on variance, would trivially fail under rescaling. This illustrates a strength of our framework over other existing algorithms for SPCA.

Below we show explicitly that $Q$ statistics are indeed robust to rescaling:
Let $\widetilde{X} = DX$ be the rescaling of $X$, where $D$ is some diagonal matrix. Let $D_{S}$ be $D$ restricted to rows and columns in $S$.
Note that $\widetilde{\Sigma}$, the covariance matrix of the rescaled data, is just $D \Sigma D$ by expanding the definition. Similarly, note
$\widetilde{\Sigma}_{2:d,1} = D_{1} D_{2:d}  \Sigma_{2:d,1}$ where $D_{2:d}$ denotes $D$ without row and column $1$.
Now, recall the term which dominated our analysis of $Q_i$ under $H_1$, $(\beta^\ast)^\top \Sigma_{2:d} \beta^\ast$,
which was equal to
\[ \Sigma_{1,2:d} \Sigma_{2:d}^{-1} \Sigma_{2:d,1} \]

We replace the covariances by their rescaled versions to obtain:
\[ \tilde{\beta^\ast}^\top \widetilde{\Sigma} \tilde{\beta^\ast}
= (D_{1} \Sigma_{1,2:d} D_{2:d}) D_{2:d}^{-1} \Sigma_{2:d}^{-1} D_{2:d}^{-1} (D_{2:d} \Sigma_{2:d,1} D_{1} ) =   D_{1}^2  (\beta^\ast)^\top \Sigma_{2:d} \beta^\ast \]

For the spiked covariance model, rescaling variances to one amount to rescaling with $D_{1} = \nicefrac{1}{(1+\theta)}$. Thus, we see that our signal strength is affected only by constant factor (assuming $\theta \le 1$).

\section{Experiments}
\label{sec:exp}
We test our algorithmic framework on randomly generated synthetic data and compare to other existing algorithms for SPCA. The code was implemented in Python using standard libraries.

We refer to our general algorithm from Section \ref{sec:algo} that uses the $Q$ statistic as SPCAvSLR.
For our SLR ``black-box,'' we use thresholded Lasso \cite{zhang2014lower}.\footnote{Thresholded Lasso is a variant of lasso where after running lasso, we keep the $k$ largest (in magnitude) coefficients to make the estimator k-sparse. Proposition 2 in \cite{zhang2014lower} shows that thresholded Lasso satisfies Condition \ref{blackcond}}. (We experimented with other SLR algorithms such as the forward-backward algorithm of \cite{foba} and CoSaMP\footnote{CoSaMP in theory requires the stronger condition of restricted isometry on the ``sensing'' matrix.} \cite{cosamp}, but results were similar and only slower.)

For more details on our experimental setup, including hyperparameter selection, see Appendix \ref{app:experiment}.

\textbf{Support recovery}
We randomly generate a spike $u \in \mathbb{R}^d$ by first choosing a random support of size $k$, and then using random signs for each coordinate (uniformity is to make sure Condition \ref{cond2} is met). Then spike is scaled appropriately with $\theta$ to build the spiked covariance matrix of our normal distribution, from which we draw samples.

We study how the performance of six algorithms vary over various values of $k$ for fixed $n$ and $d$.\footnote{We remark that while the size of this dataset might seem too small to be representative ``high-dimensional'' setting, these are representative of the usual size of dataset that these methods are usually tested on. One bottleneck is the computation of the covariance matrix.} As in the \cite{deshpande2014sparse}, our measure is the fraction of true support.
We compare SPCAvSLR with the following algorithms: diagnoal thresholding, which is a simple baseline; ``SPCA'' (ZHT \cite{zou2006sparse}) is a fast heuristic also based on the regression idea; the truncated power method of \cite{yuan2013truncated}, which is known for both strong theoretical guarantees and empirical performance; covariance thresholding, which has state-of-the-art theoretical guarantees.

We modified each algorithm to return the top $k$ most likely coordinates in the support (rather than thresholding based on a cutoff);
for algorithms that compute a candidate eigenvector, we choose the top $k$ coordinates largest in absolute value.
\begin{figure}[!ht]
\includegraphics[scale=0.34]{supprec_d625_n625.png}
\includegraphics[scale=0.34]{supprec_d1250_n625.png}
\centering
\caption{Performance of diagonal thresholding, SPCA (ZHT), truncated power method, covariance thresholding, and SPCAvSLR for support recovery at $n=d=625$ (left) and $n=625,d=1250$ (right), varying values of $k$, and $\theta=3.0$. On the horizontal axis we show $k/\sqrt{n}$; the vertical axis is the fraction of support correctly recovered. Each datapoint on the figure is averaged over 50 trials.}
\setlength{\belowcaptionskip}{-20pt}
\end{figure}

We observe that SPCAvSLR performs better than covariance thresholding and diagonal thresholding, but its performance falls short of that of the truncated power method and the heuristic algorithm of \cite{zou2006sparse}.
We suspect the playing with different SLR algorithms may slightly improve its performance.
The reason for the gap between performance of SPCAvSLR and other state of the arts algorithms despite its theoretical guarantees is open to further investigation.

\textbf{Hypothesis testing}
We generate data under two different distributions: for the spiked covariance model, we generate a spike $u$ by sampling a uniformly random direction from the $k$-dimensional unit sphere, and embedding the vector at a random subset of $k$ coordinates among $d$ coordinates; for the null, we draw from standard isotropic Gaussian. In a single trial, we draw $n$ samples from each distribution and we compute various statistics\footnote{While some of the algorithms used for support recovery in the previous section could in theory be adapated for hypothesis testing, the extensions were immediate so we do not consider them here.} (diagonal thresholding (DT), Minimal Dual Perturbation (MDP), and our $Q$ statistic, again using thresholded Lasso). We repeat for $100$ trials, and plot the resulting empirical distribution for each statistic. We observe similar performance of DT and $Q$, while MDP seems slightly more effective at distinguishing $H_0$ and $H_1$ at the same signal strength (that is, the distributions of the statistics under $H_0$ vs. $H_1$ are more well-separated).

\textbf{Rescaling variables}
As discussed in Section \ref{sec:rescale}, our algorithms are robust to rescaling the covariance matrix to the correlation matrix.
As illustrated in Figure \ref{fig:rescaling} (right), DT fails while $Q$ appears to be still effective for distinguishing hypotheses the same regime of parameters.
Other methods such as MDP and CT also appear to be robust to such rescaling (not shown).
This suggests that more modern algorithms for SPCA may be more appropriate than diagonal thresholding in practice, particularly  on instances where the relative scales of the variables may not be accurate or knowable in advance, but we still want to be able to find correlation between the variables.

\begin{figure}[!ht]
\includegraphics[scale=0.21]{n200_d500_k30_t4_alp0-167.png}
\includegraphics[scale=0.21]{MDP-n200-d500-k30-t5.png}
\includegraphics[scale=0.21]{n200_d400_k30_t10_lasso_rescaled.png}
\centering
\caption{Performance of diagonal thresholding (D), MDP, and $Q$ for hypothesis testing at $n=200,d=500,k=30,\theta=4$ (left and center). T0 denotes the statistic $T$ under $H_0$, and similarly for T1.
The effect of rescaling the covariance matrix to make variances indistinguishable is demonstrated (right).}
\label{fig:rescaling}
\end{figure}

\section{Previous work}
Here we discuss in more detail previous approaches to SPCA and how it relates to our work. Various approaches to SPCA have been designed in an extensive list of prior work. As we cannot cover all of them, we focus on works that aim to give computationally efficient (i.e. polynomial time) algorithms with provable guarantees in settings similar to ours.

These algorithms include fast, heuristic methods based on $\ell_1$ minimization \cite{jolliffe2003modified, zou2006sparse}, rigorous but slow methods based on natural semidefinite program (SDP) relaxations \cite{d2007direct, amini2009high, vu2013fantope, wang2016statistical}, iterative methods motivated by power methods for approximating eigenvectors \cite{yuan2013truncated, journee2010generalized}, non-iterative methods based on random projections \cite{gataric2017sparse}, among others.
Many iterative methods rely on initialization schemes, such as ordinary PCA or diagonal thresholding \cite{johnstone2012consistency}.

Below, we discuss the known sample bounds for support recovery and hypothesis testing.

\textbf{Support recovery}
\cite{amini2009high} analyzed both diagonal thresholding and an SDP for support recovery under the spiked covariance model.\footnote{They analyze the subcase when the spike is uniform in all $k$ coordinates.}
They showed that the SDP requires an order of $k$ fewer samples when the SDP optimal solution is rank one.
However, \cite{krauthgamer2013semidefinite} showed that the rank one condition does not happen in general, particularly in the regime approaching the information theoretic limit ($\sqrt{n} \lesssim k \lesssim \nicefrac{n}{\log d}$). This is consistent with computational lower bounds from \cite{berthet2013complexity} ($k \gtrsim \sqrt{n}$),
but a small gap remains (diagonal thresholding and SDP's succeed only up to $k \lesssim \sqrt{n / \log d}$). The above gap was closed by the covariance thresholding algorithm,
first suggested by \cite{krauthgamer2013semidefinite} and analyzed by \cite{deshpande2014sparse}, that succeeds in the regime $\sqrt{n / \log d} \lesssim k \lesssim \sqrt{n}$, although the theoretical guarantee is limited to the regime when $d/n \rightarrow \alpha$ due to relying on techniques from random matrix theory.

\textbf{Hypothesis testing}
Some works \cite{berthet2013optimal, amini2009high, d2014approximation} have focused on the problem of detection. In this case, \cite{berthet2013optimal} observed that it suffices to work with the much simpler dual of the standard SDP called Minimal Dual Perturbation (MDP). Diagonal thresholding (DT) and MDP work up to the same signal threshold $\theta$ as for support recovery, but MDP seems to outperform DT on simulated data \cite{berthet2013optimal}.
MDP works at the same signal threshold as the standard SDP relaxation for SPCA. \cite{d2014approximation} analyze a statistic based on an SDP relaxation and its approximation ratio to the optimal statistic. In the regime where $k,n$ are proportional to $d$, their statistic succeeds at a signal threshold for $\theta$ that is independent of $d$, unlike the MDP. However, their statistic is quite slow to compute; runtime is at least a high order polynomial in $d$.

\textbf{Regression based approaches}
To the best of our knowledge, our work is the first to give a general framework for SPCA that uses SLR in a \emph{black-box} fashion. \cite{zou2006sparse} uses specific algorithms for SLR such as Lasso as a subroutine, but they use a heuristic alternating minimization procedure to solve a non-convex problem, and hence lack any theoretical guarantees.
\cite{meinshausen2006high} applies a regression based approach to a restricted class of graphical models. While our regression setup is similar, their statistic is different and their analysis depends directly on the particulars of Lasso. Further, their algorithm requires extraneous conditions on the data
 \cite{cai2013sparse} also uses a reduction to linear regression for their problem of sparse subspace estimation. Their iterative algorithm depends crucially on a good initialization done by a diagonal thresholding-like pre-processing step, which fails under rescaling of the data.\footnote{See Section \ref{sec:rescale} for more discussion on rescaling.} Furthermore, their framework uses regression for the specific case of orthogonal design, whereas our design matrix can be more general as long as it satisfies a condition similar to the restricted eigenvalue condition. On the other hand, their setup allows for more general $\ell_q$-based sparsity as well as the estimation of an entire subspace as opposed to a single component. \cite{ma2013sparse} also achieves this more general setup, while still suffering from the same initialization problem.

\textbf{Sparse priors}
Finally, connections between SPCA and SLR have been noted in the probabilistic setting \cite{koyejo2014prior, khanna2015sparse}, albeit in an indirect manner: the same sparsity-inducing priors can be used for either problem. We view our work as entirely different as we focus on giving a black-box reduction. Furthermore, provable guarantees for the EM algorithm and variational methods are lacking in general, and it is not immediately obvious what signal threshold their algorithm achieves for the single spike covariance model.

\section{Conclusion}
\label{sec:conclude}
We gave a black-box reduction for reducing instances of the SPCA problem under the spiked covariance model to instances
of SLR. Given oracle access to SLR black-box meeting a certain natural condition, the reduction is shown to efficiently
solve hypothesis testing and support recovery.

Several directions are open for future work. The work in this paper remains limited to the Gaussian setting and to the single spiked covariance model. Making the results more general would make the connection made here more appealing. Also, the algorithms designed here, though simple, seem a bit wasteful in that they do not aggregate information from different statistics. Designing a more efficient estimator that makes a more efficient use of samples would be interesting. Finally, there is certainly room for improvement by tuning the choice of the SLR black-box to make the algorithm more efficient for use in practice.

\pagebreak

\newpage

{\small
\bibliographystyle{plain}