\section{Introduction}
In this era of big data, deep neural networks (DNNs) have achieved great success in machine learning research and commercial applications. When large amounts of training data are available, the capacity of DNNs can  easily be increased by adding more units or layers to extract more effective high-level features \cite{Goodfellow-et-al-2016,He_2016_CVPR,Wan_2013_ICML}. However, big networks with millions of parameters can easily overfit even the largest of datasets. It is thus crucial to regularize DNNs so that they can extract ``meaningful'' features not only from the training data, but also from the test data.

Many widely-used network regularizations are data-independent. Such techniques include  weight decay, parameter sharing, \textit{DropOut} \cite{Srivastava_2014_dropout}, \textit{DropConnect} \cite{icml2013_wan13}, etc. Intuitively, weight decay alleviates overfitting by reducing the magnitude of the  weights and the features, and \textit{DropOut} and \textit{DropConnect} can be viewed as computationally inexpensive ways to train an exponentially large ensemble of DNNs. Their effectiveness as network regularizers can be quantified by analyzing the Rademacher complexity, which provides an upper bound for the generalization error \cite{bartlett1998sample,icml2013_wan13}.

\begin{figure*}
    \centering
    \begin{subfigure}[t]{0.22\textwidth}
        \includegraphics[width=\textwidth]{feature_mnist_original}
        \caption{10,000 test data from MNIST}
        \label{fig:mnist_original}
    \end{subfigure}
    ~

    \begin{subfigure}[t]{0.22\textwidth}
        \includegraphics[width=\textwidth]{feature_mnist_wd}
        \caption{Features learned from a network regularized by weight decay}
        \label{fig:mnist_wd}
    \end{subfigure}
    ~

    \begin{subfigure}[t]{0.22\textwidth}
        \includegraphics[width=\textwidth]{feature_mnist_dropout}
        \caption{Features learned from a network regularized by \textit{DropOut}}
        \label{fig:mnist_dropout}
    \end{subfigure}
    ~

    \begin{subfigure}[t]{0.22\textwidth}
        \includegraphics[width=\textwidth]{feature_mnist_ldm}
        \caption{Features learned from a network regularized by \textit{LDMNet}}
        \label{fig:mnist_ldm}
    \end{subfigure}
    \caption{Test data of MNIST and their features learned by the same network with different regularizers. All networks are trained from the same set of 1,000 images. Data are visualized in two dimensions using PCA, and ten classes are distinguished by different colors. In (b) and (c), the features learned by weight decay and \textit{DropOut} typically sample two dimensioanl regions. Whereas in (d), features learned by \textit{LDMNet} tend to concentrate on one-dimensional and zero-dimensional manifolds (curves and points). }\label{fig:visual_mnist}
\end{figure*}

\begin{figure*}
    \centering
    \begin{subfigure}[b]{0.25\textwidth}
        \includegraphics[width=\textwidth]{feature_vgg}
        \caption{VGG-face}
        \label{fig:vgg}
    \end{subfigure}
    ~

    \hspace{-5mm}
    \begin{subfigure}[b]{0.25\textwidth}
        \includegraphics[width=\textwidth]{feature_wd}
        \caption{Weight decay}
        \label{fig:wd}
    \end{subfigure}
    ~

    \hspace{-5mm}
    \begin{subfigure}[b]{0.25\textwidth}
        \includegraphics[width=\textwidth]{feature_dropout}
        \caption{\textit{DropOut}}
        \label{fig:dropout}
    \end{subfigure}
    ~

    \hspace{-5mm}
    \begin{subfigure}[b]{0.25\textwidth}
        \includegraphics[width=\textwidth]{feature_ldm}
        \caption{\textit{LDMNet}}
        \label{fig:ldm}
    \end{subfigure}
    \caption{Features generated from VGG-face model \cite{parkhi2015deep} with or without embedding learned from the network specified in Table~\ref{tab:nirvis_net} with different regularizations. A total of five subjects are displayed, with one color per subject. VIS images are denoted as filled circles and NIR images are denoted as unfilled diamonds. All features are visualized in two dimensions using PCA. In (a), (b), and (c), features of the same subject typically form two clusters, one for NIR and the other for VIS. In (d), with \textit{LDMNet}, features of the same subject from  different modalities merge to form a single low dimensional manifold.}\label{fig:visual}
\end{figure*}

Most of the data-dependent regularizations are  motivated  by the empirical observation that data of interest typically lie close to a manifold, an assumption that has previously assisted machine learning tasks such as nonlinear embedding \cite{nonlinear_embedding}, semi-supervised labeling \cite{belkin2006manifold}, and multi-task classification \cite{evgeniou2005learning}. In the context of DNN,  data-dependent regularization techniques include the tangent distance algorithm \cite{Simard_1993,Simard_2012}, tangent prop algorithm \cite{Simard_1992}, and manifold tangent classifier \cite{Rifai_2011}. Typically, these algorithms only focus on the geometry of the input data, and do not encourage the network to produce geometrically meaningful features. Moreover, it is typically hard to explicity parametrize the underlying manifold, and some of the algorithms require human input of the tangent planes or tangent vectors \cite{Goodfellow-et-al-2016}.

Motivated by the same manifold observation, we propose a new network regularization technique that focuses on the geometry of both the input data and the output features, and build low-dimensional-manifold-regularized neural networks (\textit{LDMNet}). This is inspired by a recent algorithm in image  processing, the \textit{low dimensional manifold model} \cite{ldmm,ldmm_wgl,ldmm_scientific}. The idea of \textit{LDMNet} is that the concatenation $(\bm{x}_i,\bm{\xi}_i)$ of the input data $\bm{x}_i$ and the output features $\bm{\xi}_i$ should sample a collection of low dimensional manifolds. This idea is loosely related to the Gaussian mixture model: instead of assuming that the generating distribution of the input data is a mixture of Gaussians, we assume that the input-feature tuples $(\bm{x}_i,\bm{\xi}_i)$ are generated  by a mixture of low dimensional manifolds. To emphasize this, we explicitly penalize the loss function with a term using an elegant formula from differential geometry to compute the dimension of the underlying manifold. The resulting variational problem is then solved via alternating minimization with respect to the manifold and the network weights. The corresponding Euler-Lagrange equation is a Laplace-Beltrami equation on a point cloud, which is efficiently solved by the point integral method (PIM) \cite{pim} with $O(N)$ computational complexity, where $N$ is the size of the input data. In \textit{LDMNet}, we never have to explicitly parametrize the manifolds or derive the tangent planes, and the solution is obtained by solving only the variational problem.

In the experiments, we demonstrate two benefits of \textit{LDMNet}: First, by extracting geometrically meaningful features, \textit{LDMNet}  significantly outperforms widely-used regularization techniques such as weight decay and \textit{DropOut}. For example, Figure~\ref{fig:visual_mnist} shows the two-dimensional projections of the test data from MNIST and their features learned from 1,000 training data by the same network with different regularizers. It can be observed that the features learned by weight decay and \textit{DropOut} typically sample two-dimensional regions, whereas the features learned by \textit{LDMNet} tend to lie close to one-dimensional and zero-dimensional manifolds (curves and points). Second, in some imaging problems, \textit{LDMNet} is more likely to find the model that is subject to different illumination patterns. By regularizing the network outputs, \textit{LDMNet}  can extract common features of the same subject imaged via different modalities so that these features sample the same low dimensional manifold. This can be observed in Figure~\ref{fig:ldm}, where the features of the same subject  extracted by \textit{LDMNet} from visible  (VIS) spectrum images and  near-infrared (NIR) spectrum images merge to form a single low dimensional manifold. This significantly increases the accuracy of cross-modality face recognition. The details of the experiments will be explained in Section \ref{sec:experiments}.

\section{Model Formulation}

For simplicity of explanation, we consider a $K$-way classification using a DNN. Assume that $\{(\bm{x}_i,y_i)\}_{i=1}^N\subset \mathbb{R}^{d_1} \times \{1,\ldots,K\}$ is the labeled training set, and $\bm{\theta}$ is the collection of network weights. For every datum $\bm{x}_i$ with class label $y_i \in \{1,\ldots,K\}$, the network first learns a $d_2$-dimensional feature $\bm{\xi}_i=f_{\bm{\theta}}(\bm{x}_i)\in\mathbb{R}^{d_2}$, and then applies a softmax classifier to obtain the probability distribution of $\bm{x}_i$ over the $K$ classes. The softmax loss $\ell(f_{\bm{\theta}}(\bm{x}_i),y_i)$ is then calculated for $\bm{x}_i$ as a negative log-probability of class $y_i$. The empirical loss function $J(\bm{\theta})$ is defined as the average loss on the training set:
\begin{align}
  \label{eq:loss}
  J(\bm{\theta}) = \frac{1}{N}\sum_{i=1}^N\ell(f_{\bm{\theta}}(\bm{x}_i),y_i).
\end{align}
When the training samples are scarce, statistical learning theories predict that overfitting to the training data  will occur \cite{Vapnik_1999}. What this means is that the average loss on the testing set can still be large even if the empirical loss $J(\bm{\theta})$ is trained to be small.

\textit{LDMNet} provides an explanation and a solution for network overfitting. Empirical observation suggests that many data of interest typically sample a collection of low dimensional manifolds, i.e. $\{\bm{x}_i\}_{i=1}^N \subset \mathcal{N} =\cup_{l=1}^L\mathcal{N}_l \subset \mathbb{R}^{d_1}$. One would also expect that  the feature extractor, $f_{\bm{\theta}}$, of a good learning algorithm  be a smooth function over $\mathcal{N}$ so that small variation in $\bm{x} \in \mathcal{N}$ would not lead to dramatic change in the learned feature $\bm{\xi} = f_{\bm{\theta}}(\bm{x})\in \mathbb{R}^{d_2}$. Therefore the concatenation of the input data and output features,  $\{(\bm{x}_i,\bm{\xi}_i)\}_{i=1}^N$, should sample a collection of low dimensional manifolds $\mathcal{M} = \cup_{l=1}^L\mathcal{M}_l \subset \mathbb{R}^{d}$, where $d = d_1+d_2$, and $\mathcal{M}_l=\left\{(\bm{x}, f_{\bm{\theta}}(\bm{x}))\right\}_{x \in \mathcal{N}_l}$ is the graph of $f_{\bm{\theta}}$ over $\mathcal{N}_l$. We suggest that network overfitting occurs when $\dim(\mathcal{M}_l)$ is too large after training. Therefore, to reduce overfitting,  we explicitly use the dimensions of $\mathcal{M}_l$ as a regularizer in the following variational form:
\begin{align}
  \label{eq:ldm}
 &\min_{\bm{\theta}, \mathcal{M}} \quad J(\bm{\theta}) + \frac{\lambda}{|\mathcal{M}|} \int_{\mathcal{M}}\dim(\mathcal{M}(\bm{p}))d\bm{p}\\ \nonumber
& \text{s.t. } \quad \{(\bm{x}_i,f_{\bm{\theta}}(\bm{x}_i))\}_{i=1}^N \subset \mathcal{M},
\end{align}
where for any $\bm{p} \in \mathcal{M} = \cup_{l=1}^L\mathcal{M}_l$, $\mathcal{M}(\bm{p})$ denotes the manifold $\mathcal{M}_l$ to which $\bm{p}$ belongs, and $|\mathcal{M}| = \sum_{l=1}^L|\mathcal{M}_l|$ is the volume of $\mathcal{M}$.
The following theorem from differential geometry provides an elegant way of calculating the manifold dimension in (\ref{eq:ldm}).
\begin{thm}
\label{thm:dim}
{\normalfont{\cite{ldmm}}} Let $\mathcal{M}$ be a smooth submanifold isometrically embedded in $\mathbb{R}^d$. For any $\bm{p}=(p_i)_{i=1}^d\in \mathcal{M}$,
  \begin{align*}
    \dim(\mathcal{M})=\sum_{i=1}^d\left|\nabla_\mathcal{M}\alpha_i(\bm{p})\right|^2,
  \end{align*}
where $\alpha_i(\bm{p})=p_i$ is the coordinate function, and $\nabla_\mathcal{M}$ is the gradient operator on the manifold $\mathcal{M}$. More specifically, $\nabla_\mathcal{M} \alpha_i = \sum_{s,t=1}^{k}g^{st}\partial_t\alpha_i\partial_s$, where $k$ is the intrinsic dimension of $\mathcal{M}$, and $g^{st}$ is the inverse of the metric tensor.
\end{thm}
\iffalse
\begin{figure}[t]
  \includegraphics[width=0.9\linewidth]{network_pic}
\caption{The structure of \textit{LDMNet} built upon a general DNN, where $\bm{x}_i$ is the input, $\bm{\xi}_i$ is the output feature, $\bm{\theta}$ is the network weights, and $J(\bm{\theta})$ is the original loss function. \textit{LDMNet} minimizes the new objective function $\tilde{J}(\bm{\theta})$, which is the original loss penalized with the dimensions of the manifolds $\{\mathcal{M}_l\}_{l=1}^L$ sampled by $\{(\bm{x}_i,\bm{\xi}_i)\}_{i=1}^N$.}
\label{fig:network_pic}
\end{figure}
\fi

As a result of Theorem \ref{thm:dim}, (\ref{eq:ldm}) can be reformulated as:
\begin{align}
  \label{eq:energy_alpha}
  &\min_{\bm{\theta}, \mathcal{M}} \quad J(\bm{\theta}) + \frac{\lambda}{|\mathcal{M}|} \sum_{j=1}^d \|\nabla_\mathcal{M}\alpha_j\|_{L^2(\mathcal{M})}^2 \\ \nonumber
  &\text{s.t. } \quad \{(\bm{x}_i,f_{\bm{\theta}}(\bm{x}_i))\}_{i=1}^N \subset \mathcal{M}
\end{align}
where $\sum_{j=1}^d \|\nabla_\mathcal{M}\alpha_j\|_{L^2(\mathcal{M})}^2$ corresponds to the $L^1$ norm of the local dimension. To solve (\ref{eq:energy_alpha}), we alternate the direction of minimization with respect to $\mathcal{M}$ and $\bm{\theta}$. More specifically, given $(\bm{\theta}^{(k)},\mathcal{M}^{(k)})$ at step $k$ satisfying $\{(\bm{x}_i,f_{\bm{\theta}^{(k)}}(\bm{x}_i))\}_{i=1}^N \subset \mathcal{M}^{(k)}$, step $k+1$ consists of the following
\begin{itemize}
\item Update $\bm{\theta}^{(k+1)}$ and the perturbed coordinate functions $\bm{\alpha}^{(k+1)}=(\alpha_1^{(k+1)},\cdots,\alpha_d^{(k+1)})$ as the minimizers of (\ref{eq:update_two}) with the fixed manifold $\mathcal{M}^{(k)}$:
  \begin{align}\label{eq:update_two}
    &\min_{\bm{\theta}, \bm{\alpha}} \quad J(\bm{\theta})+\frac{\lambda}{|\mathcal{M}^{(k)}|} \sum_{j=1}^d \|\nabla_{\mathcal{M}^{(k)}}\alpha_j\|_{L^2(\mathcal{M}^{(k)})}^2\\ \nonumber
    & \text{s.t. } \quad \bm{\alpha}(\bm{x}_i,f_{\bm{\theta}^{(k)}}(\bm{x}_i)) = (\bm{x}_i,f_{\bm{\theta}}(\bm{x}_i)),\quad \forall i= 1,\ldots, N
  \end{align}
\item Update $\mathcal{M}^{(k+1)}$:
  \begin{align}
    \label{eq:update_M}
    \mathcal{M}^{(k+1)} = \bm{\alpha}^{(k+1)}(\mathcal{M}^{(k)})
  \end{align}
\end{itemize}
\begin{remark}\normalfont
As mentioned in Theorem \ref{thm:dim}, $\bm{\alpha}=(\alpha_1,\cdots,\alpha_d)$ is supposed to be the coordinate functions. In (\ref{eq:update_two}), we solve a perturbed version $\bm{\alpha}^{(k+1)}$ so that it maps the previous iterates of point cloud $\{(\bm{x}_i,f_{\bm{\theta}^{(k)}}(\bm{x}_i))\}_{i=1}^N$ and  manifold $\mathcal{M}^{(k)}$ to their corresponding updated versions. If the iteration converges to a fixed point, the consecutive iterates of the manifolds $\mathcal{M}^{(k+1)}$ and $\mathcal{M}^{(k)}$ will be very close to each other for sufficiently large $k$, and $\bm{\alpha}^{(k+1)}$ will be very close to the coordinate functions.
\end{remark}

Note that (\ref{eq:update_M}) is straightforward to implement, and (\ref{eq:update_two}) is an optimization problem with linear constraint, which can be solved via the alternating direction method of multipliers (ADMM). More specifically,
\begin{align}
  \nonumber
  & \bm{\alpha_{\bm{\xi}}}^{(k+1)} = \arg \min_{\bm{\alpha_{\bm{\xi}}}}  \sum_{j=d_1+1}^d\|\nabla_{\mathcal{M}^{(k)}}\alpha_j\|_{L^2(\mathcal{M}^{(k)})}\\    \label{eq:update_alpha}
 &+ \frac{\mu|\mathcal{M}^{(k)}|}{2\lambda N}\sum_{i=1}^N\|\bm{\alpha}_{\bm{\xi}}(\bm{x}_i,f_{\bm{\theta}^{(k)}}(\bm{x}_i))-(f_{\bm{\theta}^{(k)}}(\bm{x}_i)-Z_i^{(k)})\|^2_2.
\end{align}
\begin{align}
  \nonumber
  \bm{\theta}^{(k+1)} = &\arg \min_{\bm{\theta}} J(\bm{\theta})+ \frac{\mu}{2N}\sum_{i=1}^N\|\bm{\alpha}_{\bm{\xi}}^{(k+1)}(\bm{x}_i,f_{\bm{\theta}^{(k)}}(\bm{x}_i))\\ \label{eq:update_theta}
&-(f_{\bm{\theta}}(\bm{x}_i)-Z_i^{(k)})\|^2_2.
\end{align}
\begin{align}
  \label{eq:update_z}
  Z_i^{(k+1)} = Z_i^{(k)} + \bm{\alpha}_{\bm{\xi}}^{(k+1)}(\bm{x}_i,f_{\bm{\theta}^{(k)}}(\bm{x}_i))-f_{\bm{\theta}^{(k+1)}}(\bm{x}_i),
\end{align}
where $\bm{\alpha_\bm{\xi}}$ is defined as $\bm{\alpha} = (\bm{\alpha_{\bm{x}}},\bm{\alpha_{\bm{\xi}}}) = \left((\alpha_1,\ldots,\alpha_{d_1}),(\alpha_{d_1+1},\ldots,\alpha_{d})\right)$, and $Z_i$ is the dual variable. Note that we  need to perturb only the coordinate functions $\bm{\alpha_{\bm{\xi}}}$ corresponding to the features in (\ref{eq:update_alpha}) because the inputs $\bm{x}_i$ are given and fixed. Also note that because the gradient and the $L^2$ norm in (\ref{eq:update_alpha}) are defined on $\mathcal{M}$ instead of the projected manifold, we are \textit{not} simply minimizing the dimension of the manifold $\mathcal{M}$ projected onto the feature space. For computational efficiency, we update $\bm{\alpha},\bm{\theta}$ and $Z_i$ only once every manifold update (\ref{eq:update_M}).

Among (\ref{eq:update_alpha}),(\ref{eq:update_theta}) and (\ref{eq:update_z}), (\ref{eq:update_z}) is the easiest to implement, (\ref{eq:update_theta}) can be solved by stochastic gradient descent (SGD) with modified back propagation, and (\ref{eq:update_alpha}) can be solved by the point integral method (PIM) \cite{pim}. The detailed implementation of (\ref{eq:update_alpha}) and (\ref{eq:update_theta}) will be explained in the next section.

\section{Implementation Details and Complexity Analysis}

In this section, we present the details of the algorithmic implementation, which includes  back propagation for the $\bm{\theta}$ update (\ref{eq:update_theta}), point integral method for the $\bm{\alpha}$ update (\ref{eq:update_alpha}), and the complexity analysis.
\subsection{Back Propagation for the $\bm{\theta}$ Update}
We derive the gradient of the objective function in (\ref{eq:update_theta}). Let
\begin{align}
  \label{eq:second_term}
  E_i(\bm{\theta})=\frac{\mu}{2}\|\bm{\alpha}_{\bm{\xi}}^{(k+1)}(\bm{x}_i,f_{\bm{\theta}^{(k)}}(\bm{x}_i))-(f_{\bm{\theta}}(\bm{x}_i)-Z_i^{(k)})\|^2_2.
\end{align}
Then the objective function in (\ref{eq:update_theta}) is
\begin{align}
  \label{eq:loss_augmented}
  \tilde{J}(\bm{\theta}) = \frac{1}{N}\sum_{i=1}^N\ell(f_{\bm{\theta}}(\bm{x}_i),y_i)+\frac{1}{N}\sum_{i=1}^NE_i(\bm{\theta}).
\end{align}
Usually the back-propagation of the first term in (\ref{eq:update_theta}) is known for a given network. As for the second term, let $\bm{x}_i$ be a given datum in a mini-batch. The gradient of the second term with respect to the output layer $f_{\bm{\theta}}(\bm{x}_i)$ is:
\begin{align}
  \label{eq:gradient_E}
  \frac{\partial E_i}{\partial f_{\bm{\theta}}(\bm{x}_i)} = \mu\left(f_{\bm{\theta}}(\bm{x}_i)-Z_i^{(k)}-\bm{\alpha_{\bm{\xi}}}^{(k+1)}(\bm{x}_i,f_{\bm{\theta}^{(k)}}(\bm{x}_i))\right)
\end{align}
This means that we need to only add the extra term (\ref{eq:gradient_E}) to the original gradient, and then use the already known procedure to back-propagate the gradient. This essentially leads to no extra computational cost in the SGD updates.

\subsection{Point Integral Method for the $\bm{\alpha}$ Update}

Note that the objective funtion in (\ref{eq:update_alpha}) is decoupled with respect to $j$, and each $\alpha_j$ update can be cast into:
\begin{align}
  \label{eq:cannonical}
  \min_{u \in H^1(\mathcal{M})} \|\nabla_\mathcal{M} u\|_{L^2(\mathcal{M})}^2 + \gamma \sum_{\bm{q} \in P} |u(\bm{q}) - v(\bm{q})|^2,
\end{align}
where $u = \alpha_j$, $\mathcal{M} = \mathcal{M}^{(k)}$, $\gamma = \mu |\mathcal{M}^{(k)}|/2\lambda N$, and $P = \left\{\bm{p}_i=\left(\bm{x}_i,f_{\bm{\theta}^{(k)}}(\bm{x}_i)\right)\right\}_{i=1}^N\subset \mathcal{M}$. The Euler-Lagrange equation of (\ref{eq:cannonical}) is:
\begin{align}\nonumber
  -\Delta_\mathcal{M} u(\bm{p}) + \gamma \sum_{\bm{q}\in P}\delta(\bm{p}-\bm{q})(u(\bm{q})-v(\bm{q}))&= 0, \text{ } \bm{p} \in \mathcal{M}\\ \label{eq:EL}
\frac{\partial u}{\partial n} &= 0, \text{ } \bm{p} \in \partial \mathcal{M}
\end{align}
It is hard to discretize the Laplace-Beltrami operator $\Delta_\mathcal{M}$ and the delta function $\delta(x,y)$ on an unstructured point cloud $P$. We instead use the point integral method to solve (\ref{eq:EL}). The key observation in PIM is the following theorem:

\begin{thm} {\normalfont\cite{pim}}
  \label{thm:local-error}
If $u\in C^3(\mathcal{M})$ is a function on $\mathcal{M}$,
then
\begin{align}\nonumber
  &\left\|\int_\mathcal{M} \Delta_{\mathcal{M}} u(\bm{q})R_t(\bm{p},\bm{q})d \bm{q}- 2\int_{\partial \mathcal{M}} \frac{\partial u(\bm{q})}{\partial n} R_t(\bm{p}, \bm{q}) d \tau_{\bm{q}}\right.\\
  & +\left.\frac{1}{t}\int_{\mathcal{M}} (u(\bm{p})-u(\bm{q}))R_t(\bm{p}, \bm{q}) d \bm{q}\right\|_{L^2(\mathcal{M})}=O(t^{1/4}).
\end{align}
where $R_t$ is the normalized heat kernel:
\begin{align}
  \label{eq:Rt}
  R_t(\bm{p},\bm{q}) = C_t\exp\left(-\frac{|\bm{p}-\bm{q}|^2}{4t}\right).
\end{align}
\end{thm}

After convolving equation (\ref{eq:EL}) with the heat kernel $R_t$, we know the solution $u$ of (\ref{eq:EL}) should satisfy
\begin{align}
  \nonumber
  &-\int_\mathcal{M} \Delta_{\mathcal{M}} u (\bm{q})R_t(\bm{p},\bm{q})d\bm{q}\\  \label{eq:EL-convolve}
  + &\gamma\sum_{\bm{q}\in P}R_t(\bm{p},\bm{q})\left(u(\bm{q})-v(\bm{q})\right) = 0.
\end{align}
Combined with Theorem \ref{thm:local-error} and the  Neumann boundary condition, this implies that $u$ should approximately satisfy
\begin{align}
 \nonumber
  &\int_\mathcal{M}\left(u(\bm{p})-u(\bm{q})\right)R_t(\bm{p},\bm{q})d\bm{q}\\
\label{eq:linear-system-continuous}
+ &\gamma t \sum_{\bm{q} \in P}R_t(\bm{p},\bm{q})\left(u(\bm{q}) - v(\bm{q})\right) = 0
\end{align}

Note that (\ref{eq:linear-system-continuous}) no longer involves the gradient $\nabla_\mathcal{M}$ or the Laplace-Beltrami operator $\Delta_\mathcal{M}$. Assume that $P = \{\bm{p}_1,\ldots,\bm{p}_N\}$ samples the manifold $\mathcal{M}$ uniformly at random, then (\ref{eq:linear-system-continuous}) can be discretized as
\begin{align}
  \label{eq:linear-system-discretize}
\frac{|\mathcal{M}|}{N}\sum_{j=1}^NR_{t,ij}(u_i-u_j) +\gamma t \sum_{j=1}^N R_{t,ij} (u_j - v_j) = 0,
\end{align}
where $u_i= u(\bm{p}_i)$, and $R_{t,ij}= R_t(\bm{p}_i,\bm{p}_j)$. Combining the definition of $\gamma$ in (\ref{eq:cannonical}), we can write (\ref{eq:linear-system-discretize}) in the matrix form
\begin{align}
  \label{eq:linear-system-matrix}
\left(\bm{L}+  \frac{\mu}{\tilde{\lambda}} \bm{W}\right) \bm{u} = \frac{\mu}{\tilde{\lambda}} \bm{W} \bm{v}, \quad \tilde{\lambda} = 2\lambda/t,
\end{align}
where $\tilde{\lambda}$ can be chosen instead of $\lambda$ as the hyperparameter to be tuned, $\bm{u} = (u_1, \ldots, u_N)^T$, $\bm{W}$ is an $N\times N$ matrix
\begin{align}
  \label{eq:W}
  \bm{W}_{ij} = R_{t,ij} = \exp\left(-\frac{|\bm{p}_i-\bm{p}_j|^2}{4t}\right),
\end{align}
and $\bm{L}$ is the graph Laplacian of $\bm{W}$:
\begin{align}
  \label{eq:L}
  \bm{L}_{ii} = \sum_{j\not=i}\bm{W}_{ij}, \quad \text{and}\quad \bm{L}_{ij} = -\bm{W}_{ij}\quad \text{if} \quad i\not=j.
\end{align}
Therefore, the update of $\bm{\alpha_\bm{\xi}}$, which is cast into the canonical form (\ref{eq:cannonical}), is achieved by solving a linear system (\ref{eq:linear-system-matrix}).

\subsection{Complexity Analysis}
\label{sec:complexity}
Based on the analysis above, we present a summary of the traning for \textit{LDMNet} in Algorithm \ref{alg:ldmnet}.
\begin{algorithm}
\floatname{algorithm}{Algorithm}
\caption{\textit{LDMNet} Training}
\label{alg:ldmnet}
\begin{algorithmic
\REQUIRE Training data $\{(\bm{x}_i,y_i)\}_{i=1}^N \subset \mathbb{R}^{d_1}\times \mathbb{R}$, hyperparameters $\tilde{\lambda}$ and $\mu$, and a neural network with the weights $\bm{\theta}$ and the output layer $\bm{\xi}_i = f_{\bm{\theta}}(\bm{x}_i)\in \mathbb{R}^{d_2}$.
\ENSURE  Trained network weights $\bm{\theta}^*$.
\STATE Randomly initialize the network weights $\bm{\theta}^{(0)}$. The dual variables $Z_i^{(0)} \in \mathbb{R}^{d_2}$ are initialized to zero.
\WHILE {not converge}
\STATE 1. Compute the matrices $\bm{W}$ and $\bm{L}$ as in (\ref{eq:W}) and (\ref{eq:L}) with $\bm{p}_i = (\bm{x}_i,f_{\bm{\theta}^{(k)}}(\bm{x}_i))$.
\vspace{1mm}
\STATE 2. Update $\bm{\alpha}^{(k+1)}$ in (\ref{eq:update_alpha}): solve the linear systems (\ref{eq:linear-system-matrix}), where
\vspace{-3mm}
\begin{align}
  \label{eq:u_and_v}
  \bm{u}_i=\alpha_j(\bm{p}_i),\quad \bm{v}_i = f_{\bm{\theta}^{(k)}}(\bm{x}_i)_j-Z_{i,j}^{(k)}.
\end{align}
\vspace{-5mm}
\STATE 3. Update $\bm{\theta}^{(k+1)}$ in (\ref{eq:update_theta}): run SGD for $M$ epochs with an extra gradient term (\ref{eq:gradient_E}).
\STATE 4. Update $Z^{(k+1)}$ in (\ref{eq:update_z}).
\STATE 5. $k\leftarrow k+1$.
\ENDWHILE
\STATE $\bm{\theta}^*\leftarrow \bm{\theta}^{(k)}$.
\end{algorithmic}
\end{algorithm}

The additional computation in Algorithm \ref{alg:ldmnet} (in steps 1 and 2) comes from the update of weight matrices in (\ref{eq:W}) and solving the linear system (\ref{eq:linear-system-matrix}) from PIM once every $M$ epochs of SGD. We now explain the computational complexity of these two steps.

When $N$ is large, it is not computationally feasible to compute the pairwise distances in the entire training set. Therefore the weight matrix $\bm{W}$ is truncated to only $20$ nearest neighbors. To identify those nearest neighbors, we first organize the data points $\{\bm{p}_1,\ldots,\bm{p}_N\}\subset \mathbb{R}^d$ into a $k$-d tree \cite{kdtree}, which is a binary tree that recursively partitions a $k$-dimensional space (in our case $k = d$). Nearest neighbors can then be efficiently identified because branches can be eliminated from the search space quickly. Modern algorithms to build a balanced $k$-d tree generally at worst converge in $O(N\log N)$ time \cite{kdtree_build_1,kdtree_build_2}, and finding nearest neighbours for one query point in a balanced $k$-d tree takes $O(\log N)$ time on average \cite{kdtree_search}. Therefore the complexity of the weight update is $O(N\log N)$.

Since $\bm{W}$ and $\bm{L}$ are sparse symmetric matrices with a fixed maximum number of non-zero entries in each row, the linear system (\ref{eq:linear-system-matrix}) can be solved efficiently with the preconditioned conjugate gradients method. After restricting the number of matrix multiplications to a maximum of 50, the complexity of the $\bm{\alpha}$ update is $O(N)$.

\section{Experiments}
\label{sec:experiments}
In this section, we compare the performance of \textit{LDMNet} to widely-used network regularization techniques, weight decay and \textit{DropOut}, using the same underlying network structure. We point out that our focus is to compare the effectiveness of the regularizers, and not to investigate the state-of-the-art performances on the benchmark datasets. Therefore we typically use simple network structures and relatively small training sets, and no data augmentation or early stopping is implemented.

Unless otherwise stated, all experiments use mini-batch SGD with momentum on batches of 100 images. The momentum parameter is fixed at $0.9$. The networks are trained using a fixed learning rate $r_0$ on the first $200$ epochs, and then $r_0/10$ for another $100$ epochs.

As mentioned in Section \ref{sec:complexity}, the weight matrices $\bm{W}$ are truncated to $20$ nearest neighbors. For classification tasks, nearest neighbors can be searched within each class in the labeled training set. We also normalize the weight matrices with local scaling factors $\sigma(\bm{p})$ \cite{zelnik2005self}:
\begin{align}
  \label{eq:weight}
  w(\bm{p},\bm{q}) = \exp \left(-\frac{\|\bm{p}-\bm{q}\|^2}{\sigma(\bm{p})\sigma(\bm{q})}\right),
\end{align}
where $\sigma(\bm{p})$ is chosen as the distance between $\bm{p}$ and its $10$th nearest neighbor. This is based on the empirical analysis on choosing the parameter $t$ in \cite{pim}, and has been used in \cite{ldmm_scientific}. The weight matrices and $\bm{\alpha}$ are updated once every $M=2$ epochs of SGD.

All hyperparameters are optimized so that those reported are the best performance of each method. For \textit{LDMNet}, $\tilde{\lambda}$ defined in (\ref{eq:linear-system-matrix}) typically decreases as the training set becomes larger, whereas the paramter $\mu$ for augmented Lagrangian can be fixed to be a constant. For weight decay,
\begin{align}
  \label{eq:weight_decay}
  \min_{\bm{\theta}} J(\bm{\theta}) + w\|\bm{\theta}\|_2^2,
\end{align}
the parameter $w$ also usually decreases as the training size increases. For \textit{DropOut}, the corresponding \textit{DropOut} layer is always chosen to have a drop rate of $0.5$.

\subsection{MNIST}

The MNIST handwritten digit dataset contains approximately 60,000 training images ($28\times 28$) and 10,000 test images. Tabel \ref{tab:mnist_net} describes the network structure. The learned feature $f_{\bm{\theta}}(\bm{x}_i)$ for the training data $\bm{x}_i$ is the output of layer six, and is regularized for \textit{LDMNet} in (\ref{eq:ldm}).

\begin{table}
  \centering
  \begin{tabular}{c|c|C{4cm}}
    \hline\hline
    Layer& Type & Parameters\\
    \hline
    1 & conv & size: $5\times 5\times 1\times 20$ \newline stride: 1, pad: 0\\
    \hline
    2 & max pool & size: $2\times 2$, stride: 2, pad: 0\\
    \hline
    3 & conv & size: $5\times 5\times 20\times 50$ \newline stride: 1, pad: 0\\
    \hline
    4 & max pool & size: $2\times 2$, stride: 2, pad: 0\\
    \hline
    5 & conv & size: $4\times 4\times 50\times 500$ \newline stride: 1, pad: 0\\
    \hline
    6 & ReLu (\textit{DropOut}) & N/A\\
    \hline
    7 & fully connected & $500 \times 10$\\
    \hline
    8 & softmaxloss & N/A\\
    \hline\hline
  \end{tabular}
  \vspace{1mm}
  \caption{Network structure in the MNIST experiments. The outputs of layer 6 are the extracted features, which will be fed into the softmax classifier (layer 7 and 8).}
  \label{tab:mnist_net}
\end{table}

\begin{table}
  \centering
  \begin{tabular}{C{2cm}|C{1.5cm}|C{1.5cm}|C{1.5cm}}
    \hline\hline
    training per class & $\tilde{\lambda}$ & $\mu$ & $w$ \\
    \hline
    50 & 0.05 & 0.01 & 0.1\\
    \hline
    100 & 0.05 & 0.01 & 0.05\\
    \hline
    400 & 0.01 & 0.01 & 0.01\\
    \hline
    700 & 0.01 & 0.01 & 0.005\\
    \hline
    1000 & 0.005 & 0.01 & 0.005\\
    \hline
    3000 & 0.001 & 0.01 & 0.001\\
    \hline
    6000 & 0.001 & 0.01 & 0.001\\
    \hline\hline
  \end{tabular}
  \vspace{1mm}
  \caption{Hyperparamters used in the MNIST experiments}
  \label{tab:mnist_hyper}
\end{table}

\begin{table}
  \centering
  \begin{tabular}{C{2cm}|C{1.5cm}|C{1.5cm}|C{1.5cm}}
    \hline\hline
    training per class & weight decay & \textit{DropOut} & \textit{LDMNet} \\
    \hline
    50 & 91.32\% & 92.31\% & \textbf{95.57\%}\\
    \hline
    100 & 93.38\% & 94.05\% & \textbf{96.73\%}\\
    \hline
    400 & 97.23\% & 97.95\% &\textbf{98.41\%} \\
    \hline
    700 & 97.67\% & 98.07\% & \textbf{98.61\%}\\
    \hline
    1000 &98.06\%  & 98.71\% & \textbf{98.89\%}\\
    \hline
    3000 & 98.87\% & 99.21\% & \textbf{99.24\%}\\
    \hline
    6000 & 99.15\% & \textbf{99.41\%} & 99.39\%\\
    \hline\hline
  \end{tabular}
  \vspace{1mm}
  \caption{MNIST: testing accuracy for different regularizers}
  \label{tab:mnist_result}
\end{table}

\begin{figure}[t]
  \centering
  \begin{tabular}{cc}
    \includegraphics[width=0.48\linewidth]{mnist_class_error}&
    \includegraphics[width=0.48\linewidth]{mnist_gen_error}\\
  \end{tabular}
\caption{Comparison of the regularizers on the MNIST dataset. The first (second) figure shows the dependence of the classification (generalization) error on the size of the training set.}
\label{fig:mnist_error}
\end{figure}

While state-of-the-art methods often use the entire training set, we are interested in examining the performance of the regularization techniques with varying training sizes from 500 to 60,000. In this experiment, the initial learning rate is set to $0.001$, and the hyperparameters are reported in Table \ref{tab:mnist_hyper}. Table \ref{tab:mnist_result} displays the testing accuracy of the competing algorithms. The dependence of the classification error and generalization error (which is the difference between the softmax loss on the testing and training data) on the size of the training set is shown in Figure~\ref{fig:mnist_error}. Figure~\ref{fig:visual_mnist} provides a visual illustration of the features of the testing data learned from 1,000 training samples. It is clear to see that \textit{LDMNet} significantly outperforms weight decay and \textit{DropOut} when the training set is small, and the performance becomes broadly similar as the size of the training set reaches 60,000.

\subsection{SVHN and CIFAR-10}
\begin{table}
  \centering
  \begin{tabular}{c|c|C{4cm}}
    \hline\hline
    Layer& Type & Parameters\\
    \hline
    1 & conv & size: $5\times 5\times 3\times 96$ \newline stride: 1, pad: 2\\
    \hline
    2 & ReLu & N/A\\
    \hline
    3 & max pool & size: $3\times 3$, stride: 2, pad: 0\\
    \hline
    4 & conv & size: $5\times 5\times 96\times 128$ \newline stride: 1, pad: 2\\
    \hline
    5 & ReLu & N/A\\
    \hline
    6 & max pool & size: $3\times 3$, stride: 2, pad: 0\\
    \hline
    7 & conv & size: $4\times 4\times 128\times 256$ \newline stride: 1, pad: 0\\
    \hline
    8 & ReLu & N/A\\
    \hline
    9 & max pool & size: $3\times 3$, stride: 2, pad: 0\\
    \hline
    10 & fully connected & output: 2048\\
    \hline
    11 & ReLu (\textit{DropOut}) & N/A \\
    \hline
    12 & fully connected & output: 2048\\
    \hline
    13 & ReLu (\textit{DropOut}) & N/A\\
    \hline
    14 & fully connected & $2048\times 10$\\
    \hline
    15 & softmaxloss & N/A\\
    \hline\hline
  \end{tabular}
  \vspace{1mm}
  \caption{Network structure in the SVHN and CIFAR-10 experiments. The outputs of layer 13 are the extracted features, which will be fed into the softmax classifier (layer 14 and 15).}
  \label{tab:svhn_cifar_net}
\end{table}

\begin{table}[!t]
  \centering
  \begin{tabular}{C{1.5cm}|c|C{1.5cm}|c|C{1.5cm}}
    \hline\hline
    training& \multicolumn{2}{c|}{SVHN}&\multicolumn{2}{c}{CIFAR-10}\\
    \cline{2-5}per class & $\tilde{\lambda}$ & $w$ &  $\tilde{\lambda}$ & $w$\\ \hline
    50  & 0.1  & $10^{-6}$ & 0.01 & $5\times 10^{-4}$\\ \hline
    100 & 0.05 & $10^{-6}$ & 0.01 & $5\times 10^{-5}$\\ \hline
    400 & 0.05 & $10^{-7}$ & 0.01 & $5\times 10^{-5}$\\ \hline
    700 & 0.01 & $10^{-8}$ & 0.01 & $5\times 10^{-7}$\\ \hline\hline
  \end{tabular}
  \vspace{1mm}
  \caption{Hyperparameters used in the SVHN and CIFAR-10 experiments.}
  \label{tab:svhn_cifar_hyper}
\end{table}

\begin{table}
  \centering
  \begin{tabular}{C{2cm}|C{1.5cm}|C{1.5cm}|C{1.5cm}}
    \hline\hline
    training per class & weight decay & \textit{DropOut} & \textit{LDMNet} \\
    \hline
    50 & 71.46\% & 71.94\% & \textbf{74.64\%}\\
    \hline
    100 & 79.05\% & 79.94\% & \textbf{81.36\%}\\
    \hline
    400 & 87.38\% & 87.16\% &\textbf{88.03\%} \\
    \hline
    700 & 89.69\% & 89.83\% & \textbf{90.07\%}\\
    \hline\hline
  \end{tabular}
  \vspace{1mm}
  \caption{SVHN: testing accuracy for different regularizers}
  \label{tab:svhn_result}
\end{table}

\begin{table}
  \centering
  \begin{tabular}{C{2cm}|C{1.5cm}|C{1.5cm}|C{1.5cm}}
    \hline\hline
    training per class & weight decay & \textit{DropOut} & \textit{LDMNet} \\
    \hline
    50 & 34.70\% & 35.94\% & \textbf{41.55\%}\\
    \hline
    100 & 42.45\% & 43.18\% & \textbf{48.73\%}\\
    \hline
    400 & 56.19\% & 56.79\% &\textbf{60.08\%} \\
    \hline
    700 & 61.84\% & 62.59\% & \textbf{65.59\%}\\
    \hline
    full data & \multicolumn{2}{c|}{87.72\%} & \textbf{88.21\%}\\
    \hline\hline
  \end{tabular}
  \vspace{1mm}
  \caption{CIFAR-10: testing accuracy for different regularizers. The first four experiments use training sets varied in size from 500 to 7,000 and a simple network specified in Table~\ref{tab:svhn_cifar_net}. In the last experiment, we first train a DNN with VGG-16 architecture \cite{vgg16} from the full training data using both weight decay and \textit{DropOut}, then the DNN is fine-tuned by regularizing the output layer with \textit{LDMNet}. No data augmentation is implemented in any of the experiment.}
  \label{tab:cifar_result}
\end{table}

SVHN and CIFAR-10 are benchmark RGB image datasets, both of which contain 10 different classes. These two datasets are  more challenging than the MNIST dataset because of their weaker intraclass correlation. All algorithms use the network structure (similar to \cite{hinton2012improving}) specified in Table~\ref{tab:svhn_cifar_net}. The outputs of layer 13 are the learned features, and are regularized in \textit{LDMNet}. All algorithms start with a learning rate of $0.005$ for SVHN and $0.001$ for CIFAR-10. $\mu$ has been fixed as $0.5$ for SVHN and $1$ for CIFAR-10, and the remaining hyperparameters are reported in Table~\ref{tab:svhn_cifar_hyper}.

We report the testing accuracies of the competing regularizers in Table~\ref{tab:svhn_result} and Table~\ref{tab:cifar_result} when the number of training samples is varied from $50$ to $700$ per class. Again, it is clear to see that \textit{LDMNet} outperforms weight decay and \textit{DropOut} by a significant margin.

To demonstrate the generality of \textit{LDMNet}, we conduct another experiment on CIFAR-10 using the entire training data with a different network structure. We first train a DNN with VGG-16 architecture \cite{vgg16} on CIFAR-10 using both weight decay and \textit{DropOut} without data augmentation. Then, we fine-tune the DNN by regularizing the output layer with \textit{LDMNet}. The testing accuracies are reported in the last row of Table~\ref{tab:cifar_result}. Again, \textit{LDMNet} outperforms weight decay and \textit{DropOut}, demonstrating that \textit{LDMNet} is a general framework that can be used to improve the performance of any network structure.

\subsection{NIR-VIS Heterogeneous Face Recognition}

\begin{figure}[t]
  \centering
  \begin{tabular}{cccc}
    \includegraphics[width=0.24\linewidth]{nir_1_1}&\hspace{-4mm}
    \includegraphics[width=0.24\linewidth]{nir_1_2}&\hspace{-4mm}
    \includegraphics[width=0.24\linewidth]{nir_2_1}&\hspace{-4mm}
    \includegraphics[width=0.24\linewidth]{nir_2_2}\\
    \includegraphics[width=0.24\linewidth]{vis_1_1}&\hspace{-4mm}
    \includegraphics[width=0.24\linewidth]{vis_1_2}&\hspace{-4mm}
    \includegraphics[width=0.24\linewidth]{vis_2_1}&\hspace{-4mm}
    \includegraphics[width=0.24\linewidth]{vis_2_2}
  \end{tabular}
\caption{Sample images of two subjects from the CASIA NIR-VIS 2.0 dataset after the pre-procssing of alignment and cropping \cite{Kazemi_2014_CVPR}. Top: NIR. Bottom: VIS.}
\label{fig:nirvis_sample}
\end{figure}

Finally, we demonstrate the effectiveness of \textit{LDMNet} for NIR-VIS face recognition. The objective of the experiment is to match a probe image of a subject captured in the near-infrared spectrum (NIR) to the same subject from a gallery of visible spectrum (VIS) images. The CASIA NIR-VIS 2.0 benchmark dataset \cite{nirvis} is used to evaluate the performance. This dataset contains 17,580 NIR and VIS face images of 725 subjects. Figure~\ref{fig:nirvis_sample} shows eight sample images of two subjects after facial landmark alignment and cropping \cite{Kazemi_2014_CVPR}. Despite recent breakthroughs for VIS face recognition by training DNNs from millions of VIS images, such approach cannot be simply transferred to NIR-VIS face recognition. The reason is that, unlike VIS face images, we have only limited number of availabe NIR images. Moreover, the NIR-VIS face matching is a cross-modality comparison.

The authors in \cite{lezama2016not} introduced a way to transfer the breakthrough in VIS face recognition to the NIR spectrum. Their idea is to use a DNN pre-trained on VIS images as a feature extactor, while making two independent modifications in the input and output of  the DNN. They first modify the input by ``hallucinating'' a VIS image from the NIR sample, and then apply a low-rank embedding of the DNN features at the output. The combination of these two modifications achieves the state-of-the-art performance on cross-spectral face recognition.

\begin{table}
  \centering
  \begin{tabular}{c|c|C{3cm}}
    \hline\hline
    Layer& Type & Parameters\\
    \hline
    1 & fully connected & output:2000\\
    \hline
    2 & ReLu (\textit{DropOut})& N/A\\
    \hline
    3 & fully connected & output:2000\\
    \hline
    4 & ReLu (\textit{DropOut})& N/A\\
    \hline\hline
  \end{tabular}
  \vspace{1mm}
  \caption{Fully connected network for the NIR-VIS nonlinear feature embedding. The outputs of layer 4 are the extracted features.}
  \label{tab:nirvis_net}
\end{table}

We follow the second idea in \cite{lezama2016not}, and learn a nonlinear low dimensional manifold embedding of the output features. The intuition is that faces of the same subject in two different modalities should sample the same low dimensional feature manifold in a transformed space. In our experiment, we use the VGG-face model (downloaded at \url{http://www.robots.ox.ac.uk/~vgg/software/vgg_face/}) \cite{parkhi2015deep} as a feature extractor. The learned 4,096 dimensional features can be reduced to a 2,000 dimensional space using PCA and used directly for face matching. Meanwhile, we also put the 4,096 dimensional features into a two-layer fully connected network described in Table~\ref{tab:nirvis_net} to learn a nonlinear embedding using different regularizations. The features extracted from layer 4 are regularized in \textit{LDMNet}.

All nonlinear embeddings using the structure specified in Table~\ref{tab:nirvis_net} are trained with SGD on mini-batches of 100 images for 200 epochs. We use an exponentially decreasing  learning rate that starts at 0.1 with a decaying factor of 0.99. The hyperparameters are chosen to achieve the optimal performance on the validation set. More specifically, $\tilde{\lambda}, \mu, w$ are set to $5\times 10^{-5},5$, and $5\times 10^{-4}$ respectively.

\begin{table}
\vspace{1mm}
  \centering
  \begin{tabular}{lc}
    \hline
     & Accuracy (\%)\\
    \hline\hline
    VGG-face & $74.51 \pm 1.28$\\
    VGG-face + triplet \cite{lezama2016not}& $75.96 \pm 2.90$\\
    VGG-face + low-rank \cite{lezama2016not}& $80.69 \pm 1.02$\\
    VGG-face weight Decay & $63.87 \pm 1.33$\\
    VGG-face \textit{DropOut} & $66.97 \pm 1.31$\\
    VGG-face \textit{LDMNet} & \textbf{85.02} $\bm{\pm}$ \textbf{0.86}\\
    \hline\hline
  \end{tabular}
  \vspace{1mm}
  \caption{NIR-VIS cross-spectral rank-1 identification rate on the 10-fold CASIA NIR-VIS 2.0 benchmark. The first result is obtained by reducing the features learned from VGG-face to 2,000 dimensional space using PCA. The next two results use triplet \cite{weinberger2006distance} and low-rank embedding of the learned features, and are reported in \cite{lezama2016not}. The last three results are achieved by training the nonlinear embedding network in Table~\ref{tab:nirvis_net} with the corresponding regularizations.}
  \label{tab:nirvis_result}
\end{table}

We report the rank-1 performance score for the standard CASIA NIR-VIS 2.0 evaluation protocol in Table~\ref{tab:nirvis_result}. Because of the limited amount of training data (around 6,300 NIR and 2,500 VIS images), the fully-connected networks in Table~\ref{tab:nirvis_net} trained with weight decay and \textit{DropOut} clearly overfit the training data: they actually yield testing accuracies that are worse than using a simple PCA embedding of the features learned from VGG-face. However, the same network regularized with \textit{LDMNet} has achieved a significant 10.5\% accuracy boost (from 74.51\% to 85.02\%) to using VGG-face directly. It is also better than the results reported in \cite{lezama2016not} using the popular triplet embedding \cite{weinberger2006distance} and low-rank embedding. Figure~\ref{fig:visual} provides a visual illustration of the learned features from different regularizations. The generated features of five subjects are visualized in two dimensions using PCA, with filled circle for VIS, and unfilled diamond for NIR, and one color for each subject. Note that in Figure~\ref{fig:vgg},\ref{fig:wd},\ref{fig:dropout}, features of one subject learned directly from VGG-face or from a nonlinear embedding regularized with weight decay or \textit{DropOut} typically form two clusters, one for NIR and the other for VIS. This contrasts with the behavior in Figure~\ref{fig:ldm}, where features of the same subject from two different modalities merge to form a single low dimensional manifold.

\section{Conclusion}
We proposed a general deep neural network regularization technique \textit{LDMNet}. The intuition of \textit{LDMNet} is that the concatenation of the input data and output features should sample a collection of low dimensional manifolds, an idea that is loosely related to the Gaussian mixture model.  Unlike most data-dependent regularizations, \textit{LDMNet} focuses on the geometry of both the input data and the output features, and does not require explicit parametrization of  the underlying manifold. The dimensions of the manifolds are directly regularized in a variational form, which is solved by alternating direction of minimization with a slight increase in the computational complexity ($O(N)$ in solving the linear system, and $O(N\log N)$ in the weight update). Extensive experiments show that \textit{LDMNet} has two benefits: First, it significantly outperforms widely used regularization techniques, such as weight decay and \textit{DropOut}. Second, \textit{LDMNet} can extract common features  of the same subject imaged via different modalities so that these features sample the same low dimensional manifold, which significantly increases the accuracy of cross-spectral face recognition.

{\small
\bibliographystyle{abbrv}


\section{Introduction}
Human visual attention enables drivers to quickly identify and locate potential risks or important visual cues across the visual field, such as a darting-out pedestrian, an incursion of a nearby cyclist or a changing traffic light. Drivers' gaze behavior has been studied as a proxy for their attention. Recently, a large driver attention dataset of routine driving \cite{dreyeve2016} has been introduced and neural networks \cite{palazzi2017learning,tawari2017computational} have been trained end-to-end to estimate driver attention, mostly in lane-following and car-following situations. Nonetheless, datasets and prediction models for driver attention in rare and critical situations are still needed.

However, it is nearly impossible to collect enough driver attention data for crucial events with the conventional in-car data collection protocol, \textit{i}.\textit{e}.,~ collecting eye movements from drivers during driving. This is because the vast majority of routine driving situations consist of simple lane-following and car-following. In addition, collecting driver attention in-car has two other major drawbacks. (i) Single focus: at each moment the eye-tracker can only record one location that the driver is looking at, while the driver may be attending to multiple important objects in the scene with their covert attention, \textit{i}.\textit{e}.,~ the ability to fixate one's eyes on one object while attending to another object \cite{cavanagh2005tracking}. (ii) False positive gazes: human drivers also show eye movements to driving-irrelevant regions, such as sky, trees, and buildings \cite{palazzi2017learning}. It is challenging to separate these false positives from gazes that are dedicated to driving.

An alternative that could potentially address these concerns is showing selected driving videos to drivers in the lab and collecting their eye movements with repeated measurements while they perform a proper simulated driving task. Although this third-person driver attention collected in the lab is inevitably different from the first-person driver attention in the car, it can still potentially reveal the regions a driver should look at in that particular driving situation from a third-person perspective. These data are greatly valuable for identifying risks and driving-relevant visual cues from driving scenes. To date, a proper data collection protocol of this kind is still missing and needs to be formally introduced and tested.

Another challenge for driver attention prediction, as well as for other driving-related machine learning problems, is that the actual cost of making a particular prediction error is unknown. Attentional lapses while driving on an empty road does not cost the same as attentional lapses when a pedestrian darts out. Since current machine learning algorithms commonly rely on minimizing average prediction error, the critical moments, where the cost of making an error is high, need to be properly identified and weighted.

\begin{figure}[t]
\centering
\includegraphics[width=\linewidth]{figures/intro3.eps}
\caption{An example of input raw images ({\it left}), ground-truth human attention maps collected by us ({\it middle}), and the attention maps predicted by our model ({\it right}). The driver had to sharply stop at the green light to avoid hitting two pedestrians running the red light. The collected human attention map accurately shows the multiple regions that simultaneously demand the driver's attention. Our model correctly attends to the crossing pedestrians and does not give false alarms to other irrelevant pedestrians}
\label{fig:intro}
\end{figure}

Here, our paper offers the following novel contributions. First, in order to overcome the drawbacks of the conventional in-car driver attention collection protocol, we introduce a new protocol that uses crowd-sourced driving videos containing interesting events and makes multi-focus driver attention maps by averaging gazes collected from multiple human observers in lab with great accuracy (Fig.~\ref{fig:intro}). We will refer to this protocol as the in-lab driver attention collection protocol. We show that data collected with our protocol reliably reveal where a experienced driver should look and can serve as a substitute for data collected with the in-car protocol. We use our protocol to collect a large driver attention dataset of braking events, which is, to the best of our knowledge, the richest to-date in terms of the number of interactions with other road agents. We call this dataset Berkeley DeepDrive Attention (BDD-A) dataset and will make it publicly available. Second, we introduce Human Weighted Sampling (HWS), which uses human driver eye movements to identify which frames in the dataset are more crucial driving moments and weights the frames according to their importance levels during model training. We show that HWS improve model performance on both the entire testing set and the subset of crucial frames. Third, we propose a new driver attention prediction model trained on our dataset with HWS. The model shows sophisticated behaviors such as picking out pedestrians suddenly crossing the road without being distracted by the pedestrians safely walking in the same direction as the car (Fig.~\ref{fig:intro}). The model prediction is nearly indistinguishable from ground-truth based on human judges, and it also matches the state-of-the-art performance level when tested on an existing in-car driver attention dataset collected during driving.

\section{Related works}

\myparagraph{Image / Video Saliency Prediction: }
A large variety of the previous saliency studies explored different bottom-up feature-based models \cite{bruce2006saliency,valenti2009image,erdem2013visual,murray2011saliency,zhang2013saliency,bruce2009saliency} combining low-level features like contrast, rarity, symmetry, color, intensity and orientation, or topological structure from a scene \cite{zhang2013saliency,harel2007graph,wei2012geodesic}. Recent advances in deep learning have achieved a considerable improvement for both image saliency prediction \cite{kummerer2014deep,huang2015salicon,liu2015predicting,kummerer2016deepgaze} and video saliency prediction \cite{bazzani2016recurrent,cornia2016predicting,liu2017predicting}. These models have achieved start-of-the-art performance on visual saliency benchmarks collected mainly when human subjects were doing a free-viewing task, but models that are specifically trained for predicting the attention of drivers are still needed.

\myparagraph{Driver Attention Datasets: }
DR(eye)VE \cite{dreyeve2016} is the largest and richest existing driver attention dataset. It contains 6 hours of driving data, but the data was collected from only 74 rides, which limits the diversity of the dataset. In addition, the dataset was collected in-car and has the drawbacks we introduced earlier, including missing covert attention, false positive gaze, and limited diversity. The driver's eye movements were aggregated over a small temporal window to generate an attention map for a frame, so that multiple important regions of one scene might be annotated. But there was a trade-off between aggregation window length and gaze location accuracy, since the same object may appear in different locations in different frames. Reference \cite{fridman2016driver} is another large driver attention dataset, but only six coarse gaze regions were annotated and the exterior scene was not recorded. References \cite{simon2009alerting} and \cite{underwood2011decisions} contain accurate driver attention maps made by averaging eye movements collected from human observers in-lab with simulated driving tasks. But the stimuli were static driving scene images and the sizes of their datasets are small (40 frames and 120 frames, respectively).

\myparagraph{Driver Attention Prediction: }
Self-driving vehicle control has made notable progress in the last several years. One of major approaches is a mediated perception-based approach -- a controller depends on recognizing human-designated features, such as lane markings, pedestrians, or vehicles. Human driver’s attention provides important visual cues for driving, and thus efforts to mimic human driver’s attention have increasingly been introduced. Recently, several deep neural models have been utilized to predict where human drivers should pay attention~\cite{palazzi2017learning,tawari2017computational}. Most of existing models were trained and tested on the DR(eye)VE dataset~\cite{dreyeve2016}. While this dataset is an important contribution, it contains sparse driving activities and limited interactions with other road users. Thus it is restricted in its ability to capture diverse human attention behaviors. Models trained with this dataset tend to become vanishing point detectors, which is undesirable for modeling human attention in urban driving environment, where drivers encounter traffic lights, pedestrians, and a variety of other potential cues and obstacles. In this paper, we provide our human attention dataset as a contribution collected from a publicly available large-scale crowd-sourced driving video dataset~\cite{xu2016end}, which contains diverse driving activities and environments, including lane following, turning, switching lanes, and braking in cluttered scenes.

\section{Berkeley DeepDrive Attention (BDD-A) Dataset}
\myparagraph{Dataset Statistics: } The statistics of our dataset are summarized and compared with the largest existing dataset (DR(eye)VE) \cite{dreyeve2016} in Table~\ref{tab:dataset_statistics}. Our dataset was collected using videos selected from a publicly available, large-scale, crowd-sourced driving video dataset, BDD100k~\cite{xu2016end,yu2018bdd100k}. BDD100K contains human-demonstrated dashboard videos and time-stamped sensor measurements collected during urban driving in various weather and lighting conditions. To efficiently collect attention data for critical driving situations, we specifically selected video clips that both included braking events and took place in busy areas (see supplementary materials for technical details). We then trimmed videos to include 6.5 seconds prior to and 3.5 seconds after each braking event. It turned out that other driving actions, \textit{e}.\textit{g}.,~ turning, lane switching and accelerating, were also included. 1,232 videos (=3.5 hours) in total were collected following these procedures. Some example images from our dataset are shown in Fig. \ref{fig:examples}. Our selected videos contain a large number of different road users. We detected the objects in our videos using YOLO~\cite{redmon2017yolo9000}
On average, each video frame contained 4.4 cars and 0.3 pedestrians, multiple times more than the DR(eye)VE dataset (Table~\ref{tab:dataset_statistics}).

\begin{table}
    \caption{Comparison between driver attention datasets}
    \label{tab:dataset_statistics}
	\begin{center}
    	\resizebox{\linewidth}{!}{%
    	\begin{tabular}{@{}lccccccc@{}} \toprule
        	\multirow{2}{*}{\parbox{1.5cm}{Dataset}} & \multirow{2}{*}{\parbox{1.5cm}{\centering \# Rides}} & \multirow{1}{*}{\parbox{1.5cm}{\centering Durations}}  & \multirow{2}{*}{\parbox{1.5cm}{\centering \# Drivers}} & \multirow{2}{*}{\parbox{1.5cm}{\centering \# Gaze providers}}    & \multirow{1}{*}{\parbox{1.5cm}{\centering \# Cars}} & \multirow{1}{*}{\# Pedestrians} & \multirow{2}{*}{\parbox{2cm}{\centering \# Braking events}} \\
            & & (hours) & & & (per frame) & (per frame) & \\ \midrule
DR(eye)VE~\cite{dreyeve2016} & 74    & 6   & 8       & 8              & 1.0 & 0.04  & 464            \\
BDD-A    & 1,232 & 3.5 & 1,232   & 45             & 4.4 & 0.25  & 1,427          \\\bottomrule
        \end{tabular}}

    \end{center}
\end{table}

\myparagraph{Data Collection Procedure: } For our eye-tracking experiment, we recruited 45 participants who each had more than one year of driving experience. The participants watched the selected driving videos in the lab while performing a driving instructor task: participants were asked to imagine that they were driving instructors sitting in the copilot seat and needed to press the space key whenever they felt it necessary to correct or warn the student driver of potential dangers. Their eye movements during the task were recorded at 1000 Hz with an EyeLink 1000 desktop-mounted infrared eye tracker, used in conjunction with the Eyelink Toolbox scripts~\cite{cornelissen2002eyelink} for MATLAB. Each participant completed the task for 200 driving videos. Each driving video was viewed by at least 4 participants. The gaze patterns made by these independent participants were aggregated and smoothed to make an attention map for each frame of the stimulus video (see Fig.~\ref{fig:examples} and supplementary materials for technical details).

Psychological studies \cite{mannan1997fixation,groner1984looking} have shown that when humans look through multiple visual cues that simultaneously demand attention, the order in which humans look at those cues is highly subjective. Therefore, by aggregating gazes of independent observers, we could record multiple important visual cues in one frame. In addition, it has been shown that human drivers look at buildings, trees, flowerbeds, and other unimportant objects non-negligibly frequently \cite{dreyeve2016}. Presumably, these eye movements should be regarded as noise for driving-related machine learning purposes. By averaging the eye movements of independent observers, we were able to effectively wash out those sources of noise (see Fig.~\ref{fig:incar_inlab}B).

\myparagraph{Comparison with In-Car Attention Data: } We collected in-lab driver attention data using videos from the DR(eye)VE dataset. This allowed us to compare in-lab and in-car attention maps of each video. The DR(eye)VE videos we used were 200 randomly selected 10-second video clips, half of them containing braking events and half without braking events.

We tested how well in-car and in-lab attention maps highlighted driving-relevant objects. We used YOLO~\cite{redmon2017yolo9000} to detect the objects in the videos of our dataset. We identified three object categories that are important for driving and that had sufficient instances in the videos (car, pedestrian and cyclist). We calculated the proportion of attended objects out of total detected instances for each category for both in-lab and in-car attention maps (see supplementary materials for technical details). The results showed that in-car attention maps highlighted significantly less driving-relevant objects than in-lab attention maps (see Fig.~\ref{fig:incar_inlab}A).

\begin{figure}
\begin{center}
\includegraphics[width=\linewidth]{figures/incar_inlab.eps}
\end{center}
\caption{Comparison between in-car and in-lab driver attention maps. (A) Proportions of attended objects of different categories for in-car and in-lab driver attention maps. In-car attention maps tend to highlight significantly fewer driving-relevant objects than in-lab attention maps. (B) An example of in-car driver attention maps showing irrelevant regions. The in-lab attention map highlights the car in front and a car that suddenly backed up, while the in-car attention map highlights some regions of the building}
\label{fig:incar_inlab}
\end{figure}

The difference in the number of attended objects between the in-car and in-lab attention maps can be due to the fact that eye movements collected from a single driver do not completely indicate all the objects that demand attention in the particular driving situation. One individual's eye movements are only an approximation of their attention \cite{rizzolatti1987reorienting}, and humans can also track objects with covert attention without looking at them \cite{cavanagh2005tracking}. The difference in the number of attended objects may also reflect the difference between first-person driver attention and third-person driver attention. It may be that the human observers in our in-lab eye-tracking experiment also looked at objects that were not relevant for driving. We ran a human evaluation experiment to address this concern.

\myparagraph{Human Evaluation: } To verify that our in-lab driver attention maps highlight regions that should indeed demand drivers' attention, we conducted an online study to let humans compare in-lab and in-car driver attention maps. In each trial of the online study, participants watched one driving video clip three times: the first time with no edit, and then two more times in random order with overlaid in-lab and in-car attention maps, respectively. The participant was then asked to choose which heatmap-coded video was more similar to where a good driver would look. In total, we collected 736 trials from 32 online participants. We found that our in-lab attention maps were more often preferred by the participants than the in-car attention maps (71\% versus 29\% of all trials, statistically significant as $p$ = 1$\times$10$^{-29}$, see Table~\ref{tab:human_evaluation}). Although this result cannot suggest that in-lab driver attention maps are superior to in-car attention maps in general, it does show that the driver attention maps collected with our protocol represent where a good driver should look from a third-person perspective.

\begin{table}
    \caption{Two human evaluation studies were conducted to compare in-lab human driver attention maps with in-car human driver attention maps and attention maps predicted by our HWS model, respectively. In-car human driver attention maps were preferred in significantly less trials than the in-lab human driver attention maps. The attention maps predicted by our HWS model were not preferred in as many trials as the in-lab human driver attention maps, but they achieved significantly higher preference rate than the in-car human driver attention maps}
    \label{tab:human_evaluation}
	\begin{center}
    	\resizebox{0.5\linewidth}{!}{%
    	\begin{tabular}{@{}cclc@{}} \toprule
            & \# trials & Attention maps & Preference rate \\ \midrule
            \multirow{2}{*}{Study 1}& \multirow{2}{*}{736}& in-car human driver& 29\%\\ \cmidrule{3-4}
            & & in-lab human driver& 71\%\\ \midrule
            \multirow{2}{*}{Study 2}& \multirow{2}{*}{462}& HWS model predicted& 41\%\\ \cmidrule{3-4}
            & & in-lab human driver& 59\%\\ \bottomrule
        \end{tabular}}

    \end{center}
\end{table}

In addition, we will show in the Experiments section that in-lab attention data collected using our protocol can be used to train a model to effectively predict actual, in-car driver attention. This result proves that our dataset can also serve as a substitute for in-car driver attention data, especially in crucial situations where in-car data collection is not practical.

To summarize, compared with driver attention data collected in-car, our dataset has three clear advantages: multi-focus, little driving-irrelevant noise, and efficiently tailored to crucial driving situations.

\section{Attention Prediction Model}
\subsection{Network Configuration}

\begin{figure}[t!]
\centering
\includegraphics[width=\linewidth]{figures/structure2.png} %
\caption{An overview of our proposed model that predicts human driver’s attention from input video frame. We use AlexNet pre-trained on ImageNet as a visual feature extractor. We also use three fully convolutional layers (Conv2D) followed by a convolutional LSTM network (Conv2D LSTM)}
 \label{fig:structure}
\centering
\end{figure}

Our goal is to predict the driver attention map for a video frame given the current and previous video frames. Our model structure can be divided into a visual feature extraction module, a visual feature processing module, and a temporal processing module (Fig.~\ref{fig:structure}).

The visual feature extraction module is a pre-trained dilated fully convolutional neural network, and its weights are fixed during training. We used ImageNet pre-trained AlexNet \cite{krizhevsky2012imagenet} as our visual feature extraction module. We chose to use the features from the conv5 layer. In our experiment, the size of the input was set to $1024\times576$ pixels, and the feature map by AlexNet was upsampled to $64\times36$ pixels and then fed to the following visual feature processing module.

The visual feature processing module is a fully convolutional neural network. It consists of three convolutional layers with $1 \times 1$ kernels and a dropout layer after each convolutional layer. It further processes the visual features from the previous extraction module and reduces the dimensionality of the visual features from 256 to 8. In our experiments, we observed that without the dropout layers, the model easily got stuck in a suboptimal solution which simply predicted a central bias map, i.e. an attention map concentrated in a small area around the center of the frame.

The temporal processor is a convolutional LSTM network with a kernel size of $3\times3$ followed by a Gaussian smooth layer ($\sigma$ set to 1.5) and a softmax layer. It receives the visual features of successive video frames in sequence from the visual feature processing module and predicts an attention map for every new time step. Dropout is used for both the linear transformation of the inputs and the linear transformation of the recurrent states. We had also experimented with using an LSTM network for this module and observed that the model tended to incorrectly attend to only the central region of the video frames. The final output of this model is a probability distribution over $64\times36$ grids predicting how likely each region of the video frame is to be looked at by human drivers. Cross-entropy is chosen as the loss function to match the predicted probability distribution to the ground-truth.

\begin{figure}[t!]
\centering
\includegraphics[height=3cm]{figures/weights.pdf} %
\caption{Human Weighted Sampling: (A) For each video frame, we measure the KL divergence between the collected driver attention maps and the mean attention map for that entire video clip ($\approx$10s). We use this computed KL divergence as a weight value to sample image frames during training phase, \textit{i}.\textit{e}.,~ training a model more often with uncommon attention maps. Histograms show that more uncommon attention maps were selected for training the model, \textit{e}.\textit{g}.,~ seeing pedestrians or traffic lights is weighted more than just seeing the vanishing point of roads. (B) Normalized sampling weights as a function of KL divergence values. A normalized sampling weight value of 1 indicates that the video frame is sampled once on average during a single epoch}
 \label{fig:kl}
\centering
\end{figure}

\subsection{Human Weighted Sampling (HWS)}
Human driver attention datasets, as well as many other driving related datasets, share a common bias: the vast majority of the datasets consist of simple driving situations such as lane-following or car-following. The remaining small proportion of driving situations, such as pedestrians darting out, traffic lights changing, etc., are usually more crucial, in the sense that making errors in these moments would lead to greater cost. Therefore, ignoring this bias and simply using mean prediction error to train and test models can be misleading. In order to tackle this problem, we developed a new method that uses human gaze data to determine the importance of different frames of a driving dataset and samples the frames with higher importance more frequently during training.

In simple driving situations human drivers only need to look at the center of the road or the car in front, which can be shown by averaging the attention maps of all the frames of one driving video. When the attention map of one frame deviates greatly from the average default attention map, it is usually an important driving situation where the driver has to make eye movements to important visual cues. Therefore, the more an attention map varies from the average attention map of the video, the more important the corresponding training frame is. We used the KL divergence to measure the difference between the attention map of a particular frame and the average attention map of the video. The KL divergence determined the sampling weight of this video frame during training.

The histogram of the KL divergence of all the training video frames of our dataset is shown in Fig.~\ref{fig:kl}. As we expected, the histogram was strongly skewed to the left side. Our goal was to boost up the proportion of the frames of high KL divergence values by weighted sampling. The sampling weight was determined as a function of KL divergence (D$_{\textnormal{KL}}$) illustrated in Fig.~\ref{fig:kl}B. The middle part of this function (D$_{\textnormal{KL}}$ $\in$ [1,3]) was set to be proportional to the inverse of the histogram so that after weighted sampling the histogram of KL divergence would become flat on this range. The left part of the function (D$_{\textnormal{KL}}<$1) was set to a low constant value so that those frames would be sampled occasionally but not completely excluded. The right part of the function was set to a saturated constant value instead of monotonically increasing values in order to avoid over-fitting the model to this small proportion of data. Besides, the attention maps collected in the beginning and the end of each video clip can deviates from the average default attention map merely because the participants were distracted by the breaks between video clips. We therefore restricted the sampling weights of the first second and the last 0.5 seconds of each video to be less or equal to once per epoch. The histogram of KL divergence after weighted sampling is shown in Fig.~\ref{fig:kl}A. In our experiment, we needed to sample the training frames in continuous sequences of 6 frames. For a particular sequence, its sampling weight was equal to the sum of the sampling weights of its member frames. These sequences were sampled at probabilities proportional to the sequence sampling weights.

\section{Results and Discussion}

Here, we first provide our training and evaluation details, then we summarize the quantitative and qualitative performance comparison with existing gaze prediction models and variants of our model. To test how natural and reasonable our model prediction look to humans, we conduct a human evaluation study and summarize the results. We further test whether our model trained on in-lab driver attention data can also predict driver attention maps collected in-car.

\subsection{Training and Evaluation Details}
We made two variants of our model. One was trained with a regular regime, \textit{i}.\textit{e}.,~ equal sampling during training, and the other was trained with Human Weighted Sampling (HWS). Except for the sampling method during training, our default model and HWS model shared the same following training settings. We used 926 videos from our BDD-A dataset as the training set and 306 videos as the testing set. We downsampled the videos to $1024 \times 576$ pixels and 3Hz. After this preprocessing, we had about 30k frames in our training set and 10k frames in our testing set. We used cross-entropy between predicted attention maps and human attention maps as the training loss, along with Adam optimizer (learning rate $=0.001$, $\beta_{1}=0.9$, $\beta_{2}=0.999$, $\epsilon=1\times10^{-8}$). Each training batch contained 10 sequences and each sequence had 6 frames. The training was done for 10,000 iterations. The two models showed stabilized testing errors by iteration 10,000.

To our knowledge, \cite{palazzi2017learning} and \cite{tawari2017computational} are the two deep neural models that use dash camera videos alone to predict human driver's gaze. They demonstrated similar results and were shown to surpass other deep learning models or traditional models that predict human gaze in non-driving-specific contexts. We chose to replicate \cite{palazzi2017learning} to compare with our work because their prediction code is public. The model designed by \cite{palazzi2017learning} was trained on the DR(eye)VE dataset \cite{dreyeve2016}. We will refer to \cite{palazzi2017learning}'s model as DR(eye)VE model in the following. The training code of \cite{palazzi2017learning} is not available. We implemented code to fine-tune their model on our dataset, but the fine-tuning did not converge to any reasonable solution, potentially due to some training parameter choices that were not reported. We then tested their pre-trained model directly on our testing dataset without any training on our training dataset. Since the goal of the comparison was to test the effectiveness of the combination of model structure, training data and training paradigm as a whole, we think it is reasonable to test how well DR(eye)VE model performs on our dataset without further training. For further comparison, we fine-tuned a publicly available state-of-the-art image gaze prediction model, SALICON \cite{huang2015salicon} on our dataset. We used the open source implementation \cite{christopherleethomas2016}. We also tested our models against a baseline model that always predicts the averaged human attention map of training videos.

Kullback-Leibler divergence (KL divergence, $D_{KL}$), Pearson’s Correlation Coefficient (CC), Normalized Scanpath Saliency (NSS) and Area under ROC Curve (AUC) are four commonly used metrics for attention map prediction \cite{palazzi2017learning,tawari2017computational,bylinskii2018different}. We calculated the mean prediction errors in these four metrics on the testing set to compare the different models. In order to test how well the models perform at important moments where drivers need to watch out, we further calculated the mean prediction errors on the subset of testing frames where the attention maps deviate significantly from the average attention maps of the corresponding videos (defined as KL divergence greater than 2.0). We will refer to these frames as non-trivial frames. Our models output predicted attention maps in the size of $64 \times 36$ pixels, but the DR(eye)VE model and the SALICON outputs in bigger sizes. For a fair comparison, we scaled the DR(eye)VE model and the SALICON model's predicted attention maps into $64 \times 36$ pixels before calculating the prediction errors.

Another important evaluation criterion of driver attention models is how successfully they can attend to the objects that demand human driver's attention, e.g. the cars in front, the pedestrians that may enter the roadway, etc. Therefore, we applied the same attended object analysis described in the Berkeley DeepDrive Attention Dataset section. We used YOLO~\cite{redmon2017yolo9000} to detect the objects in the videos of our dataset. We selected object categories that are important for driving and that have enough instances in both our dataset and the DR(eye)VE dataset for comparison (car, pedestrian and cyclist). We calculated the proportions of all the detected instances of those categories that were actually attended to by humans versus the models. The technical criterion of determining attended objects was the same as described in the Berkeley DeepDrive Attention Dataset section.

\setlength{\tabcolsep}{4pt}
\begin{table}[t]
    \caption{Performance comparison of human attention prediction. Mean and 95\% bootstrapped confidence interval are reported}
    \label{tab:exp1_metrics}
	\begin{center}
    	\resizebox{1.0\linewidth}{!}{%
    	\begin{tabular}{@{}lcccccccc@{}} \toprule
        \multirow{3}{*}{} & \multicolumn{4}{c}{Entire testing set} & \multicolumn{4}{c}{Testing subset where D$_{\textnormal{KL}}$(GT, Mean) $>$ 2} \\
        & \multicolumn{2}{c}{KL divergence} & \multicolumn{2}{c}{Correlation coefficient} & \multicolumn{2}{c}{KL divergence} & \multicolumn{2}{c}{Correlation coefficient}\\ \cmidrule{2-9}
        & Mean        & 95\% CI             & Mean             & 95\% CI & Mean        & 95\% CI             & Mean             & 95\% CI                   \\ \midrule
Baseline                              & 1.50        & (1.45, 1.54)        & 0.46             & (0.44, 0.48) & 1.87        & (1.80, 1.94)        & 0.36             & (0.34, 0.37)             \\
SALICON~\cite{huang2015salicon}                               & 1.41        & (1.39, 1.44)        & 0.53             & (0.51, 0.54) & 1.76        & (1.72, 1.80)        & 0.39             & (0.37, 0.41)\\
DR(eye)VE~\cite{palazzi2017learning}                             & 1.95        & (1.87, 2.04)        & 0.50             & (0.48, 0.52)            & 2.63        & (2.51, 2.77)        & 0.35             & (0.33, 0.37)             \\
Ours (default)                        & 1.24        & (1.21, 1.28)        & 0.58             & (0.56, 0.59)             & 1.71        & (1.65, 1.79)        & 0.41             & (0.40, 0.43)             \\
Ours (HWS)                            & {\bf{1.24}}        & {\bf{(1.21, 1.27)}}        & {\bf{0.59}}             & {\bf{(0.57, 0.60)}}             & {\bf{1.67}}        & {\bf{(1.61, 1.73)}}        & {\bf{0.44}}             & {\bf{(0.42, 0.45)}}             \\ \bottomrule
        \end{tabular}}
    \end{center}
\end{table}
\setlength{\tabcolsep}{1.4pt}

\subsection{Evaluating Attention Predictor}
\myparagraph{Quantitative Analysis of Attention Prediction: } The mean prediction errors of different models are summarized in Table~\ref{tab:exp1_metrics} (measured in $D_{KL}$ and CC) and Table S1 (measured in NSS and AUC) in supplementary materials. Both of our models significantly outperformed the DR(eye)VE model, the SALICON model and the baseline model in all metrics on both the entire testing set and the subset of non-trivial frames. Our model trained with HWS was essentially trained on a dataset whose distribution was altered from the distribution of the testing set. However, our HWS model showed better results than our default model even when being tested on the whole testing set. When being tested on the subset of non-trivial frames, our HWS model outperformed our default model even more significantly. These results suggest that HWS has the power to overcoming the dataset bias and better leveraging the knowledge hidden in crucial driving moments.

The results of the attended object analysis are summarized in Fig.~\ref{fig:model_object}A. Cars turned out to be easy to identify for all models. This is consistent with the fact that a central bias of human attention is easy to learn and cars are very likely to appear in the center of the road. However, for pedestrians and cyclists, the DR(eye)VE model, SALICON model and baseline model all missed a large proportion of them compared with human attention ground-truth. Both of our models performed significantly better than all the other competing models in the categories of pedestrians and cyclists, and our HWS model matched the human attention performances the best.

Importantly, our HWS model did not simply select objects according to their categories like an object detection algorithm. Considering the category that has the highest safety priority, pedestrian, our models selectively attended to the pedestrians that were also attended to by humans. Let us refer to the pedestrians that were actually attended to by humans as the important pedestrians and the rest of them as non-important pedestrians. Among all the pedestrians detected by the object detection algorithm, the proportion of important pedestrians was $33\%$. If our HWS model were simply detecting pedestrians at a certain level and could not distinguish between important pedestrians and non-important pedestrians, the proportion of important pedestrians among the pedestrians attended to by our model should also be $33\%$. However, the actual proportion of important pedestrians that our HWS model attended to was  $48\%$ with a bootstrapped $95\%$ confidence interval of $[42\%, 55\%]$. Thus, our HWS model predicts which of the pedestrians are the ones most relevant to human drivers.

\begin{figure}[t]
\centering
\includegraphics[width=\linewidth]{figures/model_object.eps}
\caption{Analysis of attended objects for human attention and different models tested on our dataset (A) and the DR(eye)VE dataset (B). Error bars show 95\% bootstrapped confidence intervals}
\label{fig:model_object}
\end{figure}

\myparagraph{Qualitative Analysis of Attention Prediction: } Some concrete examples are shown in Figure~\ref{fig:examples} (see supplementary information for videos). These examples demonstrates some important driving scenarios: pedestrian crossing, cyclist getting very close to the vehicle and turning at a busy crossing. It can be seen from these examples that the SALICON model and the DR(eye)VE model mostly only predicted to look at the center of the road and ignored the crucial pedestrians or cyclists. In the examples of row 1, 2 and 3, both our default model and HWS model successfully attended to the important pedestrian/cyclist, and did not give false alarm for other pedestrians who were not important for the driving decision. In the challenging example shown in row 4, the driver was making a right turn and needed to yield to the crossing pedestrian. Only our HWS model successfully overcame the central bias and attended to the pedestrian appearing in a quite peripheral area in the video frame.

\begin{figure}[t]
\centering
\includegraphics[width=\linewidth]{figures/examples.pdf}
\caption{Examples of the videos in our dataset, ground-truth human attention maps and the prediction of different models. The red rectangles in the original video column highlight the pedestrians that pose a potential hazard. Row 1: the driver had the green light, but a pedestrian was about to cross the road while speaking on a phone without looking at the driver. Another pedestrian was present in the scene, but not relevant to the driving decision. Row 2: the driver had a yellow light and some pedestrians were about to enter the roadway. Another pedestrian was walking in the same direction as the car and therefore not relevant to the driving decision. Row 3: a cyclist was very close to the car. Row 4: the driver was making a right turn and needed to yield to the crossing pedestrian. Other pedestrians were also present in the scene but not relevant to the driving decision}
\label{fig:examples}
\end{figure}

\myparagraph{Human Evaluation: } To further test how natural and reasonable our HWS model's predicted attention maps look to humans, we conducted an online Turing Test. In each trial, a participant watched one driving video clip three times: the first time with no edit, and then two times in random order with the ground-truth human driver attention map and our HWS model's predicted attention map overlaid on top, respectively. The participant was then asked to choose whether the first or the second attention map video was more similar to where a good driver would look.

Note that the experiment settings and instructions were the same as the online study described in the dataset section, except that one compares model prediction against the in-lab driver attention maps, and the other compares the in-car driver attention maps against the in-lab driver attention maps. Therefore, the result of this Turing Test can be compared with the result of the previous online study. In total, we collected 462 trials from 20 participants. If our HWS model's predicted attention maps were perfect and indistinguishable from the ground-truth human driver attention maps, the participants would had to make random choices, and therefore we would expect them to choose our model prediction in about $50\%$ of the trials. If our HWS model's prediction was always wrong and unreasonable, we would expect a nearly zero chosen rate for our model prediction. Our results showed that in 41\% of all trials the participants chose our HWS model's predicted attention maps as even better than the in-lab human attention maps (see Table~\ref{tab:human_evaluation}). In the previous online study, the in-car attention maps of DR(eye)VE only achieved a chosen rate of 29\%. This result suggests that our HWS model's predicted attention maps were even more similar to where a good driver should look than the human driver attention maps collected in-car (permutation test p = $4 \times 10^{-5}$).

\subsection{Predicting In-Car Driver Attention Data:}
To further demonstrate that our model has good generalizability and that our driver attention data collected in-lab is realistic, we conducted a challenging test: we trained our model using only our in-lab driver attention data, but tested it on the DR(eye)VE dataset, an in-car driver attention dataset. Note that the DR(eye)VE dataset covers freeway driving, which is not included in our dataset due to the small density of road user interactions on freeway. The high driving speed on freeway introduces strong motion blur which is not present in our dataset videos. Furthermore, drivers need to look further ahead in high speed situations, so the main focus of driver gaze pattern shifts up as the driving speed increases. In order to adapt our model to these changes, we selected 200 ten-second-long video clips from the training set of the DR(eye)VE dataset and collected in-lab driver attention maps for those video clips (already described in the Berkeley DeepDrive Attention Dataset section). We fine-tuned our HWS model with these video clips (30 minutes in total only) and the corresponding in-lab driver attention maps, and then tested the model on the testing set of the DR(eye)VE dataset (with in-car attention maps). The mean testing errors were calculated in $D_{KL}$ and CC because the calculation of NSS and AUC requires the original fixation pixels instead of smoothed gaze maps and the original fixation pixels of the DR(EYE)VE dataset were not released. Our fine-tuned model showed a better mean value in KL Divergence and a worse mean value in CC than the DR(eye)VE model (see Table~\ref{tab:exp4_metrics}). But the 95\% bootstrapped confidence intervals for the two models in both metrics overlapped with each other. So overall we concluded that our fine-tuned model matched the performance of the DR(eye)VE model. Note that the DR(eye)VE model was trained using the DR(eye)VE dataset and represents the state-of-the-art performance on this dataset.

\setlength{\tabcolsep}{4pt}
\begin{table}[t]
    \caption{Test  results  obtained on the DR(eye)VE dataset by  the  state-of-the-art  model (DR(eye)VE) and our finetuned model. Mean and 95\% bootstrapped confidence interval are reported}
    \label{tab:exp4_metrics}
	\begin{center}
    	\resizebox{0.6\linewidth}{!}{%
    	\begin{tabular}{@{}lcccc@{}} \toprule
        \multirow{2}{*}{} & \multicolumn{2}{c}{KL divergence} & \multicolumn{2}{c}{Correlation coefficient} \\ \cmidrule{2-5}
        & Mean        & 95\% CI             & Mean             & 95\%CI                   \\ \midrule
DR(eye)VE                               & 1.76        & (1.65, 1.87)        & 0.54             & (0.51, 0.56)             \\
Ours (finetuned)                      & 1.72        & (1.66, 1.81)        & 0.51             & (0.48, 0.53)             \\ \bottomrule
        \end{tabular}}

    \end{center}
\end{table}
\setlength{\tabcolsep}{1.4pt}

We also calculated proportions of attended objects of important categories for our fine-tuned model and the DR(eye)VE model (Fig.~\ref{fig:model_object}B). Our fine-tuned model showed significantly higher proportions of attended objects in the car, pedestrian and cyclist categories and was more similar to the in-lab driver attention than the DR(eye)VE model. Note that we have shown in the Berkeley DeepDrive Attention Dataset section that humans rated the in-lab attention maps as more similar to where a good driver should look from a third-person perspective than the in-car attention maps.

\section{Conclusions}
In this paper, we introduce a new in-lab driver attention data collection protocol that overcomes drawbacks of in-car collection protocol. We contribute a human driver attention dataset which is to-date the richest and will be made public. We propose Human Weighted Sampling which can overcome common driving dataset bias and improve model performance in both the entire dataset and the subset of crucial moments. With our dataset and sampling method we contribute a novel human driver attention prediction model that can predict both in-lab and in-car driver attention data. The model demonstrates sophisticated behaviors and show prediction results that are nearly indistinguishable from ground-truth to humans.

\vspace{5mm}
\myparagraph{Acknowledgement: }
This work was supported by Berkeley DeepDrive. We thank Yuan Yuan and Victor Hiltner for assistance with driver gaze data collection, and Victor Hiltner, Kevin Li and Drew Kaul for comments and proofreading that greatly improved the manuscript.
\vspace{3mm}

\printbibliography
\end{document}