\section{Introduction}
Stochastic frontier analysis (SFA) has a vast literature, both methodological and applied, and empiricists have applied the methods to myriad industries, most notably agriculture, banking, education, healthcare, and energy. A common practice in SFA is to impose parametric assumptions on the error components, but the set of statistical tools to investigate the validity of these assumptions is still limited. This paper expands this set of tools by drawing on recently developed techniques in Extreme Value (EV) theory and by developing new diagnostic tests.

In particular, the parametric stochastic frontier model for cross-sectional data (Aigner et al. 1997) is a leading case of the error component regression model but with the unique feature that one error component ($U$) is a non-negative random variable (e.g., half-normal, exponential), while the other ($W$) is a random variable of unbounded support (e.g., normal, Laplace, Student-t). A common assumption in the stochastic frontier literature is that $W$ is drawn from a normal or Laplace distribution (both thin-tailed distributions). See Aigner et al. (1977) or Horrace and Parmeter (2018), respectively.\footnote{For other parametric specifications of the model see Li (1996), Carree (2002), Tsionas (2007), Kumbhakar et al. (2013), and Almanidis et al. (2014). } However, heavy-tailed distributions are now also being considered. For example, the findings of Wheat et al. (2019) suggest that a cost inefficiency model of highway maintenance costs in England has Student-t errors.\footnote{There are semi-parametric estimators of the model that relax the distributional assumptions on one component and estimate the density of the other using kernel deconvolution techniques. See Kneip et al. (2015), Horrace and Parmeter (2011), Cai et al. (2020), Simar et al. 2017, Hall and Simar (2002), Florens et al. (2020).}  These parametric distributions, such as normal and Student-t, display similar patterns in the middle of their supports but exhibit substantially different tail behaviors. This observation motivates and plays an essential role in our diagnostic tests, which we believe are a timely and appropriate contribution to the literature.

The key idea of our test is as follows. Assuming independence of the error components, the largest order statistics of the composed error term ($Z = W-U$) (approximately) arise from the tails of $W$, because $U$ is one-sided. Also, assuming that $W$ is in the domain of attraction (DOA) of extreme value distributions, the asymptotic distribution of the largest order statistics of $W$ is the EV distribution, which may be fully characterized (after location and scale normalization) by a single parameter that captures its tail heaviness.\footnote{The assumption that $W$ is in the DOA of extreme value distributions is not restrictive, as we shall see.} Then, likelihood ratio statistics for hypotheses on this single parameter can be derived based on the limiting EV distribution.

To be specific, consider the right tail of $W$. If the DOA assumption is satisfied, then tail behavior may be entirely characterized by a \textit{tail index}, $\xi \in \mathbb{R}$. If $\xi = 0$, then $W$ has thin tails. If $\xi > 0$, then $W$ has thick tails. Otherwise, $W$ has bounded support. Under very weak assumptions on the error components, we derive a test that the tails of $W$ are thin ($\xi=0$). We prove that this test is valid whether $Z$ is observed or appended to a regression model (as it is in the stochastic frontier model). If we assume that $U$ is also in the DOA of extreme value distributions and that $W$ is symmetric (a common assumption), we also derive a test that the (right) tail of $U$ is thinner than the left tail of $W$. If we further assume that $W$ is a member of the normal family, then we may test the hypotheses that the tails of $U$ and $W$ are both thin. Therefore, our nonparametric tests are useful diagnostic tools to help empiricists make parametric choices on the distributions of both $U$ and $W$. This is particularly important for the stochastic frontier model for cross-sectional data, where distributional assumptions on the components are typically necessary for the identification of the model's parameters.

The paper is organized as follows. The next section presents the tests. Section 3 provides a simulation study of the power and size of the test. Section 4 applies the tests to a stochastic cost function for a panel of US banks 1998-2005, revealing that the tails of $W$ are not thin. Therefore, a normal or Laplace assumption for $W$ is not justified, and perhaps a Student-t assumption may be appropriate. Section 5 concludes.

\section{Tests of Tail Behavior}
To fix ideas, we begin a review of the DOA assumption and present the test in the case where $Z$ is directly observed in Section \ref{sec: no covariate}. Then in Section \ref{sec: covariate}, we move to the case where $Z$ is appended to a regression model and has to be estimated, which covers the linear regression stochastic frontier model. Additional tests under different sets of weak assumptions are also presented.\footnote{While the analyses that follow are for cross-sectional data, they can easily be applied to panel data, as long as one is willing to assume independence in both the time and cross-sectional dimensions.}
\subsection{The case with no covariates}\label{sec: no covariate}
Consider a random sample of $Z_{i}=W_{i}-U_{i}$ for $i=1,\ldots ,n$, where $U_{i}\geq 0$ represents \textit{inefficiency}, and $W_{i}\in \mathbb{R}$ is \textit{noise} with unbounded support. We start with testing the \textit{shape} of the right tail of $W_{i}$ in a nonparametric way.

The key assumption is that the distribution of $W_{i}$ is within the domain
of attraction of EV distributions. In particular, a cumulative distribution function $F$
is in the domain of attraction of $G_{\xi }$, denoted as $F\in \mathcal{D}%
\left( G_{\xi }\right) $, if there exist constants $a_{n}>0$ and $b_{n}$
such that
\begin{equation*}
\lim_{n\rightarrow \infty }F^{n}\left( a_{n}v+b_{n}\right) =G_{\xi }\left(
v\right)
\end{equation*}%
where $G_{\xi }$ is the generalized EV distribution,
\begin{equation}
G_{\xi }(v)=\left\{
\begin{array}{l}
\exp [-(1+\xi v)^{-1/\xi }]\text{, }1+\xi v\geq 0\text{, for }\xi \neq 0 \\
\exp [-e^{-v}]\text{, }v\in
\mathbb{R}
\text{, }\xi =0%
\end{array}%
\right.  \label{def_G}
\end{equation}%
and $\xi $ is the tail index, measuring the decay rate of the tail.

The domain of attraction condition is satisfied by a large range of commonly
used distributions. If $\xi $ is positive, this condition is equivalent to
regularly varying at infinity, i.e.,
\begin{equation}
\lim_{t\rightarrow \infty }\frac{1-F\left( tv\right) }{1-F\left( t\right) }%
=v^{-1/\xi }\text{ \ \ \ \ for }v>0.  \label{RV}
\end{equation}%
This covers Pareto, Student-t\footnote{The tail index of the Student-t distribution with $\nu$ degrees of freedom is $\xi =1/\nu$.}, and F distributions, for example. The case
with $\xi =0$ covers the normal family, and the case with $\xi <0$
corresponds to distributions with a bounded support.\footnote{%
The uniform distribution has $\xi =-1$, and the triangular distribution has $%
\xi =-1/2$.} See de Haan and Ferreira (2007), Ch.1 for a complete review.

Note that the above notation is for the right tail of $W$, which can be easily adapted to the left tail by considering $-W$. For expositional simplicity, we denote $\xi_{W_{-}}$ and $\xi_{W_{+}}$ as the tail indices for the left and right tails of $W$, respectively. The same notation applies to other variables (e.g., $U$ and $Z$) introduced later.

Returning to SFA, a common assumption is that $W_{i}$ is normal
or Laplace, which implies that $\xi_{W_{+}} =0$. So our hypothesis testing problem
is as follows:%
\begin{equation} \label{eq: hypo1}
H_{0}:\xi_{W_{+}} =0\text{ against }H_{1}:\xi_{W_{+}} >0.
\end{equation} %
If the null hypothesis is rejected, we would then argue that some heavy-tailed
distribution should be used to model the noise and maybe the inefficiency as
well.

To obtain a feasible test, we argue that, since $U_{i}$ is bounded from below at zero, the largest order statistics of $%
Z_{i} $ are approximately stemming from the right tail of $W_{i}$. This is
formalized in Proposition \ref{prop evt}, which requires the following conditions. Let $Z_{n:n}\geq ,\ldots ,\geq
Z_{n:1}$ be the order statistics of $\{Z_{i}\}_{i=1}^{n}$ by descend
sorting. Denote%
\begin{equation*}
\mathbf{Z}_{+}=\left( Z_{n:n},...,Z_{n:n-k+1}\right) ^{\intercal }
\end{equation*}%
as the $k$ largest observations. From now on, we use bold letters to denote vectors.
Denote $F_{W}$ and $Q_{W}(p)=\inf \{y\in
\mathbb{R}
\left. :\right. p\leq F_{W}(y)\}$ as the CDF and the quantile function of $%
W_{i}$, respectively. Write $Q_{W}(1)$ as the right end-point of the support
of $W_{i}$. For a generic column vector $\mathbf{X}$ and scalar $c$, the notation $%
\mathbf{X}-c $ means $\mathbf{X}-(c,\ldots ,c)^{\intercal }$.

\paragraph{Assumption 1}

\begin{description}
\item (i) $\left( U_{i},W_{i}\right) ^{\intercal }$ is \textit{i.i.d.}\

\item (ii) $U_{i}$ and $W_{i}$ are independent.

\item (iii) $U_{i}\geq 0$ with $\mathbb{E}\left[ \left\vert U_{i}\right\vert %
\right] <\infty $ and $W_{i}\in
\mathbb{R}
$ with $Q_{W}(1)=\infty $.

\item (iv) $F_{W}\in \mathcal{D}(G_{\xi_{W_{+}} })$ with $\xi_{W_{+}} \geq 0$. In addition, $%
F_{W}(\cdot )$ is twice continuously differentiable with bounded
derivatives, and the density $f_{W}(\cdot )$ satisfies that $\partial
f_{W}(t)/\partial t\nearrow 0$ as $t\rightarrow \infty $ on $[c,\infty )$
for some constant $c$.
\end{description}

Assumptions 1(i)-(iii) are common in the SFA literature (see Horrace and Parmeter (2018) and the references therein). Assumption 1(iv) requires the tail of $F_{W}$ to be within the domain of
attraction of EV distributions with an infinite upper bound. Moreover, it requires
that the density derivative monotonically increases to zero. This is a mild
assumption and is satisfied by many commonly used distributions. For example,
the normal distribution is covered as seen by%
\begin{equation*}
\frac{\partial f_{W}(t)}{\partial t}\propto -t\exp (-\frac{t^{2}}{2}%
)\nearrow 0\text{ as }t\rightarrow \infty ,
\end{equation*}%
and the Pareto distribution is covered as seen by%
\begin{equation*}
\frac{\partial f_{W}(t)}{\partial t}\propto (-\alpha -1)x^{-\alpha
-2}\nearrow 0\text{ as }t\rightarrow \infty \text{ for some }\alpha >0.
\end{equation*}

Under Assumption 1, the following proposition derives the asymptotic distribution of $\mathbf{Z}_{+}$.

\begin{proposition}
\label{prop evt}
Suppose Assumption 1 holds. Then, there exist sequences of
constants $a_{n}$ and $b_{n}$ such that for any fixed $k$%
\begin{equation*}
\frac{\mathbf{Z}_{+}-b_{n}}{a_{n}}\overset{d}{\rightarrow }\mathbf{V}%
_{+}=\left( V_{1},...,V_{k}\right) ^{\intercal }
\end{equation*}%
where the joint density of $\mathbf{V}_{+}$ is given by $f_{\mathbf{V}_{+}|\xi_{W_{+}}
}(v_{1},...,v_{k})=G_{\xi_{W_{+}} }(v_{k})\prod_{i=1}^{k}g_{\xi_{W_{+}} }(v_{i})/G_{\xi_{W_{+}}
}(v_{i})$ on $v_{k}\leq v_{k-1}\leq \ldots \leq v_{1}$, and $g_{\xi_{W_{+}}
}(v)=dG_{\xi_{W_{+}} }(v)/dv$.
\end{proposition}

The proof is in Appendix A. This proposition implies that the distributions of $Z_{i}$ and $W_{i}$ share
the same (right) tail shape, which is entirely characterized by the tail index $\xi_{W_{+}} $. Such tail equivalence does not hold, however, for the left tails due to the existence of $U$. This is studied in Section \ref{sec sym} under the additional assumption that $W$ is symmetric.

If the constants $a_{n}$ and $b_{n}$ were known, $\mathbf{Z_{+}}$ is then
approximately distributed as $\mathbf{V_{+}}$, and the limiting problem is
reduced to the well-defined finite sample problem: constructing some inference
method based on one draw $\mathbf{V_{+}}$ whose density $f_{\mathbf{V_{+}}|\xi_{W_{+}} }$ is
known up to $\xi_{W_{+}} $. However, $a_{n}$ and $b_{n}$ depend on $F_{W}$ and hence are unknown \textit{a priori}.

To avoid the need for knowledge of $a_{n}$ and $b_{n}$, we consider the following
self-normalized statistic
\begin{eqnarray}
\mathbf{Z}_{+}^{\ast } &=&\frac{\mathbf{Z}_{+}-Z_{n:n-k+1}}{%
Z_{n:n}-Z_{n:n-k+1}}  \label{max_inv} \\
&=&\left( 1,\frac{Z_{n:n-1}-Z_{n:n-k+1}}{Z_{n:n}-Z_{n:n-k+1}},...,\frac{%
Z_{n:n-k+2}-Z_{n:n-k+1}}{Z_{n:n}-Z_{n:n-k+1}},0\right) ^{\intercal }.  \notag
\end{eqnarray}%
It is easy to establish that $\mathbf{Z}^{\ast }$ is maximally invariant with
respect to the group of location and scale transformations (cf., Lehmann and Romano (2005), Ch.6). In words, the
estimator constructed as a function of $\mathbf{Z}^{\ast }$ remains
unchanged if data are shifted and multiplied by any non-zero constant. This
makes senses since the tail shape should be preserved no matter how data are
linearly transformed. This invariance property allows us to construct nonparametric tests for a stochastic frontier model that is otherwise not identified without parametric assumptions on $U$ and $W$.\footnote{In particular the non-zero expectation of $U$ precludes identification of unknown parameter $\delta$ in the model $Z_i=\delta+W_i-U_i$. } As such, our tests do not reveal anything about the location or the scale of the error components.

The continuous mapping theorem and Proposition \ref{prop evt} imply that for
any fixed $k$, as $n\rightarrow \infty $,%
\begin{equation*}
\mathbf{Z}_{+}^{\ast }\overset{d}{\rightarrow }\mathbf{V}_{+}^{\ast }\equiv
\frac{\mathbf{V}_{+}-V_{k}}{V_{1}-V_{k}}.
\end{equation*}%
The CDF\ of $\mathbf{V}_{+}^{\ast }$ can be calculated via change of
variables as%
\begin{equation}
f_{\mathbf{V}_{+}^{\ast }\mathbf{|}\xi_{W_{+}} }\left( \mathbf{v}_{+}^{\ast }\right) =\Gamma
\left( k\right) \int_{0}^{b_{0}\left( \xi_{W_{+}} \right) }t^{k-2}\exp \left(
-(1+1/\xi_{W_{+}} )\sum_{i=1}^{k}\log (1+\xi_{W_{+}} v_{i}^{\ast }t)\right) dt,
\label{density_xs}
\end{equation}%
where $\mathbf{v}_{+}^{\ast }=(v_{1}^{\ast },\ldots ,v_{k}^{\ast })$, $%
b_{0}\left( \xi \right) =\infty $ if $\xi \geq 0$ and $-1/\xi $ otherwise,
and $\Gamma \left( k\right) $ is the gamma function. Note that the
invariance restriction costs two degrees of freedom since the first and last
elements of $\mathbf{V}_{+}^{\ast }$ are always 1 and 0, respectively. We
calculate this density by numerical quadrature.

Given $f_{\mathbf{V}_{+}^{\ast }\mathbf{|}\xi_{W_{+}} }$, we can construct the
generalized likelihood-ratio test for problem (\ref{eq: hypo1}). Since the alternative hypothesis is composite, we follow Andrews and Ploberger (1994) and Elliott et al. (2015) to consider the weighted average alternative
\begin{equation*}
\int f_{\mathbf{V}_{+}^{\ast}|\xi_{W_{+}}}(\cdot)w(\xi_{W_{+}} )d\xi_{W_{+}},
\end{equation*}
where $w(\cdot )$ is a weighting function that reflects the importance of
rejecting different alternative values. Then our test is constructed as \footnote{In later sections, we set $w(\cdot )$ to be the standard uniform distribution over $(0,1)$ for simplicity.}
\begin{equation}
\varphi (\mathbf{v}_{+}^{\ast })=\mathbf{1}\left[ \frac{\int f_{\mathbf{V}_{+}^{\ast
}|\xi_{W_{+}} }(\mathbf{v}_{+}^{\ast })w(\xi_{W_{+}} )d\xi_{W_{+}} }{f_{\mathbf{V}_{+}^{\ast }|0}(\mathbf{v}_{+}%
^{\ast })}>\mathrm{cv}(k,\alpha )\right] ,  \label{test}
\end{equation}%
 where the critical value $\mathrm{cv}(k,\alpha )$ depends on $k$ and the level of significance $\alpha $. We can
obtain it by simulation. By Proposition \ref{prop evt} and the continuous mapping theorem, this test controls size asymptotically as $\lim_{n\rightarrow \infty }\varphi (\mathbf{Z}_{+}^{\ast })=\alpha$.

We end this subsection by briefly discussing the choice of $k$, that is, the number of the largest order statistics used to approximate the EV distribution. On the one hand, larger $k$ means including more mid-sample observations, which induces a larger finite sample bias in the EV approximation. On the other hand, smaller $k$ provides a better asymptotic approximation but uses less sample information, leading to a lower power test. This trade-off leads to difficulty in theoretical justification of an optimal $k$ in standard EV theory literature (cf., M\"{u}ller and Wang (2017)). It is even more difficult, if at all possible, in our case, since we only observe $Z$, and not $W$. Nonetheless, our asymptotic arguments show that the test (\ref{test}) controls size for any fixed $k$, as long as $n$ is sufficiently large. Figure \ref{fig:power_curves} depicts the asymptotic power of the test (\ref{test}) with $\mathbf{V}_{+}^{\ast}$ generated from the density (\ref{density_xs}) based on 10,000 simulation draws. The test controls size for all values of $k$ by construction and has reasonably large power when $k$ exceeds 20.

With ideas fixed, we now turn to the regression version of the test, with application to SFA.

\subsection{The case with covariates: SFA}\label{sec: covariate}

Now consider the linear regression with
\begin{equation*}
Y_{i}=\mathbf{X}_{i}^{\intercal }\beta _{0}+Z_{i},
\end{equation*}%
where $Z_{i}=-U_{i}+W_{i}$ is as in the previous section, and $\beta _{0}$ is
some pseudo-true parameter in some compact parameter space. This could be a Cobb-Douglas production function (in logarithms), where $Y$ is productive output and $U$ is now called \textit{technical efficiency}, which measures distance ($U_i$) from a stochastic frontier ($\mathbf{X}_{i}^{\intercal }\beta _{0}+W_{i}$). The slopes ($\beta_0$) are marginal products of the productive inputs, $\mathbf{X}_{i}$. It could also be a stochastic cost function if we multiply $U$ by $-1$. Suppose we have some estimator, $\hat\beta$ of $\beta_0$. The following assumption is imposed to construct our diagnostic test.

\paragraph{Assumption 2}

\begin{description}
\item (i) $\left( X_{i},U_{i},W_{i}\right) ^{\intercal }$ is i.i.d.\

\item (ii) $U_{i}$ and $W_{i}$ are independent.

\item (iii) $U_{i}\geq 0$ with $\mathbb{E}\left[ \left\vert U_{i}\right\vert %
\right] <\infty $ and $W_{i}\in \mathbb{R}$ with $Q_{W}(1)=\infty $.

\item (iv) $F_{W}\in \mathcal{D}(G_{\xi_{W_{+}} })$ with $\xi_{W_{+}} \geq 0$. In addition, $%
F_{W}(\cdot )$ is twice continuously differentiable with bounded
derivatives, and the density $f_{W}(\cdot )$ satisfies that $\partial
f_{W}(t)/\partial t\nearrow 0$ as $t\rightarrow \infty $ on $[c,\infty )$
for some constant $c$.

\item (v) $\left\vert \left\vert \hat{\beta}-\beta _{0}\right\vert
\right\vert \sup_{i}\left\vert \left\vert \mathbf{X}_{i}\right\vert
\right\vert =o_{p}(n^{\xi_{W_{+}} })$, if  $\xi_{W_{+}} >0$. $\left\vert
\left\vert \hat{\beta}-\beta _{0}\right\vert \right\vert \sup_{i}\left\vert
\left\vert \mathbf{X}_{i}\right\vert \right\vert /f_{W}\left(
Q_{W}(1-1/n)\right) \left. =\right. o_{p}(1)$, otherwise.
\end{description}

Assumption 2 is similar to Assumption 1 with additional restrictions on the covariate $\mathbf{X}$. In particular, Assumption 2(v) bounds the norm of $\hat{\beta}$ and $\left\vert \left\vert
\mathbf{X}_{i}\right\vert \right\vert $. A sufficient condition when $\xi_{W_{+}} $
is positive is that $\left\vert \left\vert \hat{\beta}-\beta _{0}\right\vert
\right\vert =O_{p}(n^{-1/2})$ and $\sup_{i}\left\vert \left\vert \mathbf{X}%
_{i}\right\vert \right\vert =o_{p}(n^{1/2})$, which is easily satisfied in
many applications.\footnote{Even though $\mathbb{E}\left[ \left\vert U_{i}\right\vert \right] \ne 0$, ordinary least squares (OLS) will typically suffice for $\hat\beta$, because our test is invariant to relocation. } When $\xi $ is zero, we need slightly stronger bounds.
Straightforward calculations show that the normal distribution satisfies
Assumption 2(v) for the $\xi_{W_{+}} = 0$ case, if $\left\vert \left\vert \hat{\beta}-\beta
_{0}\right\vert \right\vert =O_{p}\left( n^{-1/2}\right) $ and $%
\sup_{i}\left\vert \left\vert \mathbf{X}_{i}\right\vert \right\vert \left.
=\right. O_{p}(n^{1/2-\varepsilon })$ for some $\varepsilon >0$. This is
seen by $1/f_{W}\left( Q_{W}\left( 1\left. -\right. 1/n\right) \right)
\left. \leq \right. O(\log (n))$ (cf.\ Example 1.1.7 in de Haan and Ferreira (2007)).

Denote $\hat{Z}_i$ as the OLS residuals and
\begin{equation*}
\mathbf{\hat{Z}}_{+}=\left( \hat{Z}_{n:n},...,\hat{Z}_{n:n-k+1}\right) ^{\intercal }
\end{equation*} the largest $k$ order statistics. Then given Assumption 2, the following proposition derives the asymptotic distribution of $\mathbf{\hat{Z}}_{+}$
\begin{proposition}
\label{prop evtx}
Suppose Assumption 2 holds. Then, there exist sequences of
constants $a_{n}$ and $b_{n}$ such that for any fixed $k$%
\begin{equation*}
\frac{\mathbf{\hat{Z}}_{+}-b_{n}}{a_{n}}\overset{d}{\rightarrow }\mathbf{V}%
_{+}
\end{equation*}%
where the joint density of $\mathbf{V}_{+}$ is the same as in Proposition %
\ref{prop evt}.
\end{proposition}

The proof is in Appendix A. Proposition \ref{prop evtx} implies that the largest order statistics of the
regression residuals satisfy the same convergence as the no-covariate case. In other words, the estimation error from the OLS becomes negligible so that the largest order statistics are stemming from the right tail of $W$ asymptotically.
This validates the construction of the test (\ref{test}) by replacing
$\mathbf{Z}_{+}^{\ast }$ with $\mathbf{\hat{Z}}_{+}^{\ast }$, where%
\begin{equation*}
\mathbf{\hat{Z}}_{+}^{\ast }=\frac{\mathbf{\hat{Z}}_{+}-\hat{Z}_{n:n-k+1}}{%
\hat{Z}_{n:n}-\hat{Z}_{n:n-k+1}}.
\end{equation*}
Proposition \ref{prop evtx} and the continuous mapping theorem, we similarly have  $\lim_{n\rightarrow \infty }\varphi (\mathbf{\hat{Z}_{+}}^{\ast })=\alpha$.

\subsection{Symmetry of noise $W$}\label{sec sym}

The previous analysis studies the right tail of $W$ (and equivalently $Z$). Suppose we assume $W$ has a symmetric distribution, then the tail indices of
both tails of $W$ become equivalent, and hence we can learn about the tail of $U$ using the left tail index
of $Z$. To this end, we make the following
additional assumption.

\paragraph*{Assumption 3}

\begin{description}
\item (i) $W_{i}$ is symmetric at zero.

\item (ii) $F_{U}\in \mathcal{D}( G_{\xi _{U_{+}}})$ with $\xi_{U_{+}} \geq0$.
\end{description}

Assumption 3(i) implies that $\xi_{W_{-}}=\xi_{W_{+}}$, and the condition that $U>0$ implies its left tail index is negative. Therefore, in this subsection only, we simply denote $\xi_{U}$ and $\xi_{W}$ as the right tail indices of $U$ and $W$, respectively. Now we can test if $U$ has a thinner or equal right tail than $W$ by specifying
the following hypothesis testing problem,%
\begin{equation}
H_{0}:\xi _{U}\leq \xi _{W}\text{ against }H_{1}:\xi _{U}>\xi _{W}\text{.}
\label{hypo uw}
\end{equation}%
Moreover, if $W$ is in the normal or Laplace family ($\xi _{W}=0$), since we limit the tail indices to be non-negative, the null hypothesis
then reduces to $\xi _{U}=\xi _{W}=0$.

Under the null hypothesis of (\ref{hypo uw}), $W$ is the leading term in $Z$
in both the left and right tails. Then the DOA
assumption for both $W$ and $U$ implies that $\xi _{Z_{-}}=\max \{\xi
_{U},\xi _{W}\}$, and Proposition \ref{prop evtx} entails $\xi
_{Z_{+}}=\xi _{W}$. Therefore, the above testing problem becomes equivalent to
\begin{equation}
H_{0}:\xi _{Z_{-}}=\xi _{Z_{+}}\text{ against }H_{1}:\xi _{Z_{-}}>\xi
_{Z_{+}}\text{.}  \label{hypo equaltail}
\end{equation}

We now construct a test for (\ref{hypo equaltail}). Define $\mathbf{\hat{Z}}%
_{-}$ as the smallest $k$ order statistics of the estimation residuals, that is,
\begin{equation*}
\mathbf{\hat{Z}}_{-}=\left( \hat{Z}_{n:1},\hat{Z}_{n,2},\ldots ,\hat{Z}%
_{n,k}\right) ^{\intercal }
\end{equation*}%
and its self-normalized analogue as%
\begin{equation*}
\mathbf{\hat{Z}}_{-}^{\ast }=\frac{\mathbf{\hat{Z}}_{-}-\hat{Z}_{n:k}}{\hat{Z%
}_{n:1}-\hat{Z}_{n:k}}.
\end{equation*}%
The following proposition establishes that $\mathbf{\hat{Z}}_{-}^{\ast }$
asymptotically has the EV distribution with tail index $\xi _{Z_{-}}$ and is
independent from $\mathbf{\hat{Z}}_{+}^{\ast }$.

\begin{proposition}
\label{prop sym}Suppose Assumptions 2 and 3 hold. Then, for any fixed $k$,%
\begin{equation*}
\binom{\mathbf{\hat{Z}}_{-}^{\ast }}{\mathbf{\hat{Z}}_{+}^{\ast }}\overset{d}%
{\rightarrow }\binom{\mathbf{V}_{-}^{\ast }}{\mathbf{V}_{+}^{\ast }}\text{
as }n\rightarrow \infty ,
\end{equation*}%
where $\mathbf{V}_{-}^{\ast }$ and $\mathbf{V}_{+}^{\ast }$ are independent
and both EV distributed with density (\ref{density_xs}) and tail indices $\xi _{Z_{-}}$ and $\xi _{Z_{+}}$, respectively.
\end{proposition}

The proof is in Appendix A. Given the above proposition, we aim to construct a generalized likelihood ratio
test for (\ref{hypo equaltail}) as follows,
\begin{equation}
\varphi_{\pm} \left( \mathbf{v}_{-}^{\ast },\mathbf{v}_{+}^{\ast }\right) =%
\mathbf{1}\left[ \frac{\int_{\{\left( \xi _{-},\xi _{+}\right) \in \Xi
^{2}:\xi _{+}<\xi _{-}\}}f_{\mathbf{V}_{-}^{\ast }|\xi _{-}}(\mathbf{v}%
_{-}^{\ast })f_{\mathbf{V}_{+}^{\ast }|\xi _{+}}\left( \mathbf{v}_{+}^{\ast
}\right) w\left( \xi _{-},\xi _{+}\right) d\xi_{-}d\xi_{+}} {\int_{\Xi }f_{\mathbf{V}%
_{-}^{\ast }|\xi }(\mathbf{v}_{-}^{\ast })f_{\mathbf{V}_{+}^{\ast }|\xi
}\left( \mathbf{v}_{+}^{\ast }\right) d\Lambda \left( \xi \right) }>\mathrm{%
cv}(k,\alpha )\right] ,  \label{test eq}
\end{equation}%
where $\Xi$ denotes the parameter space of the tail indices, and $w\left( \cdot ,\cdot \right) $ is the weighting function for the
alternative hypothesis as in (\ref{test}). We set $\Xi$ to be $[0,1)$ to cover all distributions with a finite mean and $w(\cdot)$ to be uniform over the alternative space. The weight $\Lambda \left( \cdot \right) $ can be considered as the least favorable distribution, which we discuss more now.

Note that the null hypothesis of (\ref{hypo equaltail}) is composite. We need to
control size uniformly over all $\xi _{Z_{-}}=\xi _{Z_{+}}\in \Xi $. To that end, we can transform the composite null into a simple one by considering the weighted average density with respect to the
weight $\Lambda $. Together with a suitably chosen the critical value, this test (\ref{test eq}) maintains the uniform size control. Now the problem reduces to determining an appropriate weight $\Lambda$. Elliott et al. (2015) study the generic
hypothesis testing problem where a nuisance parameter exists in the null
hypothesis. We tailor their argument for our test (\ref{test eq}) and adopt
their computational algorithm for implementation. In particular, $\Lambda
\left( \cdot \right) $ and $\mathrm{cv}(k,\alpha )$ are numerically
calculated only once by the authors instead of the empiricists who use our test. They only need to construct the order statistics $\mathbf{\hat{Z}}_{-}^{\ast }$ and $\mathbf{\hat{Z}}_{+}^{\ast }$ and numerically evaluate the density. We
provide more computational details in the Appendix and the corresponding MATLAB code in the supplemental materials. By the
continuous mapping theorem and Proposition \ref{prop sym}, for any fixed $k$%
, $\lim \sup_{n\rightarrow \infty }\mathbb{E}\left[ \varphi _{\pm }\left(
\mathbf{\hat{Z}}_{-}^{\ast },\mathbf{\hat{Z}}_{+}^{\ast }\right) \right]
\leq \alpha $ under the null hypothesis of (\ref{hypo equaltail}).

As we discussed above, the hypothesis testing problem (\ref{hypo equaltail})
simplifies to
\begin{equation*}
H_{0}:\xi _{Z_{-}}=\xi _{Z_{+}}=0\text{ against }H_{1}:\xi _{Z_{-}}>\xi
_{Z_{+}}=0\text{,}
\end{equation*}%
if $W$ is assumed to be in the normal family ($\xi_{W}=0$). Proposition \ref{prop sym}
implies $\mathbf{\hat{Z}}_{-}^{\ast }$ and $\mathbf{\hat{Z}}_{+}^{\ast }$
are asymptotically independent and both of them are EV distributed. Then accordingly, our test (\ref{test eq}%
) reduces to%
\begin{equation*}
\mathbf{1}\left[ \frac{\int_{\Xi }f_{\mathbf{V}_{-}^{\ast }|\xi _{-}}(%
\mathbf{v}_{-}^{\ast })dW\left( \xi _{-}\right) }{f_{\mathbf{V}_{-}^{\ast
}|0}(\mathbf{v}_{-}^{\ast })}>\mathrm{cv}(k,\alpha )\right] ,
\end{equation*}%
which is identical to (\ref{test}). This suggests that we can simply substitute $%
\mathbf{\hat{Z}}_{-}^{\ast }$ into (\ref{test}) for implementation.

\section{Simulation Study}

\subsection{Hypothesis testing about noise $W$}

We set $w(\cdot )$ to be the uniform weight on $\left[ 0,0.99\right] $ to
include all distributions with a finite mean and the level of significance
to be $0.05$. In Table 1, we report the small sample rejection probabilities
of the test (\ref{test}). We generate $U_{i}$ from the right half-standard
normal and the right half-Laplace(0,1) distributions and $W_{i}$ from four
distributions: standard normal, Laplace(0,1) (denoted La(0,1)) Student-t(2),
Pareto(0.5) and F(4,4). The normal and Laplace distributions correspond to
the null hypothesis, and the other three are alternative hypotheses. The
results suggest that the test (\ref{test}) has an excellent performance in size
and power. Note that when $k=50$ and $n=100$, we essentially include too
many mid-sample observations so that the EV approximation is
poor.

Now we consider the linear regression model that $Y_{i}=\mathbf{X}%
_{i}^{\intercal }\beta _{0}+Z_{i}$ with $\mathbf{X}_{i}=\left(
1,X_{2i}\right) ^{\intercal }$ and $\beta _{0}=\left( 1,1\right) ^{\intercal
}$. We assume $X_{2i}\sim \mathcal{N}\left( 0,1\right) $ and independent
from $Z_{i}$. Table 2 reports the rejection probabilities of our test (\ref%
{test}). Findings are similar to those in Table 1.

\subsection{Hypothesis testing about inefficiency $U$ and noise $W$}

Consider the hypothesis testing problem (\ref{hypo equaltail}). We implement
the test (\ref{test eq}) with the same setup as above. Table 3 reports the
rejection probabilities under the null and alternative hypotheses. We
make the following observations. First, the test controls size well
unless $k$ is too large relative to $n$, as seen in the column with $n=100$
and $k=50$. This is again because we are using too many mid-sample observations to
approximate the tail so that the EV convergence in Propositions 1-3 provides
poor approximations. Second, the test has good power properties as seen from
the last five rows. In particular, using only the largest 50 order
statistics from 1000 observations leads to the power of 0.94. Finally, the
power decreases as the alternative hypothesis becomes closer to the null, as
we move down along rows.

Now we consider the special case where $W$ is in the normal family. Then
we implement (\ref{test}) with $\mathbf{\hat{Z}}_{-}^{\ast }$ as the input.
Table 4 contains the rejection probabilities under the null and alternative
hypotheses. The rows with $F_{U}$ being half-normal or Laplace correspond to
the size under the null hypothesis, while other rows the power under the
alternative hypothesis. The new test has excellent size and power properties.

\section{Empirical illustration}

We illustrate the new method using the US bank data collected by Feng and Serletis (2009). The data are a sample of US banks covering the period from 1998 to 2005 (inclusive). After deleting banks with negative or zero input prices, we are left with a balanced panel of 6,010 banks observed annually over the 8-year period. A more detailed description of
the data may be found in Feng and Serletis (2009). Here we specify a stochastic cost function, letting $Z=W+U$, so $U \geq 0$ is cost inefficiency, and more inefficient banks have higher total costs, $Y$. Since our tests are designed for cross-sectional data, we divide the original panel data into
cross-sections (one for each year) and regress the logarithm of total bank cost on
a constant and the logarithms of six control variables, including the wage
rate for labor, the interest rate for borrowed funds, the price of physical
capital, and the amounts of consumer loans, non-consumer loans, and
securities. Since the object of interest is the cost function, we multiply the OLS residuals by $-1$ and take the smallest and the largest $k\in\{25, 50, 75, 100\}$ order
statistics, respectively, to implement the test (\ref{test}). The p-values
are reported in Table 5. Under the assumption that $W$ is symmetric\footnote{The symmetry assumption is reasonable here and is imposed in Feng and Serletis (2009).}, these small p-values suggest that $W$ has heavy tails on both sides, so a Student-t assumption (e.g., Wheat, Stead, and Greene, 2019) is more appropriate.

\section{Concluding remarks}
We derive several nonparametric tests of the tail behavior of the error components in the stochastic frontier model. The tests are easy to implement in MATLAB and are useful diagnostic tools for empiricists.

Often a first-step diagnostic tool for SFA is to calculate the skewness of the OLS residuals to see if they are properly skewed. See Waldman (1982), Simar and Wilson (2010), and Horrace and Wright (2020). If they are positively skewed, the maximum likelihood estimator of the variance of inefficiency is zero, and OLS is the maximum likelihood estimator of $\beta_0$. If they are negatively skewed, then OLS is not a stationary point in the parameter space of the likelihood, and the stochastic frontier model is well-posed. After calculating negatively skewed OLS residuals, a useful second-step diagnostic tool is to implement our nonparametric tests to understand the tail behaviors of the error component distributions and to guide parametric choices subsequently .


\section{Introduction}
Convolutional neural networks (CNNs) are the method of choice for image classification \cite{RanzatoHBL07, Alexnet, VGG, He_2016_CVPR, imagenet}. CNNs capture multi-scale contextual information in an image via subsampling the intermediate features through the network layers \cite{RanzatoHBL07,SchererMB10}. This approach is successful due to the expansion of the receptive fields, which are large enough to capture the context of the image for classification.

In this work, we consider dense prediction tasks, which are popular in computer vision \cite{taskonomy2018}. One desideratum of dense prediction tasks is to have pixel-accurate predictions, for example in semantic segmentation \cite{He_2016_CVPR,long-shelhamer-fcn-2015,Noh_2015_ICCV,U-net,Hariharan_2015_CVPR} or depth estimation \cite{eigen-nips14, eigen-iccv15}. Even small inaccuracies, such as missing a small object lying on the street, may lead to an accident of an autonomous vehicle.

\begin{figure*}[!htb]
    \centering
    \includegraphics[width=0.95\linewidth]{gfx/split-merge.pdf}
    \caption{\textbf{Split and Merge Pooling.} The illustration of split and  merge pooling layers with a window size of $2\times 2$. The advantage of splitting the input into batches is to make it possible to process each part of the input (\ie each split batch) independently. In this example, after splitting $A$ into batches in $B$, we can send each batch $B^{(i,j)}$ to a different GPU and continue the forward pass from this point onward on multiple GPUs, or process one batch at a time on a single GPU. This enables us to execute dense prediction tasks with very deep networks and for large images, while always preserving the spatial information.
    }
    \label{fig:split-merge}
\end{figure*}

Most state-of-the-art approaches for dense prediction adapt existing CNNs which were originally designed for image classification.
However, these adapted CNNs lose a vast amount of spatial information due to subsampling. As a result this reduces the prediction performance, mostly because of missing small objects or predicting coarse and inaccurate object boundaries. In order to mitigate this problem, other methods are proposed, such as progressive upsampling \cite{eigen-nips14, eigen-iccv15, long-shelhamer-fcn-2015}, deconvolution (or transpose convolutions) \cite{Noh_2015_ICCV}, skip connections \cite{U-net,Hariharan_2015_CVPR,badrinarayanan2017segnet}, utilizing multiple scales of features \cite{jafari_icra17, zhao2017pspnet}, and attention mechanism \cite{zhao2018psanet}.
Another line of work preserves the resolution by replacing subsampling layers with dilated convolution layers \cite{ChenPKMY14,YuKoltun2016,Yu2017} and widely used in other approaches \cite{zhao2017pspnet, zhao2018psanet,SunXLW19,WangSCJDZLMTWLX19}. In theory, it is possible to design a dilated convolution network without using any subsampling operation to obtain an output with the same size as the input. Unfortunately, for large input sizes the training and inference of the resulting CNN is extremely slow or even sometimes intractable due to limited amount of memory on a GPU, hence subsampling is still needed in practice.

In this work, we propose novel pooling layers called Split-Merge Pooling (SMP). The split pooling layer splits (re-arranges) a feature tensor along its spatial dimension and treats the resulting splits as individual batches (see Fig.~\ref{fig:split-merge}). The split pooling reduces the spatial size with a fixed scaling factor, related to the size of the non-overlapping pooling window. The merge layer acts as the reverse of split pooling operation, i.e., it receives the split batches and re-arranges the elements of the batches to their original locations. In our experiments, we replace all subsampling layers of standard ResNet networks by our split pooling. Finally, we merge the resulting split batches using merge pooling layers to output full resolution prediction. In this way, we do not lose any spatial information since it is preserved in the split batches.
The batches are processed independently after the split, which has the major advantage that after the split each batch can be processed on a different GPU. This enables us to perform dense prediction tasks with very large networks on large images, while always preserving the spatial information.
Furthermore, to increase training speed, we introduced the Shrink and Expand pooling layers as a batch-subsampling of the split and merge pooling.

\noindent
In summary, our contributions are as follows:
\begin{itemize}
    \item We propose a novel pooling method called Split-Merge Pooling (SMP) which enables the unique mapping of each input element to one output element, without losing any spatial information.
    \item We propose a batch-subsampling variant of SMP, termed Shrink-Expand Pooling, to make the training efficient and tractable for very deep networks.
    \item To show the effectiveness of SMP on dense prediction tasks, we choose semantic segmentation and apply SMP to ResNet networks with varying depths. We chose ResNet as our baseline since it is used as the backbone in state-of-the-art approaches. For semantic segmentation of large images from Cityscapes \cite{Cordts2016Cityscapes} and GTA-5 \cite{Richter_2016_ECCV} datasets, our SMP networks outperform the corresponding ResNet networks by a significant margin, with up to $6.8\%$ in IoU score.
    \item We even observe that a SMP version of a shallow ResNet (ResNet18) outperforms the original ResNet101 by $2.8\%$ in IoU score, although ResNet101 is 3.8 times deeper than ResNet18.
\end{itemize}

\section{Related Work}
Utilizing CNNs for image classification became very popular with the introduction of the Imagenet challenge \cite{imagenet} and after some popular network architectures such as Alexnet \cite{Alexnet}, VGG \cite{VGG} and ResNet \cite{He_2016_CVPR} emerged.
Pooling layers in CNNs were originally introduced for image classification tasks on MNIST dataset \cite{RanzatoHBL07}. The aim of pooling layers is to summarize the information over a spatial neighborhood.

\paragraph{Dense prediction tasks.} Most of the dense prediction models are based on adapted versions of image classification networks. Eigen \etal \cite{eigen-nips14, eigen-iccv15} adapt the Alexnet \cite{Alexnet} for single image depth estimation. Since the output of such networks is very coarse, they upsample the output progressively and combine it with local features from the input image. Long \etal \cite{long-shelhamer-fcn-2015} introduce the first fully covolutional network (FCN) for semantic image segmentation by adapting VGG network \cite{VGG} for this task. They map intermediate features with higher resolution to label space and combine them with the coarse prediction of the network in order to recover the missing spatial information. Noh \etal \cite{Noh_2015_ICCV} introduce more parameters in a decoder consisting of transpose convolutions (Deconvolution layers) to upsample the coarse output of the encoder. However, it is still difficult to recover the missing information from the coarse output with this approach. Furthermore, other approaches \cite{U-net, Hariharan_2015_CVPR,badrinarayanan2017segnet} use skip connections between encoder and decoder to obtain more detailed information from lower level features.

Yu \etal introduced the concept of dilated convolutions \cite{YuKoltun2016} and later developed the dilated residual network (DRN) \cite{Yu2017}, which uses dilated convolutions to preserve the spatial information. However, they still use three subsampeling layers to reduce the spatial resolution to make training  feasible, i.e. to fit the model into memory. Therefore, DRNs lose spatial information and need upsampling to obtain the full size output prediction. In contrast, our approach does not lose any spatial information. The batch-based design of our split pooling gives us the flexibility of distributed processing of batches on multiple GPUs or sequential processing of batches on one GPU during inference time (see Sec \ref{sec:split-merge}). Furthermore, for training, we learn our network weights using only one subsampled split batch to make the training faster and tractable (see Sec \ref{sec:shrink-expand}).

\begin{figure*}[!htb]
    \centering
    \includegraphics[width=\linewidth]{gfx/shrink-expand-seq.pdf}
    \caption{\textbf{Shrink and Expand Pooling example.} The illustration of two pairs of shrink and expand pooling layers with window size of $2\times 2$ and sampling location of $(1,0)$ and $(1,1)$. During training the expand pooling layers back-propagate the error of sparse elements.}
    \label{fig:shrink-expand-seq}
\end{figure*}

\section{Method}
The goal of this work is to design a set of pooling layers that preserve  spatial information, in contrast to subsampling operations such as max-pooling and striding convolutions. Additionally, the new pooling layers have to be applicable to very deep networks. Losing spatial information caused by subsampling is a prevalent problem for dense prediction tasks, such as object detection, semantic segmentation, instance segmentation, or depth estimation. Small objects in the input image can get lost completely with subsampling. Moreover, recovering precise object boundaries in an image becomes challenging due to missing spatial information. These problems, caused by losing spatial information, obviously reduced prediction accuracy. However, simply storing all the information with the original spatial resolution is not a practical option since it causes storage issues, in particular for training on large image datasets.

Our idea is to preserve all the spatial information by splitting the feature tensor into multiple downsampled batches (see Fig.~\ref{fig:split-merge}) instead of preserving it inside the original spatial resolution. In this way, we have both advantages of large receptive fields due to standard downsampling, as well as keeping all the available spatial information for obtaining precise dense predictions.

\subsection{Split-Merge Pooling}
\label{sec:split-merge}
The split pooling layer splits an input feature tensor along spatial dimensions and outputs each split as a batch. The number of output batches depends on the window size of split pooling, \eg a split pooling with a window size of $2\times 2$ splits the input into 4 batches (see Fig.~\ref{fig:split-merge}). The purpose of the split pooling is to downsample the input while preserving the whole information. Hence, the pooling window covers each element of the input once, without any gap or overlap.

The merge pooling layer acts as the inverse of split pooling, i.e., it takes the split batches and merges them into one batch. For dense prediction tasks, we apply the same number of merge pooling layers on the final output of network as the number of split pooling layers in a network. The result will be a one-to-one mapping between the input and output pixels of the network.

Formally, given an input $A_{W\times H}$, a split pooling with window size of $w\times h$ splits $A$ into $w*h$ batches
\[\mathcal{B}=\{\,B^{(k,l)}\mid 0\leq k<w \land 0\leq l<h \,\}.\]
The elements of $A$ are assigned to batch $B^{(k,l)}$ as follows
\begin{equation}
\label{eq:split}
b^{(k,l)}_{i,j} = a_{i*w+k,j*h+l}
\end{equation}
for all $0\leq i < W/w$ and $0\leq j<H/h$. The spatial size of each batch will be $W/w\times H/h$. If the input to the split pooling consist of multiple batches, the split pooling splits each input batch separately and returns all resulting splits as a set of batches.

The merge pooling performs the inverse of Eq.~\ref{eq:split}, \ie, if $\mathcal{B}$ is the input and $A$ is the output, the elements of $\mathcal{B}$ are assigned to $A$ as follows
\begin{equation}
\label{eq:merge}
a_{i*w+k,j*h+l} = b^{(k,l)}_{i,j}
\end{equation}
for all $0\leq k<w$, $0\leq l<h$, $0\leq i < W/w$ and $0\leq j<H/h$.

Although preserving the spatial information is beneficial, it increases the number of elements to store during the forward pass.

The advantage of splitting the inputs into batches is that we can distribute the computation of batches on multiple GPUs or processing them sequentially on one GPU. This is ideal for inference. However, in contrast to inference, training phase additionally requires to compute and store the gradients.
This makes the training of very deep networks intractable for large batches of inputs with high resolution images.
We handle this issue by introducing Shrink-Expand pooling layers.

\subsection{Shrink-Expand Pooling}
\label{sec:shrink-expand}
The batch-based design of the split pooling layer makes the forward process of each part (split batch) of the input independent of the other parts (split batches). Furthermore, the one-to-one mapping between the elements of input and output gives a clear path between these elements in both forward and backward direction. Giving these two properties, it is possible to train a network using a subset of split batches produced by each split pooling layer. If we reduce the size of the batch subset to one, the space complexity will be the same as the space complexity of the max pooling and striding convolutions. Also, the training time complexity stays the same.

Giving this reasoning, we introduce the shrink and expand pooling layers, which are the batch-subsampled versions of the split and merge pooling layers. The shrink pooling samples one element at a fixed location within the pooling window, and the corresponding expand pooling uses the same fixed location to perform the reverse of the shrink pooling (see Fig.~\ref{fig:shrink-expand-seq}). Hence, the output of expand pooling is sparse. In other words, the shrink pooling is the same as split pooling except it samples only one of the split batches and returns it.

The sampling location $(i,j)$ is set randomly in each forward pass during training to avoid overfitting to part of the training data. For each pair of shrink and expand pooling in the network, we sample only two numbers $(i,j)$. Fig.~\ref{fig:shrink-expand-seq} shows an example of using a sequence of shrink and expand pooling layers. During training, the expand layers only backpropagate the error of sparse valid elements.

\section{Experiments}
We evaluate the effectiveness of our approach for the semantic segmentation task. To examine the impact of our pooling layer, we modify the ResNet and then compare the performance of the modified version and the original one.

\begin{figure}[!htb]
    \centering
    \includegraphics[width=\linewidth]{gfx/resnet_adapt.pdf}
    \caption{\textbf{Applying split pooling to ResNet.} For applying split pooling to ResNet, we replace the max pooling layer with split pooling (top). Furthermore, we add a split pooling layer after batch normalization layer of convolutional blocks with the stride of $2\times2$ and set the stride to $1\times1$ (bottom).}
    \label{fig:ResNet_adapt}
\end{figure}

\subsection{Experimental Setup}
\paragraph{Baseline FCN32s.} Our experimental models are based on FCN32s \cite{long-shelhamer-fcn-2015} with a variant of ResNet \cite{He_2016_CVPR} as backbone. We adapt the ResNet for semantic segmentation task by removing the average pooling and fully connected layer and replacing them by a $1\times 1$ convolutional layer which maps the output channels of last layer (layer4) to the number of semantic classes. We refer to this model as FCN32s and use it as baseline. FCN32s predictions are 32 times smaller than the input image to the network; thus, the coarse predictions are upsampled to full resolution using bilinear interpolation.

\begin{figure}[!htb]
    \centering
    \begin{tabular}{cc}
    \includegraphics[width=0.23\textwidth]{gfx/bad_gt/img_cs_frankfurt_000001_007857} &
    \includegraphics[width=0.23\textwidth]{gfx/bad_gt/gt_cs_frankfurt_000001_007857} \\
    (a) Image & (b) Ground-truth \\
    \\
    \includegraphics[width=0.23\textwidth]{gfx/bad_gt/8_101_cs_frankfurt_000001_007857} &
    \includegraphics[width=0.23\textwidth]{gfx/bad_gt/sm_101_cs_frankfurt_000001_007857}\\
    (c) FCN8s & (d) Ours
    \end{tabular}
    \caption{\textbf{An example of a slightly inaccurate annotation in Cityscapes dataset.} Although the trunk of the tree in the image (a) is visible through bicycles, it is labeled as bicycle in ground truth (b). Our method can obtain detailed boundaries (d) while FCN8s with max pooling cannot (c). Hence, to fully validate the full potential of our method we need a pixel-accurate ground-truth, such as GTA-5.}
    \label{fig:cs-inaccurate}
\end{figure}

\paragraph{SMP-\{18, 34, 101\}.} As illustrated in Fig.~\ref{fig:ResNet_adapt}, in order to apply SMP to ResNet backbone of FCN32s, we simply replace the max pooling layer with a split pooling layer and add a split pooling layer after the batch normalization layer of convolutional blocks (Conv-BN-ReLU) with a stride of $2\times 2$ and set their strides to $1\times 1$. In our experiments, the window size of split pooling layers is $2\times 2$. We call the resulting model, according to the type of backbone ResNet, SMP-18, SMP-34 and SMP-101. The output of SMP-X consist of 1024 batches, since we always have 5 SMP layers each giving 4 batches. We merge the output split batches by performing merge pooling five times. The final result has the same resolution as the input. During training phase, we replace split and merge pooling layer by shrink and expand layers.

\paragraph{FCN8s.} Furthermore, we compare our models to the FCN8s model with original ResNet backbone. The main difference between FCN32s and FCN8s is that FCN8s has two extra $1\times 1$ convolution layers to map the features channel of layer2 and layer3 outputs of ResNet to the number of semantic classes. Then, the output of the network and these new convolutions are resized to the same size and are added together. The final output of FCN8s is 8 times smaller than the input image and should be upsampled to full resolution using bilinear interpolation.

\subsection{Implementation Details}
\paragraph{Data Augmentation.} For all the models we used random crop of size $512\times 512$ and horizontal flip data augmentation. We used a fixed seed for data augmentation for all the experiments.

\paragraph{Training.} We initialize the ResNet backbones with pretrained models from Imagenet. We optimize the parameters using Adam solver \cite{KingmaB14} with learning rate of $1e-5$ and weight decay of $5e-4$. We use the batch size of 10 for all the experiments.

\setlength{\tabcolsep}{4pt}
\begin{table}[!htb]
    \centering
    \begin{tabular}{llccc}
        \hline
         & Backbone &  Pooling & IoU & Size \\
        \hline \hline
        FCN32s & ResNet34 & MP & 64.5 &  21.29M\\
        FCN8s  & ResNet34 & MP & 67.8 &  21.30M\\
        \hline
        SMP-34 (ours) & ResNet34 & SMP & \textbf{68.8} & 21.29M\\
        \hline
        \hline
        FCN32s & ResNet101 & MP & 65.5 & 42.54M\\
        FCN8s  & ResNet101 & MP & 69.1 & 42.56M\\
        \hline
        SMP-101 (ours) & ResNet101 & SMP & \textbf{69.2} & 42.54M\\
        \hline
    \end{tabular}
    \vspace{0.1cm}
    \caption{\textbf{Cityscapes quantitative results.} Please note that the architecture of SMP-X networks is the same as the FCN32s. The FCN8s architecture has two extra $1\times 1$ convolutions.}
    \label{tab:cs}
\end{table}

\setlength{\tabcolsep}{4pt}
\begin{table}[!htb]
    \centering
    \begin{tabular}{llccc}
        \hline
         & Backbone & Pooling & IoU & Size\\
        \hline \hline
        FCN32s & ResNet18 & MP & 71.1 & 11.186M\\
        FCN8s  & ResNet18 & MP & 74.2 & 11.193M\\
        \hline
        SMP-18 (ours) & ResNet18 & SMP & \textbf{77.5} & 11.186M\\
        \hline \hline
        FCN32s & ResNet34 & MP & 73.1 & 21.29M\\
        FCN8s  & ResNet34 & MP & 76.7 & 21.30M\\
        \hline
        SMP-34 (ours) & ResNet34 & SMP & \textbf{80.2} & 21.29M\\
        \hline
        \hline
        FCN32s & ResNet101 & MP & 74.6 & 42.54M\\
        FCN8s  & ResNet101 & MP & 76.7 & 42.56M\\
        \hline
        SMP-101 (ours) & ResNet101 & SMP & \textbf{80.3} & 42.54M\\
        \hline
    \end{tabular}
    \vspace{0.1cm}
    \caption{\textbf{Quantitative results for GTA-5.} Please note that the architecture of SMP-X networks is the same as the FCN32s, while the FCN8s has two extra $1\times 1$ convolutions.}
    \label{tab:gta}
\end{table}

\setlength{\tabcolsep}{2.5pt}
\begin{table}[!htb]
    \centering
    \begin{tabular}{l ccccccc c}

         & \rot{pole}& \rot{traffic light}& \rot{traffic sign}& \rot{person}& \rot{rider}& \rot{motorcycle}& \rot{bicycle} & IoU  \\
         \hline
         FCN32s-Res34 & 35.4 & 51.5 & 63.0 & 70.4 & 50.3 & 51.8 & 68.4 & 55.8 \\
         FCN8s-Res34 & 53.5 & 57.6 & 69.9 & 77.1 & 53.6 & 51.0 & 72.2 & 62.1 \\
         \hline
         SMP-34(ours) & \textbf{60.5} & \textbf{63.7} & \textbf{74.5} & \textbf{78.8} & \textbf{54.1} & \textbf{52.9} & \textbf{73.9} & \textbf{65.5}  \\
         \hline
         \hline
         FCN32s-Res101 & 39.2 & 58.1 & 66.9 & 71.9 & 51.4 & 53.9 & 70.6 & 58.9 \\
         FCN8s-Res101 & 56.1 & 63.1 & 72.2 & 78.2 & 55.4 & 52.8 & 74.7 & 64.6 \\
         \hline
         SMP-101(ours) & \textbf{63.3} & \textbf{68.7} & \textbf{75.8} & \textbf{79.9} & \textbf{56.1} & \textbf{52.9} & \textbf{76.6} & \textbf{67.6}\\
         \hline
    \end{tabular}
    \vspace{0.1cm}
    \caption{\textbf{Performance on Cityscapes for small and thin objects.}}
    \label{tab:cs-classes}
\end{table}

\setlength{\tabcolsep}{2.5pt}
\begin{table}[!htb]
    \vspace{-0.3cm}
    \centering
    \begin{tabular}{l ccccccccccccccccccc c}

         & \rot{pole}& \rot{traffic light}& \rot{traffic sign}& \rot{person}& \rot{rider}& \rot{motorcycle}& \rot{bicycle} & IoU  \\
         \hline
         FCN32s-Res18 & 44.6 & 45.3 & 60.7 & 69.3 & 64.5 & 59.8 & 38.5 & 54.7\\
         FCN8s-Res18 & 54.0 & 52.7 & 63.8 & 74.4 & \textbf{70.4} & 66.5 & 40.3 & 60.3\\
         \hline
         SMP-18(ours) & \textbf{72.7} & \textbf{69.2} & \textbf{73.6} & \textbf{76.8} & 67.8 & \textbf{69.3} & \textbf{46.5} & \textbf{68.0} \\
         \hline
         \hline
         FCN32s-Res34 & 46.0 & 47.7 & 63.0 & 70.0 & 67.9 & 63.9 & 49.5 & 58.3 \\
         FCN8s-Res34 & 55.6 & 57.2 & 68.3 & 76.3 & \textbf{74.0} & 65.8 & 51.3 & 64.1 \\
         \hline
         SMP-34(ours) & \textbf{74.3} & \textbf{71.4} & \textbf{73.7} & \textbf{81.0} & 70.5 & \textbf{73.7} & \textbf{58.5} & \textbf{71.9}\\
         \hline
         \hline
         FCN32s-Res101 & 47.3 & 50.3 & 68.7 & 70.3 & 63.9 & 66.2 & 58.1 & 60.7\\
         FCN8s-Res101 & 57.1 & 59.0 & 70.8 & 75.5 & 67.1 & \textbf{70.6} & 62.4 & 66.1\\
         \hline
         SMP-101(ours) & \textbf{76.8} & \textbf{74.9} & \textbf{78.9} & \textbf{79.0} & \textbf{69.2} & 69.6 & \textbf{67.8} & \textbf{73.7}\\
         \hline
    \end{tabular}
    \vspace{0.1cm}
    \caption{\textbf{Performance on GTA-5 for small and thin objects}}
    \label{tab:gta-classes}
\end{table}

\subsection{Cityscapes}
The Cityscapes dataset \cite{Cordts2016Cityscapes} consists of images with the size of $2048\times 1024$. The images are annotated with 19 semantic classes. We evaluate the performance of our SMP-34 and SMP-101 on Cityascapes validation set and compare it with FCN32s and FCN8s with ResNet34 and ResNet101 backbone networks. The full size images are used for evaluation, \ie without downsampling or cropping.

\setlength{\tabcolsep}{2pt}
\begin{figure*}[!htp]
    \centering
    \begin{tabular}{lc}
        & \begin{tabularx}{\textwidth}{>{\centering\arraybackslash}X>{\centering\arraybackslash}X>{\centering\arraybackslash}X>{\centering\arraybackslash}X}
        FCN32s & FCN8s & SMP (ours) & \shortstack{Image / \\ Ground Truth}
        \end{tabularx}\\
    \rot{\hspace{0.3cm}ResNet101 \hspace{0.5cm}ResNet34} & \includegraphics[width=0.96\linewidth]{gfx/cs_qual/cs_frankfurt_000001_014406_crop} \\
    \\
    \rot{\hspace{0.3cm}ResNet101 \hspace{0.5cm}ResNet34} & \includegraphics[width=0.96\linewidth]{gfx/cs_qual/cs_munster_000152_000019_crop} \\
    \\
    \rot{\hspace{0.3cm}ResNet101 \hspace{0.5cm}ResNet34} & \includegraphics[width=0.96\linewidth]{gfx/cs_qual/cs_frankfurt_000001_011835_crop}\\
    \\
    \rot{\hspace{0.3cm}ResNet101 \hspace{0.5cm}ResNet34} & \includegraphics[width=0.96\linewidth]{gfx/cs_qual/cs_frankfurt_000001_007285_crop} \\
    \end{tabular}
    \caption{\textbf{Cityscapes qualitative results.} In each block, the top row is related to models with ResNet34 backbone and the bottom row to ResNet101. The last column of each block shows the input image (top) and the ground-truth (bottom). To enhance the visual comparison of the results, we have cropped the output labelling. Further results, also in full resolution, can be found in supplementary material.\vspace{1cm}}
    \label{fig:cs-qual}
\end{figure*}

\begin{figure*}[!htp]
    \centering
    \begin{tabular}{cc}
        & \begin{tabularx}{\textwidth}{>{\centering\arraybackslash}X>{\centering\arraybackslash}X>{\centering\arraybackslash}X>{\centering\arraybackslash}X}
        FCN32s & FCN8s & SMP (ours) & \shortstack{Image / \\ Ground Truth}
        \end{tabularx}\\
    \rot{\hspace{0.3cm}ResNet101 \hspace{0.6cm}ResNet34 \hspace{0.7cm}ResNet18} & \includegraphics[width=0.96\linewidth]{gfx/gta_qual/gta_16087_crop} \\
    \\
    \rot{\hspace{0.3cm}ResNet101 \hspace{0.6cm}ResNet34 \hspace{0.7cm}ResNet18} & \includegraphics[width=0.96\linewidth]{gfx/gta_qual/gta_14989_crop} \\
    \\
    \rot{\hspace{0.3cm}ResNet101 \hspace{0.6cm}ResNet34 \hspace{0.7cm}ResNet18} & \includegraphics[width=0.96\linewidth]{gfx/gta_qual/gta_16335_crop}
    \end{tabular}
    \caption{\textbf{GTA-5 qualitative results.} In each block, the top row is related to ResNet18 backbone, middle row to ResNet34, and bottom row to ResNet101. The last column of each block shows the input image (top) and the ground-truth (bottom). To enhance the visual comparison of the results, we have cropped the output labelling. Further results, also in full resolution, can be found in supplementary material.}
    \label{fig:gta-qual}
\end{figure*}

Table~\ref{tab:cs} summarizes the comparison of different methods with respect to their setups, performances and number of parameters (size). Each SMP model outperforms the corresponding FCN32s and FCN8s with the same ResNet backbone significantly. SMP-34 even outperforms FCN32s-ResNet101 by $3.3\%$, although the latter is almost 2 times larger.
The performance of our approach on objects from small and thin classes is reported in Table~\ref{tab:cs-classes}. As we can see in the table, SMP models outperforms their corresponding original models on small and thin objects. Both SMP models outperform their corresponding FCN32s models for \textit{pole} by $24\%$, for \textit{traffic light} by more than $10\%$, for \textit{traffic sign} by more than $9\%$, and for  \textit{person} by $8\%$. As it is shown in Fig.~\ref{fig:cs-qual}, the improvement of these classes can be noticed visually as well. The performance for the remaining objects is mostly better, or sometimes slightly worse.

Qualitative results are shown in Fig.~\ref{fig:cs-qual}. For the sake of improved visibility we have cropped the results. The full size output images can be found in supplementary material.

\subsection{GTA-5}

Cityscapes is one of the most accurately annotated semantic segmentation datasets, however it is still not pixel-accurate (see Fig.~\ref{fig:cs-inaccurate}). Obtaining pixel-accurate annotations from real data is extremely challenging and expensive. Therefore, for analysing the full potential of our method, we evaluate our method on GTA-5 \cite{Richter_2016_ECCV}, which is a synthetic dataset with the same semantic classes as Cityscapes. Since GTA-5 is a synthetic, the annotations are pixel-accurate and ideal for our purpose.
The GTA-5 dataset \cite{Richter_2016_ECCV} consist of 24,999 realistic synthetic images with pixel-accurate semantic annotations. We randomly select 500 images as validation set, which we did not use for training.

Table~\ref{tab:cs} summarizes the comparison of different methods with respect to their setups, performances and number of parameters (size). As we can see, due to the pixel-accurate annotations of GTA-5 dataset, the improvement of our proposed models, over their baselines, is more significant compared to Cityscapes. Each SMP model outperforms the corresponding FCN32s and FCN8s with the same ResNet backbone significantly. Particularly, our SMP-18 even outperforms its FCN32s-ResNet101 and FCN8s-ResNet101 counterparts, although it has 4 times fewer parameters.
The performance of our approach on objects from small and thin classes is reported in Table~\ref{tab:gta-classes}. As we can see, similarly to Cityscapes, SMP models outperforms their corresponding original models on small and thin objects. Compared to FCN32s models, our corresponding SMP models improve the categories \textit{pole} by more than $28\%$, \textit{traffic light} by more than $23\%$, \textit{traffic sign} by more than $10\%$, and \textit{person} by more than $7\%$. The improvement over these classes is also visually significant (see Fig.~\ref{fig:gta-qual}).

\subsection{Run-time Analysis}
For analyzing the time complexity of the Split-Merge pooling, we designed small networks to just focus on the proposed pooling layers instead of analyzing the time complexity of them on a particular task with specific network architecture. As it is shown in Fig.~\ref{fig:run-time-archs}, we consider three networks with (a) max pooling, (b) dilated convolution, and (c) Split pooling.
We choose to compare our proposed pooling setup (c) with dilated convolutions (b) due to the success of the dilated convolution in dense prediction tasks. Almost all state-of-the-art approaches in dense prediction tasks (such as semantic segmentation, depth estimation, optical flow estimation) are using dilated convolutions in their architectures to achieve a detailed output.

\begin{figure}[!htb]
	\centering
	\includegraphics[width=\linewidth]{gfx/runtime-setup}
	\caption[Runtime analysis setup]{Architectures used for runtime analysis. \textit{conv2} in (c) is identical to \textit{conv2} in (a) while \textit{conv2} in (c) is dilated convolution with padding $2$.}
	\label{fig:run-time-archs}
\end{figure}

In Table~\ref{tab:run-time}, we show the Giga floating-point operation (GFLOP) of each component of each setup for an input tensor of size $1\times3\times256\times256$. As we can see, the dilated convolution setup and split pooling setup have the same GFLOPs which means dilated convolution layers can be replaced with our proposed pooling layers in an arbitrary architecture without changing the complexity of the network. However, our split pooling layer has two advantages:
\begin{enumerate}
	\item faster training time using shrink-expand layers
	\item faster inference time by parallelizing the forward-computation of split layer output batches (in this example setup computation of \textit{conv2}, see Table~\ref{tab:run-time})
\end{enumerate}

\begin{table}[!htb]
	\centering
	\begin{tabular}{l|cc|c|cc||c}
		& batches & conv1 & pooling & batches & conv2 & total \\
		\hline
		(a) Max Pooling   & 1 & 0.23 & 0 & 1 & 2.42 & 2.65 \\
		(b) Dilated Conv. & 1 & 0.23 & - & 1 & 9.68 & 9.92 \\
		(c) Split Pooling & 1 & 0.23 & 0 & 4 & 9.68 & 9.92 \\
		\hline
	\end{tabular}
	\vspace{0.1cm}
	\caption[Runtime analysis]{GFLOPs of the models calculated on the input size of $1\times3\times256\times256$. Note that the number of batches are increased after split pooling.}
	\label{tab:run-time}
\end{table}

\vspace{-0.2cm}
\section{Conclusion}
We proposed a novel pooling method SMP with the goal of preserving the spatial information throughout the entire network. SMP can be used instead of any subsampling operations in a network architecture. We show that by replacing subsampling operations with SMP in ResNet, we achieved two important properties for any dense prediction task at the same time: i) the network has a large receptive field, ii) the network provides a unique mapping from input pixels to output pixels. Furthermore, the computation of a network with SMP can be distributed to multiple GPUs due to batch-based design of SMP.
We show experimentally that the resulting network outperforms the original one significantly.

\paragraph{}
\textbf{Acknowledgements} This work is funded by the DFG grant “COVMAP: Intelligente Karten mittels gemeinsamer GPS- und Videodatenanalyse” (RO 4804/2-1).

{\small
\bibliographystyle{ieee_fullname}


\section{Introduction}
\label{sec:intro}
With the advantage of low-cost and simplicity in design, more and more rolling-shutter (RS) CMOS cameras \cite{Meingast_Geometric_2005} have been used in various real-world computer vision applications. Different from their global shutter (GS) counterparts that expose all pixels simultaneously, RS cameras capture the pixels of each row (or column) at consecutive times. Therefore, when there is instantaneous motion between the scene and the camera, the captured image will show rolling shutter effects such as skew or wobble, and the level of rolling shutter effect varies with respect to the camera motion and the object motion.

\begin{figure}
\begin{minipage}[b]{0.6\linewidth}
  \centering
  \centerline{\includegraphics[width=7.25cm]{Fig-1-2-1.pdf}}
\end{minipage}
\centering
\vspace{-0.4cm}
\caption{Overview of our proposed method. From consecutive stereo RS images, we solve for the relative pose via our proposed linear algorithm. With the estimated camera motion and depth map by SGM, we could rectify the RS image.}
\vspace{-0.5cm}
\end{figure}

Currently, most of the geometric algorithms in computer vision are designed under the perspective camera model (\emph{i.e.} GS camera model). However, this model fails to handle the RS images due to the significant difference between RS and GS cameras. Therefore, geometric algorithms for RS cameras have attracted more and more attention. In \cite{Meingast_Geometric_2005}, Meingast \emph{et al} introduced the geometric model for RS cameras. Along this direction, many algorithms have been proposed to solve the corresponding multi-view geometric problems. Purkait \emph{et al}. \cite{Purkait_manhattan_2017} performed pose estimation under the Manhattan world assumption.
Oliver \emph{et al}. \cite{Saurer_RS_stereo_2013} proposed an algorithm based on plane sweeping to recover the depth map and 3D scene for multi-view RS images. Fraundorfer \emph{et al}. \cite{Fraundorfer_Minimal_2010} fused information from external sensors like IMU to estimate the relative motion while Lee \emph{et al} \cite{Lee_Gyroscope_2019} simplified the pose estimation problem by using the gyroscope.
Dai \emph{et al}. \cite{Dai_epipolar_2016} derived a 20-point linear algorithm and a 44-point linear algorithm under different RS models based on generalized epipolar geometry.

It should be noted that all of these works are designed to solve the relative pose of a monocular RS camera. In practise, stereo cameras provide extra benefits over its monocular counterpart, \emph{i.e.,} more accurate motion estimation, depth estimation and the availability of global scale. Existing relative pose estimation methods for stereo cameras such as \cite{Vakhitov_2018_line_point} are all designed for GS cameras.
Surprisingly, no previous attempt has been reported on solving the {\bf 6 DoF relative pose} problem with stereo RS cameras.
Motivated by relative pose estimation for monocular camera \cite{Dai_epipolar_2016}, we exploit the epipolar geometry to solve the relative pose from consecutive RS frames and present the first relative pose estimation algorithm for stereo RS cameras as shown in Fig.1.
Specifically, we figure out a linear solver for 6 DoF relative pose from consecutive stereo images under a continuous motion model. To solve the pose, 9 pairs of correspondences from the left and right cameras respectively are needed. The pose of each scanline on each frame can be solved under the assumption of a constant velocity motion model.

The estimated 6 DoF camera motion also enables the correction of the RS effects, which generally requires the underlying scene structures and the camera motions between scanlines or between views. Existing work with monocular cameras depend on strong assumptions to estimate the scene structure.
Lao \emph{et al}. \cite{Lao_2018_Robust} used at least 4 curves to estimate the camera's rotation speed to eliminate the RS effect.
Zhuang \emph{et al}.\cite{Zhuang_Rolling_shutter_aware_2017} introduced optical flow based on RS successive frames, which proposed an 8-point algorithm for constant motion, and 9-point algorithm for constant acceleration motion, then they achieved SfM as well as image rectification.
Vasu \emph{et al}. \cite{Vasu_Occlusion_2018} realized that the RS distortion effect is depth-dependent, and proposed an algorithm that can recover 3D scene from a set of RS distorted images.

Under our stereo RS camera configuration, we are able to estimate the scene structure without assumptions. Specifically, we show that for stereo RS images, using the semi-global matching (SGM) algorithm which is based on stereo GS image can also accurately estimate the depth map, then the depth will be used in the RS image correction to obtain a smoothly corrected image.
Experimental results on both simulated points and synthetic RS images prove the effectiveness of our relative pose algorithm and its application in RS image correction,

\section{METHOD}
\label{sec:format}
In this section, we establish the model of stereo RS cameras under consecutive frames and derive the linear solver.

\subsection{Rolling shutter camera model}
\label{ssec:subhead}
First of all, we will illustrate the differences between the GS model and the RS model in the form of image formation. GS cameras expose all pixels simultaneously, thus the projection can be expressed as:
$\lambda_{i} \mathbf{x}_{i} = \boldsymbol{\rm K}(\boldsymbol{\rm R} \mathbf{X}_{i} + \mathbf{T})$, where $\boldsymbol{\rm R}$ and $\mathbf{T}$ represent the rotation and translation of the camera in the world coordinate respectively, $\boldsymbol{\rm K}$ denotes the intrinsic calibration matrix, $\mathbf{X}_{i}$ is the 3D point which is projected to $\mathbf{x}_{i}=[u_{i},v_{i},1]^{T}$ in image plane with a depth value $\lambda_{i}$.

For RS camera, each scanline will correspond to different poses when the camera is moving during image acquisition. We assume that the camera scans in rows. So both $\boldsymbol{\rm R}$ and $\boldsymbol{\rm T}$ can be described as functions of image row $u_{i}$. The projection process under RS camera model is expressed as:
\begin{align}
\lambda_{i} \mathbf{x}_{i}= \lambda_{i}[u_{i},v_{i},1]^{T}=
\boldsymbol{\rm K}[ \boldsymbol{\rm R}(u_{i}) \mathbf{X}_{i} + \mathbf{T}(u_{i})].
\end{align}
Next, we derive our linear solver of the stereo RS cameras under consecutive frames.

\subsection{Stereo RS algorithm under consecutive images}
\label{ssec:subhead}
As illustrated in Fig.1(b), we use the continuous camera motion model in this paper.
For clarity, we denote the left and right camera at the first moment and the second moment as ${\rm I}_{1}$, ${\rm I}_{2}$, ${\rm I}_{3}$, ${\rm I}_{4}$ respectively.

We suppose that stereo cameras undergo constant velocity motion among consecutive frames. As shown in Fig.2, we define the total readout time of one frame as $T_{a}$ and the delay time between consecutive frames as $T_{b}$. The stereo cameras are rigidly connected through the baseline, so $T_{a}$ and $T_{b}$ of the left and right cameras are identical. $\{\mathbf{w}_{1},\mathbf{w}_{2},\mathbf{w}_{3},\mathbf{w}_{4}\}\in so(3)$ and $\{\mathbf{d}_{1},\mathbf{d}_{2},\mathbf{d}_{3},\mathbf{d}_{4}\}\in \mathbb{R}^{3}$ represent the rotation and translation velocities of the first row on ${\mathrm{I}_{1}}$, ${\mathrm{I}}_{2}$, ${\mathrm{I}}_{3}$ and ${\mathrm {I}}_{4}$. Furthermore, we assume that the stereo cameras move around the center of the baseline, so we set pose of the center of baseline as ${\mathbf{w}_{0} = \mathbf{0}, \mathbf{d}_{0}=\mathbf{0}}$, and the length of baseline is set as $2 \mathbf{b}$. Therefore, the pose of first row on ${\mathrm{I}}_{1}$ can be represented as ${\mathbf{w}_{1} = \mathbf{0}, \mathbf{d}_{1}=\mathbf{b}}$, and the pose of first row on ${\rm I}_{2}$ is ${\mathbf{w}_{2} = \mathbf{0}, \mathbf{d}_{2}=-\mathbf{b}}$.
\begin{figure}[htb]
\begin{minipage}[b]{0.7\linewidth}
  \centering
  \vspace{-0.3cm}
  \centerline{\includegraphics[width=7.3cm]{consecutive_RS.pdf}}
\end{minipage}
\centering
\vspace{-0.5cm}
\caption{Illustration of two consecutive frames of stereo RS cameras.}
\vspace{-0.3cm}
\end{figure}

Suppose that the relative motions between consecutive frames are $\mathbf{w} = [w_{1},w_{2},w_{3}]^{T}$ and $\mathbf{d} = [{\rm d}_{1},{\rm d}_{2},d_{3}]^{T}$, the rotation and translation velocities of each row on each frame can be computed through linear interpolation.
The pose of row $u_{i}$ on ${\rm I}_{1}$ is represented as:
\begin{align}
\mathbf{w}(u_{1}) = hu_{1}\mathbf{w},  \quad   \mathbf{d}(u_{1}) = hu_{1}\mathbf{d} + \mathbf{b},
\end{align}
where $h = \varphi /N$, $\varphi = T_{a} /(T_{a}+T_{b})$ is the readout time ratio, $N$ is the total number of rows per frame. Similar to ${\rm I}_{1}$, the pose of row $u_{3}$ on ${\rm I}_{3}$ is expressed as:
\begin{align}
&\mathbf{w}(u_{3}) = (1+hu_{3})\mathbf{w},  \quad   \mathbf{d}(u_{3}) = (1+hu_{3})\mathbf{d} +\mathbf{b}.
\end{align}
Same as before, the pose of row $u_{2}$ on ${\rm I}_{2}$ and pose of row $u_{4}$ on ${\rm I}_{4}$ are:
\begin{align}
&\mathbf{w}(u_{2}) = hu_{2}\mathbf{w},  \quad   \mathbf{d}(u_{2}) = hu_{2}\mathbf{d} - \mathbf{b},
\end{align}
and
\begin{align}
\mathbf{w}(u_{4}) = (1+hu_{4})\mathbf{w},  \quad   \mathbf{d}(u_{4}) = (1+hu_{4})\mathbf{d} - \mathbf{b}.
\end{align}
Assume that the relative rotation between consecutive frames is small, the mapping between the rotation matrix $\boldsymbol{\rm R}$ and the rotation vector $\boldsymbol{w}$ can be approximated as: $\boldsymbol{\rm R}  \simeq {\rm exp}(\mathbf{w}) =\boldsymbol{\rm I}+[\mathbf{w}]_{\times}$. Therefore, we can develop the essential matrix between ${\rm I_{1}}$ and ${\rm I_{3}}$ as:
\begin{subequations}
\begin{align}
&\boldsymbol{\rm E}_{u_{1},u_{3}} = \left[\mathbf{t}_{u_{1},u_{3}}\right]_{\times} \boldsymbol{\rm R}_{u_{1},u_{3}},\\
&\boldsymbol{\rm R}_{u_{1},u_{3}} = \boldsymbol{\rm R}_{u_{3}}\boldsymbol{\rm R}_{u_{1}}^{T}= \boldsymbol{\rm I}+(1+hu_{3}-hu_{1})[\mathbf{w}]_{\times}, \\
&\mathbf{T}_{u_{1},u_{3}} = \mathbf{d}(u_{3})- \boldsymbol{\rm R}_{u_{1},u_{3}} \mathbf{d}(u_{1}).
\end{align}
\end{subequations}

\noindent Consequently, we can reorganize the formulation to get:
\begin{align}
&\boldsymbol{\rm E}_{u_{1},u_{3}} = a_{1}[\mathbf{d}]_{\times}-a_{1}[[\mathbf{w}]_{\times} \mathbf{b}]_{\times}+{{a_{1}}^{2}}[\mathbf{d}]_{\times}[\mathbf{w}]_{\times}-
\\
\notag
&a_{1}k_{1}[[\mathbf{w}]_{\times}\mathbf{d}]_{\times}-{{a_{1}}^{2}}[[\mathbf{w}]_{\times}\mathbf{b}]_{\times}[\mathbf{w}]_{\times}-{{a_{1}}^{2}}k_{1}[[\mathbf{\{w}]_{\times}\mathbf{d}]_{\times}[\boldsymbol{w}]_{\times}.
\end{align}
where $a_{1} = k_{3}-k_{1}$,$k_{1} = hu_{1}$, $k_{3}=1+hu_{3}$. Similarly, we can also sort out the essential matrix between ${\rm I_{2}}$ and ${\rm I_{4}}$ as:
\begin{align}
&\boldsymbol{\rm E}_{u_{2},u_{4}} = a_{2}[\mathbf{d}]_{\times}+a_{2}[[\mathbf{w}]_{\times}\mathbf{b}]_{\times}+{{a_{2}}^{2}}[\mathbf{d}]_{\times}[\mathbf{w}]_{\times}-
\\
\notag
&a_{2}k_{2}[[\mathbf{w}]_{\times}\mathbf{d}]_{\times}+{{a_{2}}^{2}}[[\mathbf{w}]_{\times}\mathbf{b}]_{\times}[\mathbf{w}]_{\times}-{{a_{2}}^{2}}k_{2}[[\mathbf{w}]_{\times}\mathbf{d}]_{\times}[\boldsymbol{w}]_{\times},
\end{align}
where $a_{2} = k_{4}-k_{2}$,$k_{2}=hu_{2}$, $k_{4}=1+hu_{4}$. The epipolar constraints are expressed as:
\begin{subequations}
\begin{align}
&[u_{3},v_{3},1]\boldsymbol{\rm E}_{u_{1},u_{3}}[u_{1},v_{1},1]^{{T}}=0,\\
&[u_{4},v_{4},1]\boldsymbol{\rm E}_{u_{2},u_{4}}[u_{2},v_{2},1]^{{T}}=0.
\end{align}
\end{subequations}
By formulating Eq.(9a) and (9b) into the form of a linear equation, we can get $\boldsymbol{\rm A}_{1} \mathbf{X}=0$ and $\boldsymbol{\rm A}_{2}\mathbf{X}=0$ respectively. Combining them together yields a linear equation:
\begin{align}
\left[\begin{array}{c} \boldsymbol{\rm A}_{1}\\ \boldsymbol{\rm A}_{2}\end{array}\right] \mathbf{X}=\boldsymbol{\rm A} \mathbf{X}=0,
\end{align}
where $\mathbf{X} =[d_{1},d_{2},d_{3},w_{1},w_{3},E_{11},E_{12},E_{13},E_{21},E_{22},E_{23},\\
E_{31},E_{32},E_{33},w_{1}w_{1},w_{1}w_{2},w_{1}w_{3},w_{2}w_{3},w_{3}w_{3},S_{1},S_{2},S_{3},\\
S_{4},S_{5},S_{6},S_{7}.S_{8}]^{T}$ in which each of parameter is composed by $\mathbf{w}$ and $\mathbf{d}$, and $\boldsymbol{\rm A}_{1},\boldsymbol{\rm A}_{2}$ are matrices consisted of known parameter $k_{1},k_{2},k_{3},k_{4},b$ and coordinates of correspondences points. In Eq.(10), $\mathbf{X}$ has 27 unknowns. However, we notice that there are higher-order terms in $\boldsymbol{\rm A}$, these terms will not affect the results so we omit them. Then the vector $\mathbf{X}$ contain 19 parameters, thus a linear 18-point solver must exist to solve this formulation, \emph{i.e.} 9 pairs of correspondences for ${\rm I}_{1}$ and ${\rm I}_{3}$ and 9 pairs of correspondences for ${\rm I}_{2}$ and ${\rm I}_{4}$. Next, the right singular vector corresponding to the smallest singular value of $\boldsymbol{\rm A}$ is the solution of Eq.(10). Then, $\mathbf{w}$ and $\mathbf{d}$ can be extracted from $\mathbf{X}$.

\subsection{Rolling Shutter Correction}
\label{ssec:subhead}
In this section, we employ the Semi-Global Matching (SGM) algorithm \cite{Hirschmuller_SGM_2007} to estimate the depth map from the stereo image pair, and combine the camera pose of each row obtained in the previous section, the correction of the RS images can be completed. First, we use the pose and depth of each pixel in the RS image to back-project each pixel into 3D space, and then reproject 3D points to the image plane corresponding to the pose of the first row as follow:
\begin{subequations}
\begin{align}
&\mathbf{X}_{i} = \boldsymbol{\rm R}_{u_{i}}^{T}[\lambda_{i} {\boldsymbol{\rm K}}^{-1} {[u_{i},v_{i},1]}^{T} - \mathbf{T}_{u_{i}}],\\
&{\mathbf{x}_{i}}^{'}= \boldsymbol{\rm K}[ \boldsymbol{\rm R}_{0} \mathbf{X}_{i} + \mathbf{T}_{0}].
\end{align}
\end{subequations}
We figure out that for stereo RS images, a relatively accurate depth map can be obtained by the SGM algorithm, the depth can be further applied to the correction to get a smoothly undistorted image. Note that the stereo RS cameras enable both 6 DoF camera motion estimation and accurate depth map estimation, which  results in our improved RS correction.

\section{EXPERIMENTS}
\label{sec:pagestyle}

In this section, we evaluate our algorithm on simulated points and synthetic RS images. Note that $e_{T}={\rm {acos}}(\mathbf{d}_{est}^{T} \mathbf{d}_{gt}/(\Arrowvert \mathbf{d}_{gt} \Arrowvert \\ \Arrowvert \mathbf{d}_{est} \Arrowvert ))$ and $e_{R}={\rm {acos}}(({\rm Tr}(
\boldsymbol{\rm R}_{est}\boldsymbol{\rm R}_{gt}^{T})-1)/2)$ are used to measure the translation and rotation errors, respectively.

\subsection{Simulation Experiments}
\label{ssec:subhead}
To verify the effectiveness of our proposed algorithm, we simulate matching points that satisfy the RS projection model. In the simulation experiments, we set the size of the image as $900\times900$ with an 810 focal length. Then, we give the value of relative translation $\mathbf{d}$ and relative rotation $\mathbf{w}$. The baseline length $\mathbf{b}$ and the readout time ratio $\varphi$ are also set. When all parameters are known, the correspondences between ${\rm I}_{1}$ and ${\rm I}_{2}$ and correspondences between ${\rm I}_{3}$ and ${\rm I}_{4}$ can be generated. To ensure the creditability, each experimental value is the mean value of 300 repetitions. Moreover, in order to check the robustness and practicability of the algorithm, we do not use RANSAC or nonlinear optimization in the simulation experiments, and experimental results are compared with the experimental results of the GS algorithm \cite{Hartley_Denfense_1997}.

\noindent\textbf{Accuracy versus noise level.} First, we add different levels of random Gaussian noise to the image coordinates of correspondences to evaluate the robustness of our algorithm. Statistical results are shown in Fig.3, illustrating that although the estimation error increases with the increasing noise level, the overall increase is tiny.
\begin{figure*}
\centering
\includegraphics[width=3.7cm]{noise_level1.pdf}
\hspace{0.1ex}
\includegraphics[width=4cm]{Trans_err_t.pdf}
\hspace{0.1ex}
\includegraphics[width=3.9cm]{Trans_err_w.pdf}
\hspace{0.1ex}
\includegraphics[width=3.9cm]{Trans_err_ratio.pdf}

\vspace{-0.1cm}
\setlength{\abovecaptionskip}{0.cm}
\setlength{\belowcaptionskip}{-0.2cm}
\hspace{-0.8ex}
\centering
\includegraphics[width=3.9cm]{baseline_test.pdf}
\hspace{-0ex}
\includegraphics[width=4cm]{Rots_err_t.pdf}
\hspace{-0.5ex}
\includegraphics[width=4cm]{Rots_err_w.pdf}
\hspace{-0.5ex}
\includegraphics[width=4cm]{Rots_err_ratio.pdf}
\vspace{-0.3cm}
\caption{Quantitative evaluation for our stereo RS algorithm on simulation experiments. Results on the first column show the effect of our algorithm under different noise levels and the results with different length of baseline. The second and third columns test the estimation accuracy at different translation and rotation velocity respectively. And the last columns evaluate the effect of different readout time ratio.}
\vspace{-0.5cm}
\end{figure*}

\noindent\textbf{Accuracy versus RS velocity.} Then, we test how RS velocity will affect the accuracy of our algorithm. For a constant Gaussian noise level of $10^{-3}$, we test the performance of the increasing rotation and increasing translation respectively as shown in the second and third column in Fig.3.
It can be observed that our RS algorithm always achieves higher accuracy than the GS algorithm not only in translation estimation but also in rotation estimation with the increasing velocity. Our RS algorithm can maintain more stable and accurate results.

\noindent\textbf{Accuracy versus readout time radio.} To test and observe how the readout time ratio affects the estimation result, we vary the value of readout time ratio $\varphi$ from 0.2 to 1 according to \cite{Im_High_2015}. Simultaneously, the value of $\mathbf{d}$ and $\mathbf{w}$ are kept constant with a fixed Gaussian noise level of $10^{-3}$. Results are shown in the last column in Fig.3, the estimation error of both the GS algorithm and RS algorithm increase with the increase of $\varphi$, but GS gets higher errors. The reason for this result can be interpreted as that the RS effect becomes more evident with increasing $\varphi$.

\noindent\textbf{Accuracy versus baseline length.} Finally, we investigate how the length of the baseline can affect performances. We keep the value of $\mathbf{d}$ and $\mathbf{w}$ fixed, then increase the baseline from 0.1 to 1. As illustrated in the first column in Fig.3, the length of the baseline will not have a significant impact on the performance of results.

\begin{figure}
\vspace{-0cm}
\centering
\hspace{-1.5ex}
\subfigure[Original RS image]{
\begin{minipage}{2.7cm}
\centering
\includegraphics[width=2.45cm]{rs1.pdf}
\end{minipage}%
}%
\subfigure[Corrected image]{
\begin{minipage}{2.7cm}
\centering
\includegraphics[width=2.5cm]{rs_rec.pdf}
\end{minipage}%
}%
\subfigure[RS depth truth]{
\begin{minipage}{2.7cm}
\centering
\includegraphics[width=2.5cm]{1_rs_depth.pdf}
\end{minipage}}

\noindent\subfigure[Original RS+GS]{
\begin{minipage}{2.7cm}
\centering
\includegraphics[width=2.5cm]{RS+GS1.pdf}
\end{minipage}}%
\subfigure[Corrected+GS]{
\begin{minipage}{2.7cm}
\centering
\includegraphics[width=2.5cm]{rec_RS+gs1.pdf}
\end{minipage}}
\subfigure[Depth error]{
\begin{minipage}{2.7cm}
\centering
\includegraphics[width=3.1cm]{depth_err_r1.pdf}
\end{minipage}}
\centering
\vspace{-0.3cm}
\caption{Qualitative experiment on synthetic RS images. The experiments exhibit the results of RS image correction and the depth error estimated by SGM algorithm.}
\vspace{-0.5cm}
\end{figure}
\subsection{Experiments on synthetic RS images}
\label{ssec:subhead}
To demonstrate the qualitative results of our relative pose solver, we generate RS images by the simulator and 3D models as in \cite{Zhuang_Rolling_shutter_aware_2017}. We first set the ground truth of all parameters such as $\mathbf{d}$, $\mathbf{w}$ and $\varphi$, then we can calculate the pose of each row on each frame in Fig.2. According to the poses, the simulator will generate GS images, then synthesize RS images from these GS images by extracting the corresponding row pixels of the corresponding GS image. In experiments, we set the image resolution as $900\times900$ with a 1384.6 focal length, the magnitude of $\mathbf{d}$ is 0.3, and the magnitude of $\mathbf{w}$ is fixed to $0.4\degree$. Besides, both translation and rotation include motion on x, y, and z-axis. Then we use SIFT \cite{Lowe_SIFT_2004} to extract feature points and RANSAC to reduce the effect of noise on pose estimation. Our experimental result is reported in Fig.(4).
Fig.4(a) shows the original RS image. And in Fig.4(c).(f), we can observe the depth truth of RS image and the depth error between depth truth and depth estimated by the SGM algorithm. We notice that the left part of the depth error map is set to 0 because the SGM algorithm cannot estimate the depth value of the left area on the left image.
Fig.4(b) shows the corrected image using our method and the depth estimated by the SGM algorithm.
In Fig.4(d), we overlay the original RS image and the ground truth GS image rendered by the pose of the first row, the red area shows the distortion. And then we overlap the corrected image and the GS image as shown in Fig.4(e). Comparing Fig.4(d) and Fig.4(e), we can find that the rectified RS image almost coincides with the GS image, except for the surrounding area where the SGM algorithm cannot estimate the depth and a small area where the depth error is more than 2m.
The experiments prove that our proposed algorithm can also estimate the motion well, and the SGM algorithm can be used to estimate the depth of stereo RS images and get an accurate depth value which can be used in RS image correction.

\vspace{-0.5cm}
\section{CONCLUSION}
\label{sec:typestyle}
In this paper, we have proposed a novel linear relative pose solver for stereo rolling shutter cameras under consecutive frames using 9 pairs of correspondences for the left and right camera respectively. We showed that for RS stereo images, traditional SGM algorithm can be used to estimate the depth map and achieve an accurate result. Experimental results on synthetic RS images show that we can well correct the RS image using our proposed relative pose solver and depth estimated by the SGM algorithm.

\noindent\textbf{Acknowledgements:} This research was supported in part by the National Natural Science Foundation of China under Grants 61871325, 61420106007, and 61671387 and the National Key Research and Development Program of China under Grant 2018AAA0102803.

\bibliographystyle{IEEE