\section{Introduction} \vspace{-4mm}
Kernel methods for classification and regression (and Support Vector Machines (SVMs) in particular) require selection of a kernel. Kernel Learning (KL) algorithms such as those found in~\cite{xu2010simple,sonnenburg2010shogun,yang2011efficient} automate this task by finding the kernel, $k \in \mcl K$ which optimizes an achievable metric such as the soft margin (for classification). The set of kernels, $k \in \mcl K$, over which the algorithm can optimize, however, strongly influences the performance and robustness of the resulting classifier or predictor.

To understand how the choice of $\mcl K$ influences performance and robustness, three properties were proposed in~\cite{JMLR} to characterize the set $\mcl K$ - tractability, density, and universality. Specifically, $\mcl K$ is tractable if $\mcl K$ is convex (or, preferably, a linear variety) - implying the KL problem is solvable using, e.g.~\cite{rakotomamonjy_2008,jain_2012,lanckriet_2004,qiu2005multiple,gonen2011multiple}. The set $\mcl K$ has the density property if, for any $\epsilon>0$ and any positive kernel, $k^*$ there exists a $k \in \mcl K$ where $\norm{k-k^*}\le \epsilon$. The density property implies the kernel will perform well on untrained data (robustness or generalizability). The set $\mcl K$ has the universal property if any $k \in \mcl K$ is universal - ensuring the classifier/predictor will perform arbitrarily well on large sets of training data.

In \cite{JMLR}, the Tessellated Kernels (TKs) were shown to have all 3 properties, the first known such class of kernels. This work was based on a general framework for using positive matrices to parameterize positive kernels (as opposed to positive kernel matrices as in~\cite{lanckriet_2004,qiu2005multiple,ni2006learning}). Unfortunately, however, the algorithms proposed in~\cite{JMLR} were either based on SemiDefinite Programming (SDP) (thereby limiting the amount of training data) or used a randomized linear basis for the kernels (implying loss of density).
Thus, while the algorithms in~\cite{JMLR} outperformed all other methods (including deep learning) as measured by Test Set Accuracy (TSA), the computation times were not competitive. Furthermore, the results in~\cite{JMLR} did not encompass the problem of regression.

In this paper, we extend the TK framework proposed in~\cite{JMLR} to the problem of regression. The KL problem in regression has been studied using SDP in~\cite{qiu2005multiple,ni2006learning} and Quadratic Programming (QP) in e.g.~\cite{rakotomamonjy_2008,jain_2012}. However, neither of these previous works considered a set of kernels with both the tractability and the density property. By generalizing the Tessellated KL framework proposed in~\cite{JMLR} to the regression problem, we demonstrate significant increases in performance, as measured by Mean Square Error (MSE), and when compared to the results in~\cite{rakotomamonjy_2008,jain_2012,qiu2005multiple}.

In addition, we show that the SDP-based algorithm~\cite{JMLR} for classification, and extended here to regression, can be decomposed into primal and dual sub-problems, $OPT\_A$ and $OPT\_P$ - similar to the approach taken in~\cite{rakotomamonjy_2008,jain_2012}. Furthermore, we show that $OPT\_P$ (an SDP) admits an analytic solution using the Singular Value Decomposition (SVD) - an approach which allows us to consider higher dimensional feature spaces and more complex TKs. In addition, $OPT\_A$ is a convex QP and may be solved efficiently with achieved complexity which scales as $O(m^{2.16})$ where $m$ is the number of data points. We use a two-step algorithm on $OPT\_A$ and $OPT\_P$ and show that termination at $OPT\_A=OPT\_P$ is equivalent to global optimality. The resulting algorithm, then, does not require the use of SDP and, when applied to several standard test cases, is shown to retain the favorable TSA of~\cite{JMLR} for classification, while offering improved MSE for regression, and competitive computation times as compared to other KL and deep learning algorithms. \vspace{-3mm}
\section{An Ideal Set of Kernels for KL in Classification and Regression} \label{sec:Ideal} \vspace{-4mm}
Consider a generalized representation of the KL problem, which encompasses both classification and regression where (using the representor theorem~\cite{scholkopf2001generalized}) the learned function is of the form $f_{\alpha,k}(z) = \sum_{i=1}^m \alpha_i k(x_i,z)$.\vspace{-2mm}
\begin{equation} \label{eqn:OPT}
\min_{k  \in \mcl{K}} \min_{\alpha \in \mathbb{R}^m,b} \norm{f_{\alpha,k}}^2 + C \sum\nolimits_{i=1}^m l(f_{\alpha,k},b)_{y_i,x_i}
\end{equation}
Here $\norm{f_{\alpha,k}}_X = \sum_{i=1}^m \sum_{j=1}^m \alpha_i \alpha_j k(x_i,x_j)$ and $l(f_{\alpha,k},b)_{y_i,x_i}$ is the loss function and is defined for SVM binary classification and SVM regression as $l^c(f_{\alpha,k},b)_{y_i,x_i}$ and $l^r(f_{\alpha,k},b)_{y_i,x_i}$, respectively, where
\[
l^c(f_{\alpha,k},b)_{y_i,x_i} \hspace{-1mm}=  \max \{  0, 1 - y_i (f_{\alpha,k}(x_i) - b) \}~~\text{and}~~ l(f_{\alpha,k},b)^r_{y_i,x_i} \hspace{-1mm}=  \max \{  0, |y_i - (f_{\alpha,k}(x_i) - b)|- \epsilon \}.
\]

The properties of the classifier/predictor, $f_{\alpha,k}$, resulting from Optimization Problem~\ref{eqn:OPT} will depend on the properties of the set $\mcl K$, which is presumed to be a subset of the convex cone of all positive kernels. To understand how $\mcl K$ influences the tractability of the optimization problem and the resulting fit, we consider three properties of the set, $\mcl K$.

\subsection{Tractability} \vspace{-3mm}

We say a set of kernel functions, $\mcl K$, is tractable if it can be represented using a countable basis.
\begin{defn} \label{Tractability}
The set of kernels $\mcl K$ is \textit{\textbf{tractable}} if there exist a countable set $\{G_{i}(x,y)\}_{i}$ such that, for any $k \in \mcl K$, there exists $N_G \in \mathbb{N}$ where $k(x,y)=\sum_{i=1}^{N_G} v_i G_i(x,y)$ for some $v \in \mathbb{R}^{N_G}$.
\end{defn}
Note the $G_i(x,y)$ need not be positive kernel functions. The tractable property is required for the KL problem to be tractable using algorithms for convex optimization.

\subsection{Universality}\label{Universal_Motivation} \vspace{-3mm}
Universal kernel functions always have positive definite (full rank) kernel matrices, implying that for arbitrary data $\{y_i,x_i\}_{i=1}^m$, there exists a function $f(z) = \sum_{i=1}^m \alpha_i k(x_i,z)$, such that $f(x_j) = y_j$ for all $j= 1,..,m$. Conversely, if a kernel is not universal, then exists a data set $\{x_i,y_i\}_{i=1}^m$ such that for any $\alpha \in \mathbb{R}^m$, there exists some $j\in \{1,\cdots,m\}$ such that $f(y_j) \neq \sum_{i=1}^m \alpha_i k(x_i,x_j)$. This ensures that SVMs using universal kernels can always benefit from additional training data, whereas non-universal kernels may saturate.

\begin{defn}
A kernel $k:X \times X \rightarrow \mathbb{R}$ is said to be \textit{universal} on the compact metric space $X$ if it is continuous and there exists an inner-product space $\mcl W$ and feature map, $\Phi : X \rightarrow \mcl W$ such that $k(x,y)=\ip{\Phi(x)}{\Phi(y)}_{\mcl W}$ and where the unique Reproducing Kernel Hilbert Space (RKHS), $\mcl H:=\{f\;:\;f(x)=\ip{v}{\Phi(x)},\; v \in \mcl W\}$
with associated norm $\norm{f}_{\mcl H} :=\inf_{v}\{\norm{v}_{\mcl W}\;:\; f(x)=\ip{v}{\Phi(x)}\}$ is dense in $\mcl C(X):=\{f \,: \, X \rightarrow \mathbb{R} \;:\; $f$\, \text{ is continuous}\}$ where $\norm{f}_{\mcl C}:=\sup_{x\in X} |f(x)|$.
\end{defn}
The following definition extends the universal property to a set of kernels.
\begin{defn}
A set of kernel functions $\mcl K$ has the universal property if every kernel function $k \in \mcl K$ is universal.
\end{defn}

\subsection{Density} \vspace{-3mm}
The third property is density which distinguishes the TK class from other sets of kernel functions with the universal property. For instance consider a set containing a single Gaussian kernel function - which is clearly not ideal for kernel learning.  The set containing a single Gaussian is tractable (it has only one element) and every member of the set is universal. However, it is not dense.

Considering SVM for classification, the KL problem determines the kernel $k \in \mcl K$ for which we may obtain the maximum separation in the kernel-associated feature space. Increasing this separation distance makes the resulting classifier more robust (generalizable)~\cite{boehmke2019hands}. The density property, then, ensures that the resulting KL algorithm will be maximally robust (generalizable) in the sense of separation distance.

Likewise, considering SVMs for regression, the KL problem finds the kernel $k\in \mcl K$ which permits the ``flattest''~\cite{smola2004tutorial} function in feature space. In this case, the density property ensures that the resulting KL algorithm will be maximally robust (generalizable) in the sense of flatness.

These arguments motivate the following definition of the pointwise density property.
\begin{defn}
The set of kernels $\mcl K$ is said to be \textbf{pointwise dense} if for any positive kernel, $k^*$, any set of data $\{x_i\}_{i=1}^m$, and any $\epsilon>0$, there exists $k \in \mcl K$ such that $\norm{k(x_i,x_j)-k^*(x_i,x_j)}\le \epsilon$.
\end{defn}
\vspace{-3mm}
\section{A General Framework for Representation of Tractable Kernel Sets} \label{sec:4} \vspace{-4mm}
Here we define a framework for constructing classes of tractable positive kernel functions and illustrate this approach on the class of General Polynomial Kernels.

\begin{prop}\label{prop:kernel}
Let $N$ be any bounded measurable function $N: X \times Y \rightarrow \mathbb{R}^q$ and $P \in \mathbb{R}^{q \times q}$ be a positive semidefinite matrix $P\geq 0$. Then
\begin{equation}
k(x,y)=\int_{X} N(z,x)^T P N(z,y) dz\label{eqn:kernel}\vspace{-2mm}
\end{equation}
is a positive kernel function.
\end{prop}
The proof for Proposition~\eqref{prop:kernel} may be found in \cite{JMLR}.

\begin{lem}\label{lem:tractable}
Let $N$ be any bounded measurable function $N: X \times Y \rightarrow \mathbb{R}^q$ on compact $X$ and $Y$. Then the set of kernel functions
\begin{equation}
{\mcl K} := \left\{ k ~|~ k(x,y)=\int_{X} N(z,x)^T P N(z,y) dz, ~ P \geq 0 \right\} \qquad \text{is tractable.}\label{eqn:tractable}
\end{equation}
\end{lem}
For a given $N$, the map $P \mapsto k$ is linear. Specifically,
\[
k(x,y)=\sum\nolimits_{i=1}^q \sum\nolimits_{j=1}^q P_{i,j} G_{i,j}(x,y)~~\text{where}~~ G_{i,j}(x,y)=\int_{X}N_{i}(z,x)N_{j}(z,y)dz.
\]
and thus by Definition~\ref{Tractability} $\mcl K$ is tractable.

In Subsection~\ref{subsec:GPK} we apply this framework to obtain Generalized Polynomial Kernels. In Subsection~\ref{subsec:TK}, we use the framework to obtain the TK class.

\subsection{The Class of General Polynomial Kernels is Tractable}\label{subsec:GPK} \vspace{-3mm}
The class of General Polynomial Kernels (GPKs) is defined as the set of all polynomials, each of which is a positive kernel.
\begin{align}
 \mcl K_P &:= \{ k \in \mathbb{R}[x,y]\;:\; k \text{ is a positive kernel}  \} \label{polyKernel}
\end{align}
The GPK class is not universal, but is tractable, as per the following lemma.
\begin{lem}
$\mcl K_P$ is tractable.
\end{lem} \vspace{-4mm}

\begin{proof}
Let $Z_d: \mathbb{R}^n \rightarrow \mathbb{R}^q$ be the vector of monomials of degree $d$ or less.  From~\cite{JMLR}, we have that a polynomial $k$ of degree $2d$ is a positive polynomial kernel if and only if there exists some $P\geq 0$ such that $k(x,y)=Z_d(x)^T P Z_d(y)$. Now for any finite-dimensional subset of $\mcl K_P$, let $d$ be the maximum degree over this subset and define $N(z,y)=Z_d(y)$. Then Lemma~\ref{lem:tractable} implies that $\mcl K_P$ is tractable.
\end{proof}

This lemma implies that a representation of the form of Equation~\eqref{eqn:kernel} is necessary and sufficient for a GPK to be positive. For convenience, we denote the set of GPK kernels of degree $d$ or less as follows~\cite{recht2006convex}.
\begin{align}
 \mcl K_P^d &:= \{ k\;:\; k(x,y) = Z_d(x)^T P Z_d(y) \; : \; P \geq 0  \} \label{polyKernel}
\end{align}
\vspace{-3mm}

\section{Tessellated Kernels: Tractable, Dense and Universal} \label{sec:TK} \vspace{-4mm}
In this section, we define the class of TK kernels and show it is tractable, dense, and universal.
\subsection{Tessellated Kernels}\label{subsec:TK} \vspace{-3mm}
Again, let $Z_d: \mathbb{R}^n  \times \mathbb{R}^n \rightarrow \mathbb{R}^q$ be the vector of monomials of degree $d$. Define $\mathbf I$, the indicator function for the positive orthant, and the following choice of $N: \mathbb{R}^n \times \mathbb{R}^n \rightarrow \mathbb{R}^{2q}$ as \vspace{-2mm}
\begin{equation}
\mathbf I(z) = \begin{cases}
    1       & \quad z \ge 0\\
    0  & \quad \text{otherwise,}\\
\end{cases} \quad ~~ \text{and} ~~ \quad N^d_{T}(z,x) = \bmat{Z_d(z,x)\mathbf I(z-x) \\ Z_d(z,x)\mathbf I(x-z) } \label{eqn:N}
\end{equation}
\vspace{-2mm}
where $z\ge 0$ means $z_i \ge 0$ for all $i$.
We now define the set of TK kernels for $a<b \in \mathbb{R}^n$ as
\[
\mcl K^d_T:=\left\{k\;:\; k(x,y) = \int_{a}^b N_T^d(z,x)^T \; P\;  N^d_{T}(z,y) dz,\; P\ge 0\right\},~ \mcl K_T:=\{k\;:\; k\in \mcl K_T^d,\; d \in \mathbb{N}\}.
\]
Kernels in the TK class are ``Tessellated'' in the sense that each datapoint defines a vertex which bisects each dimension of the domain of the resulting classifier/predictor - resulting in a tessellated partition of the feature space.

\subsection{The Set of TK Kernels is Tractable} \vspace{-3mm}
The class of TK kernels is prima facie in the form of Eqn.~\eqref{eqn:tractable} in Lemma~\ref{lem:tractable} and hence is tractable.

However, we will expand on this result by specifying the basis for the set of TK kernels, which will then be used in Section~\ref{sec:2step}.

\begin{corr}\label{corr}
Suppose that for $a<b \in \mathbb{R}^n$, and $d \in\mathbb{N}$. We define the finite set $D_d:=\{(\delta,\lambda)\in \mathbb{N}^{2n} : \norm{(\delta,\lambda)}_1 \le d\}$.  Let $\{ [\delta_i, \gamma_i] \}_{i=1}^q \subseteq D_d$ be some ordering of $D_d$ and define $Z_d(x,z)_j = x^{\delta_j} z^{\gamma_j}$ where $z^{\delta_j} x^{\gamma_j}:=\prod_{i=1}^n z_i^{\delta_j,i}x_i^{\gamma_j,i}$. Now let $k$ be as defined in Eqn.~\eqref{eqn:kernel} for some $P>0$ and where $N$ is as defined in Eqn.~\eqref{eqn:N}. If we partition $P = \bmat{Q & R\\ R^T & S}$ then we have,\vspace{-2mm}
\[
k(x,y)=\sum\nolimits_{i,j=1}^q Q_{i,j} g_{i,j}(x,y) + R_{i,j}t_{i,j}(x,y) + R^T_{i,j}t_{i,j}(y,x)  + S_{i,j}  h_{i,j}(x,y)
\]
where $g_{i,j},t_{i,j},h_{i,j}:\mathbb{R}^{2n}\rightarrow \mathbb{R}$ are defined as
\begin{align}\label{g}
g_{i,j}(x,y) &:= x^{\delta_i}y^{\delta_j} T(p^*(x,y),b,\gamma_{i,j} + \mathbf{1} ), \; t_{i,j}(x,y) := x^{\delta_{i}}y^{\delta_{j}} T(x,b,\gamma_{i,j} + \mathbf{1}  ) - g_{i,j}(x,y), \; \text{and}\notag \\
 h_{i,j}(x,y) &:= x^{\delta_{i}}y^{\delta_{j}} T(a,b,\gamma_i + \gamma_j + \mathbf{1}  ) - g_{i,j}(x,y) - t_{i,j}(x,y) - t_{i,j}(y,x), \notag
\end{align}
where $\mathbf{1} \in \mathbb{N}^n$ is the vector of ones, $p^*:\mathbb{R}^{2n}\rightarrow \mathbb{R}^n$ is defined elementwise as
$p^*(x,y)_i = \max \{x_i,y_i \}$, and $T:\mathbb{R}^n \times \mathbb{R}^n \times \mathbb{N}^n\rightarrow \mathbb{R}$ is defined as \[T(x,y,\zeta) = \prod\nolimits_{j=1}^n \left(\hspace{2mm} \frac{y_j^{\zeta_j}}{\zeta_j}-\frac{x_j^{\zeta_j}}{\zeta_j}\right).
\]
\vspace{-5mm}
\end{corr}
The proof of Corollary~\ref{corr} can be found in~\cite{JMLR}.
\subsection{The TK Class is Dense} \vspace{-3mm}
The density property differentiates the set of TK kernels from other sets of kernel functions (e.g. a linear combination of Gaussian kernels of fixed bandwidths).

From \cite{JMLR} we have that the set of TK kernels satisfies the pointwise density property.
\begin{thm} \label{thm:TessellatedOptimal}
For any kernel matrix $K^*$ and any finite set $\{x_i\}_{i=1}^m$, there exists a $d \in\mathbb{N}$ and $k \in\mcl K_T^d$ such that if $K_{i,j} = k(x_i,x_j)$, then $K=K^*$.
\end{thm}

\begin{wrapfigure}{R}{0.32\textwidth}
\vspace{-1mm}
 \begin{center}
\includegraphics[trim= 20 0 50 20, clip, width=.32\textwidth]{Infnorm.eps}
  \end{center}\vspace{-3mm}
        \caption{The achieved objective ($\norm{K-K^*}_\infty$) of Optimization Problem~\ref{opt:infnorm} for TK and GPK of degree $d$; and $m$ Gaussians with bandwidths in $[.01,10]$.  The number of bandwidths is selected so that the number of decision variables match in the Gaussian and TK cases.
        }\label{fig:density}
\vspace{-7mm}
  \end{wrapfigure}

In~\cite{JMLR} an analytical solution, $K^*$, was found for the optimal trace-constrained kernel matrix that maximized the separation distance between two classes of points in the feature space. It was shown in this work that when $\{y_i\}$ has an equal number of positive and negative labels, $K^*$ contains an equal number of positive and negative elements - illustrating the importance of using kernels which are not pointwise positive (Gaussians are pointwise positive).

To illustrate the density property, then, we show how optimal GPK and TK kernels yield kernel matrices which approximate the analytic solution, $K^*$, of the optimal kernel matrix problem for a given set of data $\{x_i\}$ and labels $\{ y_i \}$, while Gaussian kernels do not. Specifically, we consider the following optimization problem. \vspace{-2mm}
\begin{equation}
\min_{k \in \mcl K} \norm{K-K^*}_\infty~~~ s.t. ~~K_{i,j} = k(x_i,x_j) \label{opt:infnorm}
\end{equation}
In these problems, the sets $\mcl K$ will be: $\mcl K_G^\gamma$ - the sum of $N$ Gaussians with bandwidths $\gamma_i$; $\mcl K_P^d$ - the GPKs of degree $d$; and $\mcl K_T^d$  - the TK kernels of degree $d$. More precisely, for bandwidths $\gamma\in \mathbb{R}^{N}$, we define $\mcl K^\gamma_G := \left\{k\;:\; k(x,y) = \sum_{i = 1}^N \mu_i \text{e}^\frac{||x-y||_2^2}{\gamma_i} \; : \; \mu_i > 0  \right\}.$

Consider a spiral data set with 20 samples, using equal numbers of positive and negative labels.
Fig.~\ref{fig:density} shows the achieved objective value of Problem~\eqref{opt:infnorm} for $\mcl K^\gamma_G$, $\mcl K_P^d$, and $\mcl K^d_T$ as a function of the number of bandwidths (top $x$ axis - $N$ in $\mcl K^\gamma_G$), polynomial degree (bottom $x$ axis - $d$ in $\mcl K_P^d$, and $\mcl K^d_T$). The $x$-axes of the plots are scaled to show equal numbers of decision variables. As expected,
the case $\mcl K=\mcl K_\gamma^G$ saturates with an objective value significantly larger than the lower bound.  The cases $\mcl K=\mcl K_P^d$ and $\mcl K=\mcl K^d_T$, meanwhile have almost no error at degree $d=7$.\vspace{-3mm}

\subsection{TK Kernels are Universal} \vspace{-3mm}
Finally we discuss the universality property of the class of TK kernels which ensures that every TK function can fit the training data well.

The following theorem from \cite{JMLR} shows that any TK kernel with $P > 0$ is necessarily universal.
\begin{thm} \label{thm:universal}
Suppose $k$ is as defined in Eqn.~\eqref{eqn:kernel} for some $P>0$, $d \in\mathbb{N}$ and $N$ as defined in Eqn.~\eqref{eqn:N}. Then $k$ is universal for $a<b \in \mathbb{R}^n$.
\end{thm}
This theorem implies that even if we use the subset of TK kernels defined by $d=0$, this subset is still universal.  \vspace{-3mm}

\section{A New Algorithm for KL in Classification and Regression using TKs}\label{sec:2step} \vspace{-4mm}
In this section, we express the KL optimization problem for both classification and regression and break this optimization problem into two sub-problems which allow us to express the problem in primal and dual form. For convenience, we define the feasible sets for the sub-problems as
\[
\mcl X:=\{P \in \mathbb{R}^{q \times q}\;:\; \text{trace}(P) = q,\; P > 0\}\vspace{-2mm}
\]
\[
\mcl Y_c:=\{\alpha \in \mathbb{R}^m\; : \; \sum\limits_{i=1}^m \alpha_iy_i = 0,\; 0 \leq \alpha_i \leq C \},\quad
\mcl Y_r:=\{\alpha \in \mathbb{R}^m\;:\; \sum\limits_{i=1}^m \alpha_i = 0, \; \alpha_i \in [-C, C]\}.
\]
The common part of the objective is
\begin{equation} \label{OalphaP}
O(\alpha,P):= -\frac{1}{2} \sum\nolimits_{i=1}^m\sum\nolimits_{j=1}^m \alpha_i \alpha_j \int_{a}^b N_T^d(z,x_i)^T P N^d_{T}(z,y_j) dz,
\end{equation}
while the unique parts of the objective are
\[
\kappa_c(\alpha):=\sum\nolimits_{i=1}^m \alpha_i\qquad\text{and}\qquad
\kappa_r(\alpha):=-\epsilon \sum\nolimits_{i=1}^m |\alpha_i| + \sum\nolimits_{i=1}^m y_i\alpha_i.
\]
Then the KL optimization problem ($OPT$) for TK kernels ($\odot$ being elementwise multiplication) is as follows for classification and regression, respectively.
\begin{align*
OPT:=\min_{P \in \mcl X} \max_{~\alpha \in \mcl Y_c}  \quad O(\alpha \odot y,P) + \kappa_c(\alpha), \quad \text{and} \quad OPT:=\min_{P \in \mcl X} \max_{~\alpha \in \mcl Y_r}  \quad O(\alpha,P) + \kappa_r(\alpha).
\end{align*}

\textbf{Primal Formulation:} We can formulate the primal problem ($OPT_P$) as
\begin{align}  \label{KernelSVC}
OPT_P=\min_{P \in \mcl X} \max_{\alpha \in \mcl Y} O(\alpha, P)+\kappa_r(\alpha)\; \text{or}\;O(\alpha\odot y, P)+\kappa_c(\alpha) =\min_{P \in \mcl X} ~~ OPT\_A(P)
\end{align}
where for classification and regression, respectively,
\[
OPT\_A(P):=\max_{\alpha \in \mcl Y_c} \;\; O(\alpha\odot y, P)+\kappa_c(\alpha), \quad \text{and} \quad OPT\_A(P):=\max_{\alpha \in \mcl Y_r}\;\; O(\alpha, P)+\kappa_r(\alpha).
\]

\textbf{Dual Formulation:} Alternatively, we have the dual formulation ($OPT_D$).
\begin{align}  \label{KernelSVC}
OPT_D=\max_{\alpha \in \mcl Y} OPT\_P( \alpha)
\end{align}
where $\mcl Y=\mcl Y_c$ for classification and $\mcl Y=\mcl Y_r$ regression. Likewise, for classification and regression, respectively,
\[
OPT\_{P}(\alpha):=\min_{P\in \mcl X}  O(\alpha \odot y, P)+\kappa_c(\alpha) \quad \text{and} \quad OPT\_{P}(\alpha):=\min_{P\in \mcl X}  O(\alpha, P) +\kappa_r(\alpha).
\]

\begin{lem}
For $\alpha \in \mcl Y$, $P \in \mcl X$, $OPT\_A(P)=OPT\_P(\alpha)$  if and only if: $\{\alpha,P\}$ solve $OPT$; $P$ solves $OPT_P$; and $\alpha$ solves $OPT_D$.
\end{lem} \vspace{-4mm}
\begin{proof}
For any minmax optimization problem with objective function $\phi$, we have
\[d^* = \max_{\alpha \in \mcl Y} \min_{P \in \mcl X} \phi(P,\alpha) \leq   \min_{P \in \mcl X}\max_{\alpha \in \mcl Y}  \phi(P,\alpha) = p^*, \]
and strong duality holds ($p^*-d^* = 0$) if $\mcl X$ and $\mcl Y$ are both convex and one is compact, $\phi(\cdot,\alpha)$ is convex for every $\alpha \in \mcl Y$ and $\phi(P,\cdot)$ is concave for every $P \in \mcl X$, and the function $\phi$ is continuous \cite{fan1953minimax}. In our case, these conditions hold for both classification and regression where $\phi(P,\alpha)=O(\alpha, P)+\kappa_r(\alpha)\; \text{or}\;O(\alpha\odot y, P)+\kappa_c(\alpha)$. Hence if $\alpha^*$ solves $OPT\_P$ and $P^*$ solves $OPT\_A$, then $\{\alpha^*, P^*\}$ solves $OPT$ and\vspace{-1.5mm}
\[
OPT\_P( \alpha^*)=\max_{\alpha \in \mcl Y} OPT\_P( \alpha)=\min_{P \in \mcl X} ~~ OPT\_A(P)=OPT\_A(P^*).\vspace{-2.5mm}
\]
Conversely, suppose $\alpha \in \mcl Y$, $P \in \mcl X$, then \vspace{1mm}
\begin{align*}\small
OPT\_P(\alpha) \le \max_{\alpha \in \mcl Y}\; OPT\_P(\alpha) &= OPT\_P( \alpha^*)\\[-2mm]
& = OPT\_A(P^*)= \min_{P \in \mcl X} ~~ OPT\_A(P) \le OPT\_A(P).\\[-6mm]
\end{align*}
Hence if $OPT\_A(P)=OPT\_P(\alpha)$, then $OPT\_A(P)=OPT\_A(P^*)=OPT\_P(\alpha^*)=OPT\_P(\alpha)$ and hence $P$ and $\alpha$ solve $OPT\_A$ and $OPT\_P$, respectively.
\end{proof}\vspace{-3mm}
We propose Algorithm~\ref{TKL} as a two-step iterative algorithm for solving Optimization Problem~\eqref{KernelSVC}.
\begin{wrapfigure}{R}{0.5\textwidth}
\vspace{-15mm}
\begin{minipage}{0.5\textwidth}
\begin{algorithm}[H]
\begin{algorithmic}
\STATE \texttt{Initialize} $P=I$;
\WHILE{$OPT\_P(\alpha_k)-OPT\_A(P_k)>\epsilon$}
\STATE $\alpha_{k+1}= OPT\_A(P_k)$
\STATE $P_{k+1}= \frac{P_K + t OPT\_P(\alpha_{k+1})}{1+t}$ (select $t$ using line search)
\STATE $k= k+1$
\ENDWHILE
\end{algorithmic}
\caption{Two Step TKL} \label{TKL}
\end{algorithm}
 \end{minipage}
\vspace{-4mm}
  \end{wrapfigure}

\begin{figure}[t]
    \centering
    \begin{subfigure}[t]{0.23\textwidth}
        \centering
\includegraphics[trim= 20 0 50 20, clip, width=.95\textwidth]{GC3.png}
        \caption{An image from Google Maps of a section of the Grand Canyon corresponding to (36.04, -112.05) latitude and (36.25, -112.3) longitude.}
    \end{subfigure}%
    ~
    \begin{subfigure}[t]{0.23\textwidth}
        \centering
\includegraphics[trim= 20 0 50 20, clip, width=0.95\textwidth]{grandCanyon3DData.eps}
        \caption{Elevation data ($m=750$) from \cite{becker2009global} for a section of the Grand Canyon between (36.04, -112.05) latitude and (36.25, -112.3) longitude.}
    \end{subfigure}
        ~
    \begin{subfigure}[t]{0.23\textwidth}
        \centering
\includegraphics[trim= 20 0 50 20, clip, width=0.95\textwidth]{grandCanyon3DGaussian.eps}
        \caption{Predictor using a hand-tuned Gaussian kernel trained on the elevation data in (b). The Gaussian predictor poorly represents the sharp edge at the north and south rim.}
    \end{subfigure}
    ~
        \begin{subfigure}[t]{0.23\textwidth}
        \centering
\includegraphics[trim= 20 0 50 20, clip, width=0.95\textwidth]{grandCanyon3D.eps}
        \caption{Predictor from Algorithm~\ref{TKL} trained on the elevation data in (b). The TK predictor accurately represent the north and south rims of the canyon.}
    \end{subfigure}
    \caption{Subfigure (a) shows an 3D representation of the section of the Grand Canyon to be fitted. In (b) we plot elevation data of this section of the Grand Canyon. In (c) we plot the predictor for a hand-tuned Gaussian kernel. In (d) we plot the predictor from Algorithm~\ref{TKL} where $d=2$. }\label{fig:GC} \vspace{-5mm}
\end{figure}
\vspace{-2mm}
\subsection{Solving $OPT\_A(P)$ } \vspace{-3mm}
For a given $P>0$, $OPT\_A(P)$ is a Quadratic Program (QP). General purpose QP solvers as applied to this problem have a worst-case complexity which scales as $O(m^3)$~\cite{ye1989extension} where $m$ is the number of data points. This computational complexity may be improved, however, by noting that the problem formulation is compatible with the representation defined in~\cite{LibSVM} for QPs derived from SVM. In this case, the algorithm in LibSVM~\cite{LibSVM} can reduce the computational burden somewhat. This improved performance is illustrated in Figure~\ref{fig:kernel_complexity} where we observe the achieved complexity scales as $O(m^{2.1})$. Note that for the 2-step algorithm proposed in this manuscript, solving the QP in $OPT\_A(P)$ is significantly slower that solving the Singular Value Decomposition (SVD) required for $OPT\_P(\alpha)$, which is defined in the following subsection. However, the achieved complexity of $O(m^{2.1})$ is also significantly faster than solving the large SDP, as described in~\cite{lanckriet_2004},~\cite{qiu2005multiple}, and~\cite{JMLR}. This complexity comparison will be further discussed in Section~\ref{sec:scalability}.\vspace{-3mm}

\subsection{Solving $OPT\_P(\alpha)$} \vspace{-3mm}
For a given $\alpha$, $OPT\_P(\alpha)$ is an SDP. Fortunately, however, this SDP is structured so as to admit an analytic solution using the SVD.
To solve $OPT\_P(\alpha)$ we minimize $O(\alpha,P)$ from Eq.~\eqref{OalphaP}
which, as per Corollary~\ref{corr}, is linear in $P$ and can be formulated as
\begin{align}  \label{opt:1}
OPT\_P(\alpha):=\min_{\substack{P \in \mathbb{R}^{q \times q}\\ \text{trace}(P) = q \\ P > 0}}  O(\alpha, P) := \min_{\substack{P \in \mathbb{R}^{q \times q}\\ \text{trace}(P) = q \\ P > 0}}  \text{trace}(C(\alpha)^T P)
\end{align}
where\vspace{-5mm}
\[C_{i,j}(\alpha) = \sum_{k,l = 1}^m (\alpha_k y_k) G_{i,j}(x_k,x_l) (\alpha_l y_l), \qquad G_{i,j}(x,y) := \begin{cases} g_{i,j}(x,y) & \text{if}~ i \leq \frac{q}{2}, j \leq \frac{q}{2} \\ t_{i,j}(x,y) & \text{if}~ i \leq \frac{q}{2}, j > \frac{q}{2} \\ t_{i,j}(y,x) & \text{if} ~i > \frac{q}{2}, j \leq \frac{q}{2} \\ h_{i,j}(x,y) & \text{if}~ i > \frac{q}{2}, j > \frac{q}{2} \end{cases}
\]
and $g$, $t$ and $h$ can be found in Corollary~\ref{corr}.

The following theorem gives an analytic solution for $OPT\_P$ using the SVD.
\begin{thm}
Let $C=V\Sigma V^T$ be the SVD of symmetric $C \in \mathbb{R}^{q \times q}$ and $v$ be the right singular vector corresponding to the minimum singular value of $C$. Then $P^* = q vv^T$ solves $OPT\_P$.
\end{thm} \vspace{-4mm}
\begin{proof}
Recall $OPT\_P$ has the form $\underset{P \in \mathbb{R}^{q \times q}}{\min}  \text{trace}(C^T P)~~\text{s.t.}~P\geq 0,~\text{trace}(P) = q$. \\
Denote the minimum singular value of $C$ as $\sigma_{\min}(C)$.  Then for any feasible $P\in \mcl X$, by~\cite{fang1994inequalities} we have
\[\text{trace}(C^TP) \geq \sigma_{\min}(C) trace(P) = \sigma_{\min}(C) q. \]
Now consider $P = q vv^T \in \mathbb{R}^{q \times q}$. $P$ is feasible since $P\ge 0$, and $\text{trace}(P) = q$. Furthermore,
\begin{align*}
\text{trace}(C P) ~=~q\,\text{trace}(V\Sigma V^T vv^T)~=~q\,\text{trace}(v^TV \Sigma V^T v) =~q\,\sigma_{\min}(C)
\end{align*}
as desired.\vspace{-2mm}
\end{proof}
Note that the size of the SVD problem in $OPT\_P(\alpha)$ is $q^2$, which increases with the number of features, which is typically relatively small. As a result, we observe that the $OPT\_P$ step of Algorithm~\ref{TKL} is typically less computationally intense than the $OPT\_A$ step.
\vspace{-3mm}

\begin{figure}[t]
    \centering
    \begin{subfigure}[t]{0.23\textwidth}
        \centering
\includegraphics[trim= 20 0 50 20, clip, width=0.8\textwidth]{complexityBinary.eps}
        \caption{Numerical complexity analysis of TKL for classification versus $m$.}
    \end{subfigure}%
    ~
    \begin{subfigure}[t]{0.23\textwidth}
        \centering
\includegraphics[trim= 20 0 50 20, clip, width=0.8\textwidth]{complexityRegression.eps}
        \caption{Numerical complexity analysis of TKL for regression versus $m$.}
    \end{subfigure}
        ~
    \begin{subfigure}[t]{0.23\textwidth}
        \centering
\includegraphics[trim= 20 0 50 20, clip, width=0.8\textwidth]{complexityQBinary.eps}
        \caption{Numerical complexity analysis of TKL for classification versus $q$.}
    \end{subfigure}
    ~
        \begin{subfigure}[t]{0.23\textwidth}
        \centering
\includegraphics[trim= 20 0 50 20, clip, width=0.8\textwidth]{complexityQRegression.eps}
        \caption{Numerical complexity analysis of TKL for regression versus $q$.}
    \end{subfigure}
    \caption{In (a) and (b) we plot log scale plots of the time taken to optimize TKL as the number of inputs change for $P \in \mathbb{R}^{q\times q}$.  The line of best linear fit is plotted as a dotted line. In (c) and (d) we plot log scale plots of the time taken to optimize TKL as the value of $q$ changes for four different values of $m$. }\label{fig:kernel_complexity} \vspace{-6mm}
\end{figure}

\section{Complexity and Scalability of the New TK Kernel Learning Algorithm}\label{sec:scalability} \vspace{-4mm}
We consider the computational complexity of Algorithm~\ref{TKL}.
If we define the number of data points used to learn the TK kernel function as $m$ and the size of $P$ as $q \times q$, then we find experimentally that the complexity of Algorithm~\ref{TKL} scales as approximately $O(m^{2.16}q^{2.23})$ for classification and $O(m^{2.24}q^{3.59})$ for regression as can be seen in Fig.~\ref{fig:kernel_complexity}. These results are lower with respect to $m$ than the value of $O(m^{2.6}q^{1.9})$ reported in~\cite{JMLR} for binary classification. The values for classification and regression are both estimated using the data set: Combined Cycle Power Plant (CCPP) in ~\cite{tufekci2014prediction,kaya2012local}, containing 4 features and $m = 9568$ samples. In the case of classification, labels with value greater than or equal to the median of the were relabeled as $1$, and those less than the median were relabeled as $-1$. Note that to study scalability in $q$, we varied the number of features in the dataset - thereby incrementing the size of the matrix $P\in\mathbb{R}^{q\times q}$.

Aside from improved scalability, the overall time required for Algorithm~\ref{TKL} is significantly reduced when compared with the algorithm in~\cite{JMLR}, improving by two orders of magnitude in some cases. This is illustrated for classification using four data sets in Table~\ref{time-table}. This improved complexity is likely due to the lower overhead associated with QP and the SVD.
\vspace{-3mm}

\begin{table*}[t]
  \caption{We report the mean computation time (in seconds), along with standard deviation, for 30 trials comparing the SDP algorithm in~\cite{JMLR} and the new TKL algorithm on several data sets.  All tests are run on a computer with an Intel i7-5960X CPU at 3.00 GHz with 128 Gb of RAM.}
  \label{time-table}
  \centering
\begin{tabular}{  c |  c c c c } \hline
Method & Liver \cite{UCI}   & Cancer \cite{mangasarian1990pattern}  &  Heart \cite{UCI} & Pima \cite{UCI} \\ \hline
SDP       &  95.75 $\pm$ 2.68 &  636.17 $\pm$ 25.43 & 221.67 $\pm$ 29.63 & 1211.66 $\pm$ 27.01 \\
TKL       &  1.10  $\pm$ 0.24 &  8.20 $\pm$ 0.36 & 3.35 $\pm$ 0.26 & 12.66 $\pm$ 0.44 \\ \hline
\end{tabular}
\end{table*}

\section{Accuracy of the New TK Kernel Learning Algorithm for Regression}\label{sec:8} \vspace{-4mm}
As expected, for classification, the accuracy of the new TK kernel learning algorithm (TKL) is identical to the analysis in~\cite{JMLR}.

For regression, we evaluate the accuracy of TKL when compared to other state of the art machine learning algorithms. Because the set of TK kernels is dense, for classification (as shown in~\cite{JMLR}), TKL outperforms all existing algorithm with respect to TSA. For regression, the appropriate metric is Mean Square Error (MSE). The algorithms used in our comparison are as follows.

\noindent\textbf{[TKL]} Algorithm~\ref{TKL} with $d=1$, $\epsilon = .1$ and we scale the data so that $x_i \in [0,1]^n$, and then select $[a,b] = [0-\delta,1+\delta]^n$, where $\delta>0$ and $C$ are chosen by 5-fold cross-validation;

\noindent\textbf{[SimpleMKL]} We use SimpleMKL \cite{rakotomamonjy_2008} with a standard selection of Gaussian and polynomial kernels with bandwidths arbitrarily chosen between .5 and 10 and polynomial degrees one through three - yielding approximately $13(n+1)$ kernels. We set $\epsilon = .1$ as in TKL and $C$ is chosen by 5-fold cross-validation;

\noindent\textbf{[Neural Net]} We use a 3 layer neural network with 50 hidden layers using MATLABs (\texttt{feedforwardnet}) implementation and stopped learning after the error in a validation set decreased sequentially 50 times.

In Table~\ref{table:Reg}, we see the average MSE on the test set for these three approaches as applied to randomly selected regression benchmark data sets where $n$ is the dimension of the data, $m$ is the number of training data and $m_t$ is the number of testing data points.  In all cases except Forest, [TKL] had both a lower (or comparable) computation time and MSE than both SimpleMKL and Neural Net.
In all cases, the MSE for TKL was significantly lower - illustrating the importance of the density property.

To further illustrate the importance of density property and the TKL framework for practical regression problems, we used elevation data from~\cite{becker2009global} to learn a TK kernel and associated SVM predictor representing the surface of the Grand Canyon in Arizona. This data set is particularly challenging due to the variety of geographical features. The result of the TKL algorithm can be seen in Figure~\ref{fig:GC}(d).
\vspace{-3mm}

\begin{table*}[t]
  \caption{Mean Square Error comparison for algorithms [TKL], [SimpleMKL] and [Neural Net].  In the data set column $m$ is the number of points in the training data set and $n$ is the number of features.  All tests are run on a computer with an Intel i7-5960X CPU at 3.00 GHz with 128 Gb of RAM.}
  \label{table:Reg}
  \centering \small
\begin{tabular}{  c |  c  c c || c |  c  c c } \hline
Data Set & Method  &   Error   & Time  & Data Set & Method        &   Error                & Time \\ \hline
CCPP \cite{tufekci2014prediction,kaya2012local}        &  TKL    &  9.70   & 1463.8  & Abalone \cite{UCI}             &  TKL    & 3.43    &  522.5  \\
$n$ = 4, $m$ = 8000  &  SimpleMKL      &    13.77      & 26097.1  & $n$ = 8, $m$ = 4000  &  SimpleMKL      &  4.28  &  1185.3 \\
$m_t$ = 1568       &  Neural Net     &    15.00       & 850.4 & $m_t$ = 177       &   Neural Net  &          8.72 & 483.4 \\ \hline
Airfoil  \cite{UCI}     &  TKL    &  1.46    &  92.1  & Forest  \cite{cortez2007data}            &  TKL    &  2.05    &  7.6  \\
$n$ = 5, $m$ = 1300 &  SimpleMKL      &    3.63      &  1025.0 & $n$ = 10, $m$ = 457  &  SimpleMKL      &    2.07      &  0.8 \\
$m_t$ = 203       &   Neural Net  &     4.28      & 61.3   & $m_t$ = 50       &   Neural Net  &     6.40      & 117.7 \\ \hline
\end{tabular}
\vspace{-5mm}
\end{table*}

\section{Conclusion} \vspace{-4mm}
We have extended the TK kernel learning framework to regression problems and proposed a faster algorithm for TK kernel learning which can be used for both classification and regression. The set of TK kernels is tractable, dense, and universal - implying that KL algorithms based on TK kernels are more robust - resulting in higher TSA for classification and lower MSE for regression. These three properties, combined with the improved computational complexity of the new algorithm, has resulted in a kernel learning framework which achieves both lower MSE and computation time when compared to both SimpleMKL and neural networks.

\newpage

\subsubsection*{Broader Impacts}
While machine learning algorithm have become very accurate in recent years, they perform poorly when faced with changes in the underlying process. As evidenced by Covid19, predictive models based on ML algorithms can be brittle~\cite{mims_2020}. The density property of the TK class ensures that the models generated using the algorithms described in this manuscript will be more robust to such changes in environment. Naturally, however, over-reliance on predictive models, without understanding of the process, can lead to negative outcomes, even if the models are robust.

\bibliographystyle{plain}


\section{Introduction}
\label{sec:introduction}
In recent years, vision-based crowd analysis has been extensively researched, due to its wide applications in crowd management, traffic control, urban planning, and surveillance. As one of the most important applications, crowd counting has been studied extensively \cite{Zhang_2016_CVPR_MCNN,Idrees_2018_ECCV_CL,Ma_2019_ICCV_BL}. With the recent progress of deep learning techniques, the performance of crowd counting has been significantly elevated \cite{Li_2018_CVPR_CSRNet,Ma_2019_ICCV_BL,Yan_2019_ICCV_PGCNet,Liu_2019_ICCV_DSSINet}. The convolutional neural network (CNN)-based methods have demonstrated excellent performance on the task of counting dense crowds in images.
Most of CNN-based methods first estimate the density map via deep neural networks and then calculate the counts \cite{Zhang_2016_CVPR_MCNN,Li_2018_CVPR_CSRNet,Cao_2018_ECCV_SA,Wan_2019_ICCV_adaptive_map,Yan_2019_ICCV_PGCNet}. In specific, the concept of {density map}, where the integral (sum) over any sub-region equals the number of objects in that region, was first proposed in \cite{NIPS2010_4043_map}.
Since the existing crowd counting benchmarks provide the point annotation for each crowd image, in which each point is located on the head of a person in the crowd. To train a CNN-based crowd counting model, the point annotations need to be converted to a density map in advance.
Lempitsky \textit{et al.} \cite{NIPS2010_4043_map} propose to use a normalized 2D Gaussian kernel to convert the point annotations to a ground-truth density map. Typically, the Gaussian kernel size is fixed while converting point annotations. But this trivial approach degrades the counting performance, since the scales of individuals in the crowd image may vary greatly. To produce better ground-truth density map, Zhang \textit{et al.} \cite{zhang2015cross} provides a manual estimation of the perspective maps for the crowd images. But it is laborious to provide the accurate perspective information for all the image captured in various scenarios. Zhang \textit{et al.} \cite{Zhang_2016_CVPR_MCNN} introduce the geometry-adaptive kernel to create ground-truth density maps. They assume that the crowd is evenly distributed and thus they can estimate the kernel size by the average distance between each point annotations and its nearest neighbors.
Generally, geometry-adaptive kernels is efficient for estimating the kernel of dense point annotations, but it is inaccurate when the crowd is not distributed evenly.
As we know, with different Gaussian kernel sizes, the ground-truth density maps converted from the point annotations can be sharper with smaller kernel sizes or smoother with larger kernel sizes (as shown in Fig.~\ref{fig:density_map}(b) and (c)). Empirically, it is easier for a CNN model to fit the smoother crowd density maps rather than the sharper ones. This is probably because the sharper ground-truth density maps contain a larger amount of zero values than the non-zero ones, which makes the network difficult to fit the ground-truth.
Nevertheless, the individual information contained in the smoother density maps is relatively vague compared with the sharper ones. So the performance of the crowd counting model will thus be degraded and the prior density map generation approaches can hardly handle this problem.

To mitigate this problem, our work attempts to tackle it in several aspects.
First, we propose a parametric perspective-aware density map generation approach. We assume that the individuals on the same horizontal line (or the same row) of the image are from the similar distance away from the camera, so the kernel sizes for the point annotations on the same row should be the same.
Thus, to determine the kernel size for each row of the image, we introduce a metric, \textit{effective density}, used to measure the density of the highly aggregated segments on each row. By linearly mapping the effective densities to the manually-defined kernel size range, we can easily produce perspective-aware density maps by assigning a small number of parameters as the ground-truths of model training.
Second, we propose a simple CNN-based architecture featuring with two output branches which are supervised by a multi-task loss with low-resolution and high-resolution ground-truth density maps. With the supervision of high- and low-resolution density maps, our model is able to generate the crowd distribution with the relatively high dimension (i.e., $1/4$ of the input size).

Last but not least, we propose an iterative distillation optimization algorithm for progressively enhancing the performance of our network model. As mentioned, although our network can regress the crowd distribution, its performance is still constrained by the quality of the density maps. To benefit learning from more accurate yet hard-to-learn density maps (i.e., the sharp density map generated by small Gaussian kernel), we propose to iteratively distill the network with the previously trained identical network. As known, the distillation techniques have been previously used to compress network or improve network capability \cite{hinton2015distilling,yang2019snapshot,furlanello2018born}. Here, we employ it to strengthen the capability of our network, when the training objective of the crowd counting model becomes more and more challenging. Particularly, during distillation, our proposed parametric density map generation approach can be iteratively utilized to generate sharper ground-truth density maps as the training objectives. In experiments, we demonstrate that our perspective-aware density map generation method is better than the prior generation techniques. Besides, we show that, although our network architecture is simple, our approach can still obtain the state-of-the-art performance compared with other methods.

In summary, our contributions are below:
\begin{itemize}
	\item We introduce a parameteric perspective-aware density map generation method to generate ground-truth density maps from crowd point annotations, so as to train a crowd counting model that can estimate the crowd density maps with relatively high spatial dimension.
	\item We present a novel iterative distillation algorithm to enhance our model while progressively reducing the Gaussian kernel sizes of the ground-truth density maps, which can further improve the counting performance.
	\item In experiments, we show that our simple network architecture can reach the state-of-the-art performance compared against the latest approaches in public benchmarks.
\end{itemize}

In the following, we first review the related works in Sec.~\ref{sec:rw}. Then, we elaborate the methodology in Sec.~\ref{sec:method}. Lastly, in Sec.~\ref{sec:exp}, we evaluate our proposed approach in public benchmarks.

\section{Related work}
\label{sec:rw}

In this section, we survey the most related works on density map generation and knowledge distillation.

\subsection{Density map generation for crowd counting}

The perspective distortion in images is one of the main factors that affect the accurate counting of the crowd. Due to the shooting scene, angle, terrain distribution and other factors, the perspective of each picture is different. It is a natural idea to combine perspective information to improve the ground-truth density map. By manually measuring the height of people in different positions in the picture, \cite{zhang2015cross} calculated the perspective of 108 scenes in WorldExpo'10 dataset , and used the perspective to generate an appropriate ground\_truth density map for the dataset.
However, most datasets contain various scenes (such as ShangHaiTech~\cite{Zhang_2016_CVPR_MCNN}, UCF\_QNRF~\cite{Idrees_2018_ECCV_CL}), so manually estimating the perspective of each image is very laborious. Thus, Zhang \textit{et al.} \cite{Zhang_2016_CVPR_MCNN} proposed the adaptive geometric estimation technique. The Gaussian covariance of each marker point is estimated by calculating its average distance from the surrounding points, thereby generating a ground-truth density map that more closely matches the actual distribution of the crowd.
Shi \textit{et al.} \cite{shi2019counting} proposed a non-uniform kernel estimation. They assumes that the crowd is unevenly distributed throughout the crowd image.
They first used {adaptive geometric estimation technique} to estimate the covariance of all points,
and then calculated the average covariance of the local neighborhood of each point, which can thus make the Gaussian distribution of a single point will not become too large or too small.
Compared with these methods, we propose a new method to produce a ground-truth density map, by assuming the points on each row of the image have the Gaussian kernels with the similar size and estimating the perspective of the entire image via rough density maps without any supervision.

\subsection{Knowledge distillation}
Knowlege distillation is a deep learning training technique that trains a network model (i.e., the student model) to mimic the behavior of a independently trained model (i.e., the teacher model). It was originally used for model compression when the student model is light-weighted. However, some recent studies~\cite{furlanello2018born,yang2019snapshot} have found that the distillation between the models with the identical network structure can obtain better classification results.
Particularly, Yang \textit{et al.} \cite{yang2019snapshot} proposed the snapshot distillation method, which divides the training process into two stages. In the first stage, the network performs normal training. Then use the first-stage model as a teacher model to guide the student model for the second-stage training. Inspired by these works, we propose our distillation algorithm for the task of crowd counting.

\section{Our Proposed Method}
\label{sec:method}

In this section, we first introduce our deep network architecture. To progressively strengthen the performance of our network, we propose a distillation based optimization approach as well as our parametric density map generation method.

\subsection{Network architecture}

In order to train a crowd counting model, we follow the principle of the prior frameworks (e.g.~\cite{zhang2015cross}), in which the deep neural network model aims to generate a density map from the input crowd image and its crowd count can be measured by accumulating the values of the entire density map.

As shown in Fig.~\ref{fig:network}, we adopt a vanilla deep convolutional neural network (CNN) as the backbone network to estimate the count of a crowd image, which consists of a downsampling module and an upsampling module. The downsampling module adopts the first 10 convolutional layers of the pretrained model VGG-16 \cite{simonyan2014very}, which extracts deep visual feature representation from the input crowd image. Then, two transposed convolutional layers are applied to upsample the spatial dimension of the feature maps to $1/4$ of the input dimension. As depicted in Table~\ref{tab:resolution}, most prior counting models produce low-resolution density map (i.e., less than 1/16 of the original size), yet producing a high-resolution density map benefits many downstream applications such as crowd analysis.
The reason we adopt such a simple network architecture is to demonstrate the performance of our proposed density map generation and distillation method.

\begin{figure}
	\centering
	\includegraphics[width=1\linewidth]{images/framework.png}
	\caption{Illustration of our crowd counting network model. As shown, our framework is composed of feature extraction, upsampling, and multi-task loss. Particularly, the feature representation of the input crowd image can be extracted via pretrained convolutional layers. Through upsampling layers, our network reconstructs the extracted features to density maps. We introduce a multi-task loss to simultaneously supervise the high-resolution and low-resolution density map generation.
		The high resolution ground-truth is $1/4$ of the input image, while the low resolution one is $1/16$ of the input image. }
	\label{fig:network}

\end{figure}

\begin{table}[t]
	\setlength{\tabcolsep}{6pt}
	\centering
	\caption{The ratio of the dimensions of the output density maps and the input image for different crowd counting methods. }
	\begin{tabular}{c|c|c}
		\toprule
		Method         & Year        & Ratio  \\
		\midrule
		MBTTBF-SCFB \cite{Sindagi_2019_ICCV_MBTTBF}     & 2019 &1/256  \\
		CSRNet \cite{Li_2018_CVPR_CSRNet}                  &2018  &1/64  \\
		ADCrowdNet \cite{Liu_2019_CVPR_ADCrowd}  & 2019 &1/64  \\
		CSRNet+PACNN \cite{Shi_2019_CVPR_PACNN} & 2019&1/64  \\
		CAN \cite{Liu_2019_CVPR_CAN}             &2019 &1/64  \\
		BL \cite{Ma_2019_ICCV_BL}                    & 2019    &1/64  \\

		PACNN \cite{Shi_2019_CVPR_PACNN}       & 2019 &1/64  \\

		MCNN \cite{Zhang_2016_CVPR_MCNN}       & 2016 &1/16   \\

		Ours                   & - &1/4  \\
		\bottomrule
	\end{tabular}

	\label{tab:resolution}

\end{table}

However, it is challenging for a vanilla deep neural network to produce a high-resolution density map.
To make the network prediction robust, we introduce a multi-task loss by employing two separate upsampling branches for low-resolution and high-resolution supervision. As illustrated in Fig.~\ref{fig:network}, the losses of these two branches are supervised by the ground-truth density maps of two different scales, respectively. In particular, the low-resolution branch consists of two convolutional layers.
Hence, the network $\mathcal{F}$ is trained via a standard $L_2$ loss:
\begin{align}
	\mathcal{L}= \frac{1}{N} \sum_{i=1}^N &(\|\mathcal{F}_{HR}(X_i)-G_{HR}\|^2_2\nonumber\\
	&+ \|\mathcal{F}_{LR}(X_i)-G_{LR}\|^2_2),\label{eq:l2loss}
\end{align}
where $X_i$ denotes the $i^{th}$ input crowd image in the $N$ training images. $\mathcal{F}_{HR}$ and $\mathcal{F}_{LR}$ are the high-resolution and low-resolution outputs of our network. $G_{HR}$ and $G_{LR}$ represent the high-resolution and low-resolution ground-truth density maps, respectively.

\subsection{Parametric perspective-aware density map}

The supervision of our crowd counting network requires the ground-truth crowd density map.
Existing crowd counting benchmarks provide the point annotation for each crowd image, in which each point annotation represents the position of a person in the crowd. To obtain the crowd density map, following \cite{NIPS2010_4043_map}, prior crowd counting works convert the point annotation of a crowd image into a density map by applying a Gaussian kernel over each point.
Most prior methods apply the Gaussian kernel with a fixed kernel size \cite{NIPS2010_4043_map,Cao_2018_ECCV_SA}, or the geometry-adaptive kernel size \cite{8099912_switch,Li_2018_CVPR_CSRNet,Liu_2019_CVPR_CAN,Liu_2019_ICCV_DSSINet,Ma_2019_ICCV_BL,xu2019learn,xiong2019open,shen2018crowd}. These density conversion techniques do not consider the perspective information or non-uniform crowd distribution of the crowd images.
Based on the recent findings \cite{Wan_2019_ICCV_adaptive_map}, the generated density map may affect the performance of the model, so these trival density map generation methods may not help us achieve the satisfactory results. Thus, we propose a spatial-aware parametric method for generating adaptive density maps.

Our proposed density conversion approach is based a simple assumption that the majority of persons located on the same row of the image have a similar distance away from the camera (e.g. Fig.~\ref{fig:density_map}(a)).
This assumption may not be suited for certain crowd scenes, but it empirically works well for most scenes. Thus, subject to the assumption, each row of the density map should be corresponded to a Gaussian kernel with the same size, which can approximate the scale variance caused by the perspective view of the crowd image.

To estimate the Gaussian kernel size for each row, we first apply a fixed large Gaussian kernel size (e.g., the standard deviation of Gaussian kernel $\sigma=25$) to produce a rough yet smooth density map $\bar{G}$ as a prior (see Fig.~\ref{fig:density_map}(b)). Due to the perspective effect of the image, the density concerns with not only the spatial distribution but also the distance away from camera. For instance, given two persons sitting side-by-side, they appear to be closer when they are farther away from camera and vice versa. Thus, we can approximately find out the relation between the density and the distance from camera for each row.
We propose to calculate the \textit{effective density} for the $i$-th row of $\bar{G}$:
\begin{align}
	D(i) &= \frac{1}{M}\sum_j \delta[\bar{G}(i,j) > \epsilon] \cdot \bar{G}(i,j),\\
	&\text{s.t. } M=\sum_j \delta[\bar{G}(i,j) > \epsilon],
\end{align}
where $\bar{G}(i,j)$ refers to the density value at the location $(i,j)$ of the image. $\delta[\cdot]$ denotes Dirac delta function, which equals to $1$ when the condition in the bracket is satisfied and $0$ otherwise. $\epsilon$ refers to the manually defined threshold.

\begin{figure}
	\centering
	\setlength{\tabcolsep}{1pt}
	\begin{tabular}{cccc}
		\includegraphics[width=0.30\linewidth]{images/img.jpg}&
		\includegraphics[width=0.30\linewidth]{images/k1.jpg}&
		\includegraphics[width=0.30\linewidth]{images/newmap.jpg}&	\includegraphics[trim=2.71cm 0 2.71cm 0, clip, width=0.06\linewidth]{images/k2.jpg}\\
		(a) & (b) & (c) &
	\end{tabular}
	\caption{(a) The input crowd image; (b) The prior density map $\bar{G}$ generated by a fixed and large Gaussian kernel size; (c) The density map generated by our proposed method. The rightmost color bar indicates the kernel size for each row.}
	\label{fig:density_map}

\end{figure}

In essence, the \textit{effective density} measures the average density of most dense segment on each row of the crowd image while filtering out the less dense part. Thus, we can determine the kernel size for each row according to its $D(\cdot)$. To accomplish this, we apply a linear mapping between the maximum and minimum values of $D(\cdot)$ (i.e. $D_{max}$ and $D_{min}$) and
the largest and smallest kernel size (i.e., the standard deviations $\sigma_{max}$ and $\sigma_{min}$), respectively. Specifically, $\sigma_{max}$ and $\sigma_{min}$ are manually determined, while $D_{max}$ and $D_{min}$ are measured over the entire training set.
Hence, the Gaussian kernel size of the $i^{th}$ row can be computed as:
\begin{align}
	\sigma(i) = \alpha D^{-1}(i) + \beta,
\end{align}
where
\begin{align}
	\beta &= \frac{D^{-1}_{min}\sigma_{min}-D^{-1}_{max}\sigma_{max}}{D^{-1}_{min}-D^{-1}_{max}},\\
	\alpha&= \frac{\sigma_{max}-\sigma_{min}}{D^{-1}_{min} - D^{-1}_{max}}.
\end{align}

Thus, after assigning the kernel size for each row, we can obtain the density map determined by $\sigma_{max}$ and $\sigma_{min}$, denoted as $G(\sigma_{min},\sigma_{max})$. We illustrate an example of our parametric method in Fig.~\ref{fig:density_map}(c).
In practice, $\sigma_{min}$ is often fixed at a small value (e.g. $2.5$), while we mainly tune the value of $\sigma_{max}$ to achieve various density maps. Therefore, for the sake of simplicity, the notation of $G(\sigma_{min},\sigma_{max})$ can be simplified as $G(\sigma_{max})$.
In the extreme case where $\sigma_{max}$ equals to $\sigma_{min}$, the generated density map is the same as the one generated by a universally-fixed kernel size.

As described in the following section, our proposed density map generation technique can benefit the training of crowd counting model by providing a perspective-aware and parametric density map.

\subsection{Distillation based crowd counting}
\label{Sec:distillation}

As shown in Fig.~\ref{fig:density_map}(b) and (c), with a larger kernel size (i.e. larger $\sigma_{max}$ for $G(\sigma_{max})$), the density map will appear to be smoother, and vice versa.
On one hand, for a deep neural network that regresses the density maps, the smoother density maps are easier to learn than the sharper ones.
On the other hand, training the regression network from such smoother density maps may degrade the results of crowd counting, since the smoother density maps blur the individual information of the crowd image. In practice, proper parameters are often empirically chosen to trade off the counting accuracy and model training.

\subsubsection{Crowd counting network distillation}

In order to break the bond caused by density maps, we propose a distillation based optimization method that enables to progressively reduce the difficulty of learning from sharper density maps so as to obtain a better solution. Knowledge distillation has been proposed in \cite{hinton2015distilling}, which has been applied for network compression. Specificially, a lightweight model is often applied to learn the behavior of a large network. The recent studies find that the identical network structures can benefit from the distillation \cite{furlanello2018born,yang2019snapshot}. Here, we employ the similar optimization strategy to iteratively train our network.

First, we train our network according to Eq.~\ref{eq:l2loss} that can be simply expressed as:
\begin{align}
	\mathcal{L}^{(t=0)}=  \frac{1}{N}{\sum_{i=1}^N\|\mathcal{F}(X_i)-G\|^2}.\label{eq:init}
\end{align}
Next, we treat the trained network $\mathcal{F}$ as the teacher model (denoted as $\mathcal{F}_T$) and we leverage it to train an identical network (i.e., the student model $\mathcal{F}_S$) from scratch. Thus, the trained student model at this stage can be treated as the teacher model of the next stage, i.e. $\mathcal{F}^{(t)}_T \leftarrow \mathcal{F}^{(t-1)}_S$, where $t$ denotes the timestamp of the training stage. Hence, the training loss of the $t^{th}$ stage in distillation can be expressed as follows:
\begin{align}
	\mathcal{L}^{(t)}= & \frac{1}{N}{\sum_{i=1}^N\|\mathcal{F}^{(t)}_S(X_i)-G\|^2} + \nonumber\\
	&  \frac{\lambda}{N} {\sum_{i=1}^N\|\mathcal{F}^{(t)}_S(X_i)-\mathcal{F}^{(t)}_T(X_i)\|^2},\label{eq:distill}
\end{align}
where $\lambda$ denotes the constant balance weight. The first term is a standard $L2$ loss, while the second term aims to align the outputs of the teacher model and the student model, which forces the student model to approach the behavior of the teacher model.

\begin{figure}
	\centering
	\includegraphics[width=0.5\textwidth]{images/distillation.png}
	\caption{Illustration of the distillation-based optimization for our crowd counting model. As shown, we first train a teacher model and freeze its weights. Then, we utilize the teacher model to train another model (i.e., the student model) with identical structure supervised by the low-resolution and high-resolution ground-truth density maps.}
	\label{fig:distill}

\end{figure}

The distillation based optimization can progressively optimize the network. However, as we mentioned, the inappropriately computed ground-truth density maps will degrade the counting performance. Furthermore, the static ground-truth density maps (i.e., $G$ remains the same in distillation) limit the optimal performance discovered by the distillation. Thus, in the following section, we introduce a simple yet effective method in order to produce adaptive density maps.

\subsubsection{Density-aware distillation}

To further strengthen the network performance, utilizing our parameteric perspective-aware density map generation technique, we improve the distillation method. According to Eq.~\ref{eq:distill}, the student model $\mathcal{F}_S$ will learn the behavior of the teacher model $\mathcal{F}_T$. But, with the static density map $G$ as ground-truth, the distillation may not further improve the performance. However, it is difficult to directly train a network from sharp density maps generated with a small kernel size, which often leads to a poor performance. With distillation, the model is able to progressively adapt to the sharper and sharper density map.
Particularly, in each stage of distillation, we slightly increase the difficulty of the task by introducing the density maps with a smaller kernel size.
In specific, on every stage of distillation, the training loss is modified as below:
\begin{align}
	\mathcal{L}^{(t)}= & \frac{1}{N}{\sum_{i=1}^N\|\mathcal{F}^{(t)}_S(X_i)-G^{(t)}(\sigma^{(t)}_{max})\|^2} \nonumber\\
	&  +\frac{\lambda}{N} {\sum_{i=1}^N\|\mathcal{F}^{(t)}_S(X_i)-\mathcal{F}^{(t-1)}_T(X_i)\|^2},\label{eq:tsoptim}
\end{align}
where $\sigma^{(t)}_{max}$ refers to the maximal Gaussian kernel size at the $t^{th}$ training stage. After each stage, $\sigma^{(t)}_{max}$ will be updated, i.e., $\sigma^{(t)}_{max}=w \sigma^{(t-1)}_{max}$, where $w$ is a constant within the value range of $[0,1)$. Our algorithm is summarized in Algorithm~\ref{algo}. By iteratively distilling the model while reducing the kernel size to generate sharper density maps for network to learn, the performance of the student model can iteratively be strengthened.

\begin{algorithm}
	\SetAlgoLined

	\KwIn{$\text{Network }\mathcal{F}, \text{Input images }X_i (i=1\cdots N)$, $\text{Maximal kernel size }\sigma^{(0)}_{max}$}
	Compute density map $G$ using $\sigma^{(0)}_{max}$ (i.e. $G^{(t=0)}$)\;
	Train the initial model $\mathcal{F}$ supervised by $G$ (Eq.~\ref{eq:init})\;
	Assign the teacher model $\mathcal{F}_T^{(t=0)} \leftarrow \mathcal{F}$\;
	\For{t=1:T}{
		$\sigma^{(t)}_{max} \leftarrow w\sigma^{(t-1)}_{max}$\;
		Produce density map $G^{(t)}$ using $\sigma^{(t)}_{max}$\;
		Set up a new student model $\mathcal{F}_S^{(t)}$\;
		Train $\mathcal{F}_S^{(t)}$ using $G^{(t)}$ and $\mathcal{F}_T^{(t-1)}$ (Eq.~\ref{eq:tsoptim})\;
		$\mathcal{F}_T^{(t)} \leftarrow \mathcal{F}_S^{(t)}$\;
	}
	\caption{Iterative distillation of crowd counting}
	\label{algo}
\end{algorithm}

\section{Experimental Results}
\label{sec:exp}

In this section, we first introduce the datasets used for evaluation and the metrics, as well as the implementation details. Then, we conduct the comparison experiments with the state-of-the-art methods in public benchmarks. Last, we perform the ablation study to investigate the density map generation, the distillation algorithm, and the multi-task loss of our model.

\subsection{Datasets and evaluation metrics}

We evaluate our approach on public crowd counting benchmarks: ShanghaiTech Part A/B~\cite{Zhang_2016_CVPR_MCNN}, UCF-QNRF~\cite{Idrees_2018_ECCV_CL}, and UCSD~\cite{chan2008privacy}.
Particularly, ShanghaiTech dataset has a total of 1198 crowd images, including Part A and Part B. Part A contains 482 images, and Part B contains 716 images. The QNRF dataset has 1535 images with average resolution at $2013\times 2902$ and a total of 1.25 million annotations, in which there are 1201 images for training and 334 images for testing. Besides, the UCSD dataset have 2000 images with the resolution of $238\times 158$ with relatively smaller density.

Following prior works, MSE and MAE are used as metrics, which are defined as follows:
\begin{align}
	MAE &= \frac{1}{N} \sum_{i=1}^N\|E_i-GT_i\|_1,\\
	MSE &= \frac{1}{N} \sqrt{\sum_{i=1}^N\|E_i-GT_i\|_2^2},
\end{align}
where $N$ is the number of test images, $E_i$ indicates the estimated density map count of the $i$-th image and $GT_i$ indicates the ground-truth count of the $i$-th image.

\subsection{Implementation details}

In practice, we set $\sigma_{min}$ as 2.5 and the initial $\sigma_{max}$ as 25. In the first timestep, the $\sigma_{max}$ is set as 20 and $w$ as 0.5.
On the initialization stage, the learning rate is set as $5\times 10^{-6}$ and the momentum 0.95. Since we do not normalize the input image, the batch size is set as 1. Other training settings differ across the benchmarks, so we will elaborate them below, respectively.

\begin{table}
	\setlength{\tabcolsep}{4pt}
	\centering

	\caption{The comparison results on the ShanghaiTech Part A (SHA) and Part B (SHB). }
	\begin{tabular}{c|c|c|c|c|c}
		\toprule
		\multirow{2}{*}{Method}  & \multirow{2}{*}{Year}& \multicolumn{2}{c|}{SHA}& \multicolumn{2}{c}{SHB} \\

		\cmidrule{3-6}
		& &  MAE& MSE & MAE & MSE \\
		\midrule
		\midrule
		MCNN~\cite{Zhang_2016_CVPR_MCNN}  & 2016 & 110.2 & 173.2 & 26.4 & 41.3\\
		Switching CNN~\cite{8099912_switch} &2017 & 90.4 & 135.0 & 21.6& 33.4 \\
		SANet~\cite{Cao_2018_ECCV_SA} & 2018& 67.0 & 104.5 & 8.4& 13.6\\
		CSRNet~\cite{Li_2018_CVPR_CSRNet}    &2018& 68.2 & 115.0 & 10.6 & 16.0 \\
		ic-CNN~\cite{ranjan2018iterative} & 2018 & 69.8& 117.3& 10.4 & 16.7\\
		CSRNet+PACNN~\cite{Shi_2019_CVPR_PACNN}   &2019             &62.4  &102.0 & 7.6 & 11.8\\
		ADCrowdNet~\cite{Liu_2019_CVPR_ADCrowd} &2019&63.2  &98.9 & 7.6 & 13.9 \\
		PACNN~\cite{Shi_2019_CVPR_PACNN} & 2019& 66.3 & 106.4 & 8.9& 13.5\\
		CAN~\cite{Liu_2019_CVPR_CAN}   & 2019             &62.3  & 100.0 &7.8 &12.2 \\
		BL~\cite{Ma_2019_ICCV_BL}            & 2019              &62.8  & 101.8 &7.7&12.7 \\
		TEDnet~\cite{jiang2019crowd} & 2019 & 64.2 & 109.1 & 8.2 & 12.8\\
		HA-CCN~\cite{sindagi2019ha} & 2019& 62.9 & \textbf{94.9} & 8.1 & 13.4\\

		Ours                       & - &\textbf{61.1} &104.7 &\textbf{7.5}&\textbf{12.0} \\
		\bottomrule

	\end{tabular}

	\label{Tab:Shanghaitech}

\end{table}

\begin{table}
	\setlength{\tabcolsep}{6pt}
	\centering

	\caption{Comparison of our approach with other state-of-the-art methods on UCF-QNRF.}
	\begin{tabular}{c|c|c|c}

		\toprule
		Method & Year & MAE& MSE\\
		\midrule
		\midrule
		MCNN~\cite{Zhang_2016_CVPR_MCNN} & 2016 & 277.0 & 426.0\\
		Switching CNN~\cite{8099912_switch} & 2017 & 228.0 & 445.0\\
		CL~\cite{Idrees_2018_ECCV_CL}           & 2018    &132.0   &191.0 \\
		HA-CCN~\cite{sindagi2019ha} & 2019 & 118.1 & 180.4\\
		RANet~\cite{zhang2019relational} & 2019 & 111.0 & 190.0\\
		CAN~\cite{Liu_2019_CVPR_CAN}             & 2019   &107.0   &183.0 \\
		TEDnet~\cite{jiang2019crowd} & 2019 & 113.0 & 188.0\\
		SPN+L2SM~\cite{xu2019learn} & 2019 & 104.7 & 173.6\\
		S-DCNet~\cite{xiong2019open} &	2019 & 104.4	&	176.1\\
		SFCN~\cite{wang2019learning}	& 2019 & 102.0	&	171.4\\
		DSSINet~\cite{Liu_2019_ICCV_DSSINet}     & 2019       &99.1  &\textbf{159.2}  \\
		MBTTBF-SCFB~\cite{Sindagi_2019_ICCV_MBTTBF}    & 2019    &97.5  &165.2 \\

		Ours             & -           &\textbf{92.9} &\textbf{159.2}  \\
		\bottomrule
	\end{tabular}

	\label{Tab:UCF-QNRF}

\end{table}

\begin{table}
	\setlength{\tabcolsep}{6pt}
	\centering

	\caption{The comparison results on UCSD. }
	\begin{tabular}{c|c|c|c}
		\toprule
		Method & Year & MAE& MSE\\
		\midrule
		\midrule
		MCNN~\cite{Zhang_2016_CVPR_MCNN} & 2016 & 1.07 & 1.35\\
		Switching CNN~\cite{8099912_switch} & 2017 & 1.62 & 2.10\\
		ConvLSTM~\cite{xiong2017spatiotemporal} & 2017 & 1.30 & 1.79\\
		BSAD~\cite{huang2017body} & 2017 & 1.00 & 1.40\\
		ACSCP~\cite{shen2018crowd} & 2018 & 1.04 & 1.35\\
		CSRNet~\cite{Li_2018_CVPR_CSRNet}   & 2018 &1.16  &1.47  \\
		SANet~\cite{Cao_2018_ECCV_SA}             & 2018         &1.02  &1.29  \\  SANet+SPANet~\cite{Cheng_2019_ICCV_SPANet}  & 2019              &1.00  &1.28 \\
		ADCrowdNet~\cite{Liu_2019_CVPR_ADCrowd}         & 2019    &\textbf{0.98}  &1.25 \\

		Ours         &     -          &\textbf{0.98}  &  \textbf{1.24}\\
		\bottomrule
	\end{tabular}

	\label{Tab:UCSD}
	\vspace{-0.4cm}

\end{table}

\begin{figure*}
	\centering
	\setlength{\tabcolsep}{1pt}
	\begin{tabular}{ccccc}
		Input image & Bayesian~\cite{Ma_2019_ICCV_BL} & CSRNet~\cite{Li_2018_CVPR_CSRNet} & Ours & Ground-truth\\
		\includegraphics[width=0.2\textwidth]{images/dif_method/IMG_20.jpg}&
		\includegraphics[width=0.2\textwidth]{images/dif_method/bayesianIMG_20_pre_464_4188232421875.jpg}&
		\includegraphics[width=0.2\textwidth]{images/dif_method/csrnetIMG_20_pre_481_46649169921875.jpg}&
		\includegraphics[width=0.2\textwidth]{images/dif_method/d4_IMG_20_pre_460_40301513671875.jpg}&
		\includegraphics[width=0.2\textwidth]{images/dif_method/d4_IMG_20_gt_460_0.jpg}\\
		Estimated count & 464.4 & 481.5 & 460.4 & 460\\
		\includegraphics[width=0.2\textwidth]{images/dif_method/IMG_22.jpg}&
		\includegraphics[width=0.2\textwidth]{images/dif_method/bayesianIMG_22_pre_269_8780517578125.jpg}&
		\includegraphics[width=0.2\textwidth]{images/dif_method/csrnetIMG_22_pre_274_3505859375.jpg}&
		\includegraphics[width=0.2\textwidth]{images/dif_method/d4_IMG_22_pre_231_74530029296875.jpg}&
		\includegraphics[width=0.2\textwidth]{images/dif_method/d4_IMG_22_gt_211_99996948242188.jpg}\\
		Estimated count & 269.9 & 274.4 & 231.7 & 212\\
		\includegraphics[width=0.2\textwidth]{images/dif_method/IMG_27.jpg}&
		\includegraphics[width=0.2\textwidth]{images/dif_method/bayesianIMG_27_pre_797_762939453125.jpg}&
		\includegraphics[width=0.2\textwidth]{images/dif_method/csrnetIMG_27_pre_917_490966796875.jpg}&
		\includegraphics[width=0.2\textwidth]{images/dif_method/d4_IMG_27_pre_758_097900390625.jpg}&
		\includegraphics[width=0.2\textwidth]{images/dif_method/d4_IMG_27_gt_759_9999389648438.jpg}\\
		Estimated count & 797.8 & 917.5 & 758.1 & 760\\
		\includegraphics[width=0.2\textwidth]{images/dif_method/IMG_28.jpg}&
		\includegraphics[width=0.2\textwidth]{images/dif_method/bayesianIMG_28_pre_542_6561279296875.jpg}&
		\includegraphics[width=0.2\textwidth]{images/dif_method/csrnetIMG_28_pre_556_3040771484375.jpg}&
		\includegraphics[width=0.2\textwidth]{images/dif_method/d4_IMG_28_pre_440_2284851074219.jpg}&
		\includegraphics[width=0.2\textwidth]{images/dif_method/d4_IMG_28_gt_422_9998474121094.jpg}\\
		Estimated count & 542.7 & 556.3 & 440.2 & 423\\
		\includegraphics[width=0.2\textwidth]{images/dif_method/IMG_31.jpg}&
		\includegraphics[width=0.2\textwidth]{images/dif_method/bayesianIMG_31_pre_260_5260314941406.jpg}&
		\includegraphics[width=0.2\textwidth]{images/dif_method/csrnetIMG_31_pre_255_705078125.jpg}&
		\includegraphics[width=0.2\textwidth]{images/dif_method/d4_IMG_31_pre_252_24856567382812.jpg}&
		\includegraphics[width=0.2\textwidth]{images/dif_method/d4_IMG_31_gt_249_99996948242188.jpg}\\
		Estimated count & 260.5 & 255.7 & 252.2 & 250\\
		\includegraphics[width=0.2\textwidth]{images/dif_method/IMG_46.jpg}&
		\includegraphics[width=0.2\textwidth]{images/dif_method/bayesianIMG_46_pre_404_7515869140625.jpg}&
		\includegraphics[width=0.2\textwidth]{images/dif_method/csrnetIMG_46_pre_466_1134033203125.jpg}&
		\includegraphics[width=0.2\textwidth]{images/dif_method/d4_IMG_46_pre_389_73870849609375.jpg}&
		\includegraphics[width=0.2\textwidth]{images/dif_method/d4_IMG_46_gt_380_99993896484375.jpg}\\
		Estimated count & 404.8 & 466.1 & 389.7 & 381
	\end{tabular}
	\caption{Illustration of comparison examples against two of the state-of-the-art methods \cite{Li_2018_CVPR_CSRNet,Ma_2019_ICCV_BL}. }
	\label{fig:comparison}
\end{figure*}

\subsection{Comparison results}

In this section, we illustrate the quantitative results on different benchmarks including ShanghaiTech Part A and Part B~\cite{Zhang_2016_CVPR_MCNN}, UCF-QNRF~\cite{Idrees_2018_ECCV_CL}, and UCSD~\cite{chan2008privacy}. We demonstrate that the counting performance achieved by our simple network architecture can be better or comparable to the state-of-the-art models.

\textbf{ShanghaiTech Part A/B.}
To train our model, we follow the general training protocal. When training the initial counting model, we train the model for 100 epochs, and the learning rate is decayed by 10 times after training for 60, 80, and 90 epochs, respectively. During distillation, the learning rate will be reduced by 10 times after training for 20, 60, and 80 epochs.
We compared our method with some recent works, including MCNN~\cite{Zhang_2016_CVPR_MCNN}, Switching CNN~\cite{8099912_switch}, SANet~\cite{Cao_2018_ECCV_SA}, CSRNet~\cite{Li_2018_CVPR_CSRNet}, ic-CNN~\cite{ranjan2018iterative}, PACNN~\cite{Shi_2019_CVPR_PACNN}, ADCrowdNet~\cite{Liu_2019_CVPR_ADCrowd}, CAN~\cite{Liu_2019_CVPR_CAN}, BL~\cite{Ma_2019_ICCV_BL}, TEDnet~\cite{jiang2019crowd}, HA-CCN~\cite{sindagi2019ha}. The results are shown in Table~\ref*{Tab:Shanghaitech}. Although our model is based on a vanilla CNN architecture, our performance is comparable to the state-of-the-art works. For reference, our network architecture is similar to that of CSRNet~\cite{Li_2018_CVPR_CSRNet}, except that our network does not incorporate dilation convolution layers. As observed, depite our simple architecture, our result is significantly better than that of CSRNet.

\textbf{UCF-QNRF:}
To train and valiate our model in this dataset, we scale the long edge of each crowd image to 1080 pixels, while maintaining the origin aspect ratio. For data augmentation, we crop 4 patches from each image and each patch is 1/4 of the original image size. The learning rate is set as $10^{-6}$ and we decay the learning rate by 10 times at the $630^{th}$, $650^{th}$, $660^{th}$ epoch, respectively. During the distillation stage, the learning rate is set to $10^{-6}$ as well, and the learning rate will decay by 10 times after 80, 100, 110 epochs, respectively.
The comparison results of our method against the state-of-the-art methods, including MCNN~\cite{Zhang_2016_CVPR_MCNN}, MCNN~\cite{Zhang_2016_CVPR_MCNN}, Switching CNN~\cite{8099912_switch}, CL~\cite{Idrees_2018_ECCV_CL}, HA-CCN~\cite{sindagi2019ha}, RANet~\cite{zhang2019relational}, CAN~\cite{Liu_2019_CVPR_CAN}, TEDnet~\cite{jiang2019crowd}, SPN+L2SM~\cite{xu2019learn}, S-DCNet~\cite{xiong2019open}, SFCN~\cite{wang2019learning}, DSSINet~\cite{Liu_2019_ICCV_DSSINet}, MBTTBF-SCFB~\cite{Sindagi_2019_ICCV_MBTTBF}, are shown in Table~\ref{Tab:UCF-QNRF}. As observed, in the term of MAE, our model is demonstrated superior to the latest methods (e.g.~\cite{Sindagi_2019_ICCV_MBTTBF}) for at least more than $4.7\%$ gain.

\textbf{UCSD:}
During training, we upscale the images by two times and we leverage the provided regions of interest.
We first set the Gaussian $\sigma$ as 7 to generate the ground-truth density maps, and the set Gaussian $\sigma$ as 5 and 3 to produce density maps for distillation, respectively. Other settings follow the general training protocal.
As depicted in Table~\ref{Tab:UCSD}, compared with MCNN~\cite{Zhang_2016_CVPR_MCNN}, Switching CNN~\cite{8099912_switch}, ConvLSTM~\cite{xiong2017spatiotemporal}, BSAD~\cite{huang2017body}, ACSCP~\cite{shen2018crowd}, CSRNet~\cite{Li_2018_CVPR_CSRNet}, SANet~\cite{Cao_2018_ECCV_SA}, SPANet~\cite{Cheng_2019_ICCV_SPANet}, ADCrowdNet~\cite{Liu_2019_CVPR_ADCrowd}, our model shows the state-of-the-art performance. Since each image of this dataset contains a small number of people, the difference amongst the comparison methods are not significant. Even so, our approach can still be comparable to the latest methods (e.g.~\cite{Liu_2019_CVPR_ADCrowd}).

\textbf{Qualitative results:} In Fig.~\ref{fig:comparison}, we show several examples in which our results are compared against those of two representative crowd counting approaches, Bayesian loss based method \cite{Ma_2019_ICCV_BL} and CSRNet \cite{Li_2018_CVPR_CSRNet}. Perceptually, since our model generates the density maps with a relatively higher resolution, the visual quality of our produced density maps is much better and similar to the ground-truth.

\begin{figure*}
	\centering
	\setlength{\tabcolsep}{2pt}
	\begin{tabular}{cccccc}
		Input image & $G^{(t=0)}(25)$ & $G^{(t=1)}(20)$ & $G^{(t=2)}(10)$ & $G^{(t=3)}(5)$ & Ground-truth\\
		\includegraphics[width=0.16\textwidth]{images/dismap/IMG_72.jpg}&
		\includegraphics[width=0.16\textwidth]{images/dismap/line_IMG_72_pre_697.jpg}&
		\includegraphics[width=0.16\textwidth]{images/dismap/d2_IMG_72_pre_686.jpg}&
		\includegraphics[width=0.16\textwidth]{images/dismap/d4_IMG_72_pre_678.jpg}&
		\includegraphics[width=0.16\textwidth]{images/dismap/d5_IMG_72_pre_679.jpg}&
		\includegraphics[width=0.16\textwidth]{images/dismap/line_IMG_72_gt_589.jpg}\\
		Estimated count & 697 & 686 & 678 & 679 & 589
	\end{tabular}
	\caption{Illustration of an example that progressively enhances the counting performance via multi-step distillation. }
	\label{fig:dif_map}
\end{figure*}

\subsection{Ablation study}
We perform the ablation study on the dataset, ShanghaiTech Part A~\cite{Zhang_2016_CVPR_MCNN}.
we thoroughly dissect and delve into the structure of our method, including our generated ground-truth density maps, the multi-task loss of our applied network, and our distillation method.

\begin{figure*}
	\centering
	\includegraphics[width=1\linewidth]{images/kernel-comparison.png}
	\caption{Illustration of examples generated by different density map generation methods. }
	\label{fig:dif_kernel}

\end{figure*}

\begin{table}
	\setlength{\tabcolsep}{4pt}
	\centering
	\caption{Ablation study on density map generation, multi-task loss, and distillation.
	}
	\begin{tabular}{c|c||c|c}
		\toprule
		Density map generation & MAE& Distillation & MAE \\
		\midrule
		+ Fixed & 67.91 & + $G^{(t=0)}(25)$ & 62.55 \\
		+ Adaptive & 66.25
		&+ $G^{(t=1)}(20)$       &61.15\\

		+ Non-uniform & 64.58 &  + $G^{(t=2)}(10)$       &\textbf{61.09} \\
		+ Ours & \textbf{62.55}  &+ $G^{(t=3)}(5)$       &61.79\\
		\cline{3-4}
		& & + $G^{(t=0)}(10)$ & 62.72\\
		\midrule
		Single loss & 68.86 & Multi-task loss & 66.25\\

		\bottomrule
	\end{tabular}

	\label{tab:ablation}

\end{table}

\textbf{Density map generation.} We apply different strategies, including fixed kernel size (``+Fixed"), geometric adaptive kernel size (``+Adaptive"), non-uniform kernel estimation (``+Non-uniform"), and our proposed perspective-aware kernel estimation, to generate ``ground-truth" density maps from point annotations to supervise our model. For comparison, the crowd counting network model is trained without distillation.
The results are depicted in the left part of Table~\ref{tab:ablation}. As observed, our generated density maps can benefit the performance of crowd counting network, which significantly outperforms others by at least 2 points in the term of MAE.
In Fig.~\ref{fig:dif_kernel}, we illustrate several examples to show the difference of our proposed density map generator against previous methods. As observed, with the default settings, our approach can produce sharper density map than prior methods, which is able to deliver more accurate spatial information to the counting model and thus leads to better estimated counts.

\textbf{Distillation.} In the right part of Table~\ref{tab:ablation}, we iteratively run our distillation algorithm and measure the performance of the model on each stage. As observed, the distillation progressively promotes the counting performance of the model from 62.55 to 61.09, which is comparable to the state-of-the-art methods. But, when the kernel size shrinks to 5, the performance cannot further improve, due to the limit of distillation.
After each distillation, the model will learn new supervision information on the original basis, optimize the performance of the model, and as the Gaussian kernel decreases, the density map will learn to have sharper position information. For reference, we also apply the density map generated by $G^{(t=0)}(10)$ to train our model from scratch. As observed, its result is worse than our model trained from $G^{(t=0)}(25)$.

We illustrate an example of distillation in Fig.~\ref{fig:dif_map}. Particularly, we observe that, after distilling the model for two time steps (i.e., from $G^{(t=0)}$ to $G^{(t=2)}$), the results of the model are obviously improved, which may be due to the fact that the kernels of the density distribution near the camera becomes smaller and clearer. After distilling our model for more than 2 time steps, the performance hardly gains further improvements.

\textbf{Multi-task loss.} In the bottom row of Table~\ref{tab:ablation}, our baseline is a vanilla CNN-based model without the extra branch for supervising low-resolution density map (i.e., Single loss). Both networks are trained using the density maps produced by the geometric-adaptive kernel. As observed, the performance in the term of MAE (i.e., 68.86) is significantly worse than the one with multi-task loss (i.e., 66.25).

\section{Conclusion and Future Works}
In this paper, we propose a perspective-aware density map generation method that is able to adaptively produce ground-truth density maps from point annotations to train crowd counting model to accomplish superior performance than prior density map generation techniques. Besides, leveraging our density map generation method, we propose an iterative distillation algorithm to progressively enhance our model with identical network structures, without significantly sacrificing the dimension of the output density maps. In experiments, we demonstrate that, with our simple convolutional neural network architecture strengthened by our proposed training algorithm, our model is able to outperform or be comparable with the state-of-the-art methods.

Although our model can obtain satisfactory performance with a simple architecture, we have not validated our approach can be adapted to more complex network models. As the future work, we will explore the possibility of transferring our algorithm to advanced network models. On the other hand, distilling network is time consuming and there are many hyper-parameters for tuning. We will investigate more efficient algorithm to accomplish this purpose.

\ifCLASSOPTIONcaptionsoff
  \newpage
\fi

\bibliographystyle{IEEEtran}