\section{Introduction}

\quad The Fisher Information Matrix (FIM) summarizes the amount of information about data related to parameters of interest, bearing in mind that the central role of data analysis is to extract information from data. We can see that the FIM plays a critical role in the theory and practice of statistical modeling, including system identification \cite{b1}, parameter estimation \cite{b2}, information theory \cite{b3}, and other areas. Let $\pmb{Z}$ be the collection of $n$ random vectors $[\pmb{z}_1,\ldots,\pmb{z}_n]$. In the following discussion, assume that $p(\pmb{Z}|\pmb{\theta})$ is the density function for the model of interest, where $\pmb{\theta}$ represents the unknown $p\times 1$ parameter vector, and $L(\pmb{\theta}|\pmb{Z})$ is the log-likelihood function of $p(\pmb{Z}|\pmb{\theta})$. Then the $p\times p$ FIM is defined as:
$$\pmb{F}(\pmb{\theta})=\mathbb{E}\bigg[\frac{\partial L}{\partial \pmb{\theta}}\frac{\partial L}{\partial \pmb{\theta}^T}\bigg].$$

In general, it is difficult to compute the expectation of a product of nonlinear score functions with multiple parameters. On the other hand, if the second partial derivatives of $L$ exists, i.e., the Hessian matrix $\pmb{H}(\pmb{Z}|\pmb{\theta})$, then under specific regularity conditions \cite{b4}, the FIM can be represented as:
$$\pmb{F}(\pmb{\theta})=-\mathbb{E}\bigg[\frac{\partial^2 L}{\partial \pmb{\theta}\partial \pmb{\theta}^T}\bigg]=-\mathbb{E}[\pmb{H}(\pmb{\theta})],$$
which is generally easier to calculate compared with the basic definitional form.

In that there are no closed forms for $\pmb{F}(\pmb{\theta})$ in many practical problems and given the importance of the FIM, one way is to use averages of the estimated Hessian matrix based on pseudo data as the Monte Carlo estimation of the FIM. Moreover, some research studies present methods to estimate the FIM in general or specific cases (e.g., \cite{b5}, \cite{b6}). In particular, we can estimate the FIM via measurements of a Hessian matrix based on the simultaneous perturbation approach. Spall in \cite{b5} discusses a basic resampling-based approach. Later several papers present improved approaches, including a feedback-based approach in \cite{b7} and a method where information on some of the FIM is available in \cite{b8}.  In this paper, we introduce another enhanced Monte Carlo method: the independent perturbation approach. Ref. \cite{b7} also mentions the independent perturbation approach, but we do a more thorough analysis here.

In Section 2, we demonstrate theories of the basic resampling-based method and the independent perturbation method. Simultaneously, we show that the independent perturbation method reduces the variance of entries of the estimated FIM by from $O(1/N)$ to $O(1/(nN))$, where $n$ is the sample size of the data and $N$ is a measure of the Monte Carlo averaging. Section 3 illustrates the efficiency of the enhanced method through two numerical examples, respectively. The signal-plus-noise example verifies the theoretical ratio about variance reduction and compares the performance at each diagonal entries of the estimated FIM for the two methods. The mixture Gaussian example utilizes the relative error based on the spectral norm to show the obvious improvement of the independent perturbation method. Section 4 gives a conclusion and ideas about future work.

\section{Theoretical Analysis}
\subsection{Standard Method}

\quad We summarize aspects of the basic method in \cite{b5} below. Assume that $\pmb{Z}_\mathrm{pseudo}(i)$ is a set of $n$ independent random vectors generated by the Monte Carlo simulation that follow the distribution with parameters $\pmb{\theta}$, and $\hat{\pmb{H}}_{k|i}$ is the $k$th estimate of
$\pmb{H}(\pmb{\theta})$ at $\pmb{Z}_\mathrm{pseudo}(i)$. Additionally, suppose that $\pmb{\Delta}_{k|i}=\left[\Delta_{k1|i},\Delta_{k2|i},\ldots,\Delta_{kp|i}\right]^T$ is a random vector with a zero mean where the scalar elements are independent and identical symmetrically distributed random variables that are uniformly bounded and satisfy $\mathbb{E}\left(\left|1/\Delta_{kj|i}\right|\right)<\infty$, where $j=1,\ldots,p$, and $\pmb{\Delta}_{k|i}^{-1}\equiv[\Delta_{k1|i}^{-1},\Delta_{k2|i}^{-1},\ldots,\Delta_{kp|i}^{-1}]^T$.

Here we generate an estimate of $\pmb{H}(\pmb{\theta})$ by simultaneous perturbation (SPSA) methods:
\begin{equation}
	\hat{\pmb{H}}_{k|i}=\frac{1}{2}\left[\frac{\delta \pmb{g}_{k|i}}{2c}(\pmb{\Delta}_{k|i}^{-1})^T+\left(\frac{\delta \pmb{g}_{k|i}}{2c}(\pmb{\Delta}_{k|i}^{-1})^T\right)^T\right],
\end{equation}
where $\delta \pmb{g}_{k|i}=\pmb{g}(\pmb{\theta}+c\pmb{\Delta}_{k|i}|\pmb{Z}_\mathrm{pseudo}(i))-\pmb{g}(\pmb{\theta}-c\pmb{\Delta}_{k|i}|\pmb{Z}_\mathrm{pseudo}(i))$, and $\pmb{g}(\cdot)$ is the (estimated if there is not enough information) gradient function of $L$, and $c>0$ is a small constant.

The Monte Carlo-based estimate of $\pmb{F}(\pmb{\theta})$ in \cite{b5}, denoted $\bar{\pmb{F}}_{M,N}(\pmb{\theta})$, is
\begin{equation}
	\bar{\pmb{F}}_{M,N}(\pmb{\theta})\equiv-\frac{1}{N}\sum_{i=1}^N\frac{1}{M}\sum_{k=1}^M\hat{\pmb{H}}_{k|i}.
\end{equation}
The first “inner” average calculates Hessian estimates at a given $\pmb{Z}_\mathrm{pseudo}(i)$ ($i = 1, 2, \ldots, N$) based on $M$ values of $\hat{\pmb{H}}_{k|i}$, and the second “outer” average sums these sample means of Hessian estimates across $N$ values of pseudo data.

To get the most intuitive comparison, we set $M=N=1$ without loss of generality since (2) indicates that Hessian estimates do not depend on $N$. Additionally, while $\hat{\pmb{H}}_{k|i} (k=1,\ldots,M)$ are estimated based on the same pseudo data with fixed $i$, \cite{b10} shows $M=1$ is the optimal choice if the $n$ vectors entering each $\pmb{Z}_\mathrm{pseudo}(i)$ are mutually independent. Then the variance of the $jj$th entry of the estimate $\bar{\pmb{F}}_{M,N}(\pmb{\theta})$ is
\begin{equation}
	\mathrm{var}\left\{\left[\bar{\pmb{F}}_{M,N}\right]_{jj}\right\}=\mathrm{var}\left(\hat{H}_{jj}\right),
\end{equation}
where $\hat{H}_{jj}$ denotes the $jj$th entry of $\hat{\pmb{H}}=\hat{\pmb{H}}_{1|1}$. We can use the calculation at $M=N=1$ to readily extend to arbitrary matrix. Let $O_{(\cdot)}(c^2)$ denote a random “big-$O$” term, where the subscript represents the relevant randomness. For example, $O_{\pmb{Z},\pmb{\Delta}_1}(c^2)$ denotes a random “big-$O$” term depends on $\pmb{Z}_\mathrm{pseudo}(1)$ and $\pmb{\Delta}_1$ such that $O_{\pmb{Z},\pmb{\Delta}_1}(c^2)/c^2$ is bounded almost surely as $c\rightarrow0$. Then, by \cite{b8}, the $jj$th entry of $\hat{\pmb{H}}$ is
\begin{equation}
	\hat{{H}}_{jj}={H}_{jj}+\sum_{l\neq j}{H}_{jl}\frac{\Delta_{1l}}{\Delta_{1j}}+O_{\pmb{Z},\pmb{\Delta}_1}(c^2).
\end{equation}
Assume that
\begin{equation}
	\mathbb{E}(\Delta_{1l}/\Delta_{1j})=0,\; \mathrm{var}(\Delta_{1l}/\Delta_{1j})=v;
\end{equation}
where $v$ is a constant, then, given the independence of the $\left\{\Delta_{1j}\right\}$, and the fact that $\mathbb{E}({H}_{jl})=-{F}_{jl}$ ($j,l=1,2,\ldots,p$), we have:
\begin{align*}
	\mathrm{var}(\hat{{H}}_{jj})_{\text{basic}}=&\mathbb{E}\big([{H}_{jj}+\sum_{l\neq j}{H}_{jl}\frac{\Delta_{1l}}{\Delta_{1j}}+O_{\pmb{Z},\pmb{\Delta}_1}(c^2)]^2\big)-\big(\mathbb{E}[{H}_{jj}+\sum_{l\neq j}{H}_{jl}\frac{\Delta_{1l}}{\Delta_{1j}}+O_{\pmb{Z},\pmb{\Delta}_1}(c^2)]\big)^2\\
	=&\mathbb{E}\big({H}_{jj}^2+2{H}_{jj}\sum_{l\neq j}{H}_{jl}\frac{\Delta_{1l}}{\Delta_{1j}}+[\sum_{l\neq j}{H}_{jl}\frac{\Delta_{1l}}{\Delta_{1j}}]^2+O_{\pmb{Z},\pmb{\Delta}_1}(c^2)\big)-\big(\mathbb{E}\big[{H}_{jj}\big]+O(c^2)\big)^2\\
	=&\mathbb{E}({H}_{jj}^2)+\sum_{l\neq j}\mathbb{E}\bigg[{H}_{jl}\frac{\Delta_{1l}}{\Delta_{1j}}\bigg]^2-\big(\mathbb{E}\big[{H}_{jj}\big]\big)^2+O(c^2)\\
	=&\mathbb{E}({H}_{jj}^2)+v\sum_{l\neq j}\mathbb{E}({H}_{jl}^2)-{F}_{jj}^2+O(c^2)\\
	=&\mathrm{var}({H}_{jj})+v\sum_{l\neq j}\left[\mathrm{var}({H}_{jl})+{F}_{jl}^2\right]+O(c^2)\\
	=&\mathrm{var}({H}_{jj})+v\sum_{l\neq j}\mathrm{var}({H}_{jl})+v\sum_{l\neq j}{F}_{jl}^2+O(c^2).
\end{align*}

\subsection{Implementation with Independent Perturbations}

\quad The estimation of $\pmb{F}(\pmb{\theta})$ can be enhanced if the $n$ vectors entering each $\pmb{Z}_\mathrm{pseudo}(i)$ are mutually independent. Note that the $n$ vectors do not have to be identically distributed --- just independent. By the independence, we have
$$\pmb{F}(\pmb{\theta})=\sum_{t=1}^n\pmb{F}^{(t)},$$
where $\pmb{F}^{(t)}$ denotes the FIM for each $t$.

In particular, the variance of the entries of each Hessian estimate $\hat{\pmb{H}}_{k|i}$ in the standard method can be reduced by decomposing
$\hat{\pmb{H}}_{k|i}$ into a sum of $n$ independent estimates $\hat{\pmb{H}}_{k|i}^{(t)}$, each corresponding to one vector in the pseudo data. To distinguish the way that $\pmb{\Delta}_{k|i}$ is generated in the basic method and the independent perturbation method, $\pmb{\Delta}_I$ is used here to represent the perturbation term $\pmb{\Delta}_{1|1}$ in the independent perturbation method, where $\pmb{\Delta}_I=\sum_{t=1}^n\pmb{\Delta}_I^{(t)}$. Then a separate perturbation vector $\pmb{\Delta}_I^{(t)}$ can be applied to each independent estimate, which results in variance reduction in the estimate $\bar{\pmb{F}}_{M,N}$. Hence, based on assumption (5), similar to the above calculation process, we have:
\begin{align*}
	&\mathrm{var}(\hat{{H}}_{jj})_{\text{indep}}=\mathrm{var}(\sum_{t=1}^n\hat{{H}}_{jj}^{(t)})=\sum_{t=1}^n\mathrm{var}(\hat{{H}}_{jj}^{(t)})\\
	=&\sum_{t=1}^n\{\mathbb{E}([{H}_{jj}^{(t)}+\sum_{l\neq j}{H}_{jl}^{(t)}\frac{\Delta_{Il}^{(t)}}{\Delta_{Ij}^{(t)}}+O_{\pmb{Z},\pmb{\Delta}_I}(c^2)]^2)-(\mathbb{E}[{H}_{jj}^{(t)}+\sum_{l\neq j}{H}_{jl}^{(t)}\frac{\Delta_{Il}^{(t)}}{\Delta_{Ij}^{(t)}}+O_{\pmb{Z},\pmb{\Delta}_I}(c^2)])^2\}\\
	=&\sum_{t=1}^n\{\mathbb{E}[({H}_{jj}^{(t)})^2]+v\sum_{l\neq j}\mathbb{E}[({H}_{jl}^{(t)})^2]-({F}_{jj}^{(t)})^2+O(c^2)\}\\
	=&\sum_{t=1}^n\{\mathrm{var}\left({H}_{jj}^{(t)}\right)+v\sum_{l\neq j}[\mathrm{var}({H}_{jl}^{(t)})+({F}_{jl}^{(t)})^2]+O(c^2)\}\\
	=&\sum_{t=1}^n\{\mathrm{var}({H}_{jj}^{(t)})+v\sum_{l\neq j}\mathrm{var}({H}_{jl}^{(t)})\}+v\sum_{t=1}^n[\sum_{l\neq j}({F}_{jl}^{(t)})^2]+O(c^2).
\end{align*}
The variance of the $j$th diagonal entry of $\hat{\pmb{H}}$ generated from the standard method can be written as:
\begin{align*}
	\mathrm{var}(\hat{{H}}_{jj})_{\text{basic}}=&\mathrm{var}(\sum_{t=1}^n{H}_{jj}^{(t)})+v\sum_{l\neq j}\mathrm{var}(\sum_{t=1}^n{H}_{jl}^{(t)})+v\sum_{l\neq j}[\sum_{t=1}^n{F}_{jl}^{(t)}]^2+O(c^2)\\
	=&\sum_{t=1}^n\{\mathrm{var}\left({H}_{jj}^{(t)}\right)+v\sum_{l\neq j}\mathrm{var}({H}_{jl}^{(t)})\}+v\sum_{l\neq j}[\sum_{t=1}^n{F}_{jl}^{(t)}]^2+O(c^2).
\end{align*}
Therefore, the difference of the variance of the diagonal entries in the estimate of $\pmb{F}(\pmb{\theta})$ between the standard method and the independent perturbation method is
$$\mathrm{var}(\hat{{H}}_{jj})_{\text{basic}}-\mathrm{var}(\hat{{H}}_{jj})_{\text{indep}}=2v\sum_{l\neq j}\sum_{t_1<t_2}{F}_{jl}^{(t_1)}{F}_{jl}^{(t_2)}+O(c^2).$$

\subsection{Theoretical Comparison via Variance Reduction}

\quad Denoting $A=\sum_{t=1}^n\{\mathrm{var}({H}_{jj}^{(t)})+v\sum_{l\neq j}\mathrm{var}({H}_{jl}^{(t)})\}$, the result above shows that the variance difference between these two method is independent of $A$. That is, $\mathrm{var}(\hat{{H}}_{jj})_{\text{indep}}$ and $\mathrm{var}(\hat{{H}}_{jj})_{\text{basic}}$ contain the same components that are made up of the variance of entries of the true Hessian matrix. Hence, we should mainly consider the difference between $v\sum_{t=1}^n[\sum_{l\neq j}({F}_{jl}^{(t)})^2]$ and $v\sum_{l\neq j}[\sum_{t=1}^n{F}_{jl}^{(t)}]^2$. However, the result above cannot directly reflect the degree of variance reduction. Let us consider the ratio $\mathrm{var}(\hat{{H}}_{jj})_{\text{indep}}$ to $ \mathrm{var}(\hat{{H}}_{jj})_{\text{basic}}$, which can be written as:
\begin{equation}
	\frac{\text{var}(\hat{{H}}_{jj})_{\text{indep}}}{\text{var}(\hat{{H}}_{jj})_{\text{basic}}}=\frac{A+v\sum_{t=1}^n\left[\sum_{l\neq j}({F}_{jl}^{(t)})^2\right]+O(c^2)}{A+v\sum_{l\neq j}[\sum_{t=1}^n{F}_{jl}^{(t)}]^2+O(c^2)}.
\end{equation}
Usually $c$ is small (for example, 0.0001), so let us focus on the second part: off-diagonal entries of the true Fisher matrix. With a fixed $j$, let us for convenience denote entries ${F}_{jl}^{(t)}=f_l^{(t)}$, where $l\neq j$. Under some assumptions, the ratio  of the second part  in (6), $v\sum_{t=1}^n[\sum_{l\neq j}({F}_{jl}^{(t)})^2]/v\sum_{l\neq j}[\sum_{t=1}^n{F}_{jl}^{(t)}]^2$, is $O(1/n)$. That is, assume that:
\begin{enumerate}
	\item There are at least 2 non-zero off-diagonal entries in each true Fisher matrix $\pmb{F}^{(t)}$.
	\item For every $j=1,2,\ldots,p$, the non-zero sequence $\{f_l^{(t)}\}_{t=1}^n$ is limited within $[f_{l}^{(\text{min})},f_{l}^{(\text{max})}]$ while fixing $l$, where $f_{l}^{(\text{min})}$ is the lower bound of non-zero entries in $\{|f_l^{(t)}|\}_{t=1}^n$, and $f_{l}^{(\text{max})}$ is the upper bound of $\{|f_l^{(t)}|\}_{t=1}^n$. For convenience, denote $r_{1_{j}}=\sum_{l\neq j}^p(f_{l}^{(\text{max})}/{f_{l}^{(\text{min})}})^2$.
	\item Denote $f_l^{(+)}$ as the sum of all positive elements in the sequence, and $f_l^{(-)}$ as the sum of all negative elements in the sequence, the ratio is $r_{2_{j}}=(f_l^{(+)}+|f_l^{(-)}|)^2/{(f_l^{(+)}-|f_l^{(-)}|)^2}$ such that $r_{1_{j}}r_{2_{j}}=O(1)$ as $n$ goes to infinity.
\end{enumerate}

Based on the above conditions, we have:

\textbf{Proposition 1:} Suppose that conditions (1), (2), and (3) hold for all $j$. Then the independent perturbations above reduce the variance of the diagonal entries in the estimate of $\pmb{F}(\pmb{\theta})$ by $O(1/n)+O(c^2)$.

The details about the proof have been shown in \cite{b10}, so we get:
\begin{align*}
	\frac{\text{var}(\hat{{H}}_{jj})_{\text{indep}}}{\text{var}(\hat{{H}}_{jj})_{\text{basic}}}
	\leq&\frac{A+\frac{\text{C}}{n}B_{\text{basic}}+O(c^2)}{A+B_{\text{basic}}+O(c^2)},
\end{align*}
where $C=r_{1_{j}}r_{2_{j}}$ is a constant, $A=\sum_{t=1}^n\{\mathrm{var}({H}_{jj}^{(t)})+v\sum_{l\neq j}\mathrm{var}({H}_{jl}^{(t)})\}=O(n)$, and $B_{\text{basic}}=v\sum_{l\neq j}^p[\sum_{t=1}^nf_{l}^{(t)}]^2\leq n^2v\sum_{l\neq j}^p[f_{l}^{(\text{max})}]^2=O(n^2)$. Therefore,
\begin{equation}
	\frac{\text{var}\left(\hat{{H}}_{jj}\right)_{\text{indep}}}{\text{var}\left(\hat{{H}}_{jj}\right)_{\text{basic}}}=O\bigg(\frac{1}{n}\bigg)+O(c^2),
\end{equation}
where $O(c^2)$ is often so small that it can be ignored (see Example 1 in Section 3).

In addition, similar to Section 4.2 in \cite{b5}, \cite{b10} shows that in the independent case, $M>1$ is considered only when it is costly to generate new pseudodata vectors. This is despite the fact that we mentioned earlier that the average of $M$ estimated Hessian matrices, which is based on the same pseudodata vector $\pmb{Z}_{\text{pseudo}}$, is usually used to estimate a Fisher information matrix.

That is, in general, if it is allowed to compute $C=MN$ Hessian matrix estimates, we can maximize the accuracy of $\bar{\pmb{F}}_{M,N}$ when each estimate $\hat{\pmb{H}}_{k|i}$ is generated based on a new $\pmb{Z}_{\text{pseudo}}$ (i.e., $M=1$).

\textbf{Proposition 2:} Assume that elements $\pmb{\Delta}_{k|i}^{(t)}$ and $\pmb{Z}^{(t)}_{\text{pseudo}}(i)$ are mutually independent, where $k=1,\ldots,M$, $i=1,\ldots,N$, and $t=1,\ldots,n$. For a fixed $C=MN$, the variance of each entry in $\bar{\pmb{F}}_{M,N}(\pmb{\theta})$ is minimized when $M=1$.

\section{Numerical Examples}

\quad Here we introduce two examples, which are also used in \cite{b11}. The first one is under the signal-plus noise setting. We calculate the variance ratio to verify the theoretical result and analyze the distinction among the variance ratio for different diagonal entries. The second one compares the performance for the two methods when the pseudo data are mixture Gaussian distributed. At the same time, we show the improvement of the enhanced method from time analysis in both examples.

Both of the mixture distribution and the signal-plus-noise settings arise regularly in the literature. Mixture problems are thoroughly reviewed in \cite{b12} and \cite{b13}, with applications in information theory in \cite{b14} and bivariate quantile estimation in \cite{b15}. The signal-plus-noise problem with non-identical noise distributions arises in practical problems where measurements are collected with varying quality of information across the sample. Some practical implications are discussed in \cite{b16} relative to the initial conditions in a Kalman filter model and \cite{b17} in the context of outlier analysis.

\subsection{Multivariate Normal with Signal-Plus-Noise}

\quad Assume that the $\pmb{z}_i$ are independently and normally distributed with mean $\pmb{\mu}$ and covariance $\pmb{\Sigma}+\pmb{P}_i$ for all $i$, where $\pmb{\mu}$ and $\pmb{\Sigma}$ are unknown parameters to be estimated. It can be interpreted that $\pmb{z}_i$ is the observed values of the $N(\pmb{\mu},\pmb{\Sigma})$ distributed signal with independent $N(\pmb{0},\pmb{P}_i)$ distributed noise. This setting has been widely discussed; for example, \cite{b16} considers estimating the initial mean vector and covariance matrix of a state-space model and \cite{b18} demonstrates parameter estimation of random-coefficient linear models.

Let us consider a 3-dimensional case: $\text{dim}(\pmb{z}_i)=3$, $n=30$, $\pmb{\mu}=\pmb{0}$, $\pmb{\Sigma}$ is a matrix with all 0.5's except 2's on the diagonal (which is a little different from that in \cite{b10}), and $\pmb{P}_i=\sqrt{i}\pmb{U}^T\pmb{U}$, where $\pmb{U}$ is a $3\times3$ matrix with entries generated from uniform (0,1) distribution. $\pmb{\theta}$ is the collection of unique elements in $\pmb{\mu}$ and $\pmb{\Sigma}$, so the FIM is $p\times p$ where $p=9$. Reference \cite{b16} gives the gradient function of log-likelihood and the analytical form of the FIM in this setting.

In particular, assume that components of the perturbation vector $\pmb{\Delta}_{k|i}$ follow Bernoulli $\pm1$ distribution. In both methods, we estimate the Hessian matrix by using the gradient of the log-likelihood function under $M=N=1$. The ratio values for estimated variance of the diagonal entries of the estimated FIM from the standard method to that from the independent perturbation method based on 2 million replicates are shown in Table 1. Furthermore, we change $n$ from 30 to 200 to observe the performance.

\begin{table}[htbp]
	\caption{Variance of $\bar{\pmb{F}}_{M,N}(\pmb{\theta})_{jj}$}
	\begin{center}
		\begin{tabular}{|c|c|c|c|}
			\hline
			\textbf{Diagonal}& \textbf{Independent} & \textbf{Standard} & \textbf{Ratio}\\
			\hline
			$j=1$ &8.24  &60.12 & 0.14\\
			\hline
			$j=2$ &13.32 &56.58 &0.23\\
			\hline
			$j=3$ & 4.88 & 24.48& 0.20\\
			\hline
			$j=4$ & 5.00 & 	11.18 &0.45\\
			\hline
			$j=5$ &23.49 & 39.13 & 0.60\\
			\hline
			$j=6$ &14.85 &26.69 & 0.56\\
			\hline
			$j=7$ &8.69 & 13.84 & 0.63\\
			\hline
			$j=8$ &12.91 & 19.70 & 0.65\\
			\hline
			$j=9$ &3.04 &4.98 & 0.61\\
			\hline
		\end{tabular}
		\label{tab1}
	\end{center}
\end{table}

Table 1 indicates that there is an obvious reduction of the variance of entries in the estimated FIM over the enhanced method compared with that of the standard method for every $j$, which matches our theoretical conclusion in Section 2.

\begin{table}[htbp]
	\caption{Time cost/$seconds$}
	\begin{center}
		\begin{tabular}{|c|c|c|c|}
			\hline
			&\textbf{$n=30$}& \textbf{$n=100$}&\textbf{$n=200$}\\
			\hline
			Independent Method &8726&	28914&	57066			\\
			\hline
			Standard Method &5259&	16978&	33541\\
			\hline
			Ratio & 0.6027&	0.5872	&0.5878\\
			\hline
		\end{tabular}
		\label{tab1}
	\end{center}
\end{table}

On the other hand, researchers usually need to trade off accuracy and time cost in experiments. Results with high accuracy generally take more time to obtain, but the payoff may not be worthwhile if the computational cost is too high. Therefore, under the same setting as above, Figures 1 and 2, and Table 2 show the performance of the standard method and the enhanced method as $n$ is equal to 30, 100, and 200 based on the variance ratio $\text{var}(\hat{H}_{jj})_{\text{indep}}/\text{var}(\hat{H}_{jj})_{\text{basic}}$ and time cost. To better compare the variance ratio with $O(1/n)$ as $n$ grows, we add the reference lines $3.5/n$ and $13/n$ in Figures 1 and 2, separately. Figure 1 and 2 reveal that the change of variance ratio as $n$ grows is similar to the trend of $O(1/n)$. That is, the curves follow the overall slope given by the reference line. But the time cost ratio in Table 2 is stable over the range of $n$ considered.

\begin{figure}[H]
	\centering
	\includegraphics[width=9.5cm,height=7.3cm]{mu1}
	\caption{Variance ratio for $j=1,2,3$}
\end{figure}

\begin{figure}[H]
	\centering
	\includegraphics[width=9.5cm,height=7.3cm]{sigma1}
	\caption{Variance ratio for $j=4,\ldots,9$}
\end{figure}

 In addition, under specific regularity conditions \cite{b4}, the $jj$th entry of $\pmb{F}(\pmb{\theta})$ is the expectation of $\partial^2 L/\partial {\theta}_j^2$, where ${\theta}_j(j=1,\ldots,9)$ corresponds to three elements in $\pmb{\mu}$ and six unique entries in $\pmb{\Sigma}$ respectively. We can see that the estimates of diagonal entries of $\pmb{F}(\pmb{\theta})$ corresponding to $\pmb{\mu}$ part are more stable than those corresponding to the $\pmb{\Sigma}$ part. Furthermore, when $n=30$, the variance ratio values in Figure 1 are significantly lower than 0.6027, which is the time cost ratio in Table 2. Part of the variance ratio values in Figure 2, however, are higher than 0.6027.

Note that compared with \cite{b10}, we only change the diagonal values of the true covariance matrix $\pmb{\Sigma}$. But in \cite{b10}, all variance ratio values are much lower than 0.9586, which is the time cost ratio as $n=30$. Here, the performance of variance ratios that corresponds to the $\pmb{\mu}$ part is much better than that corresponding to the $\pmb{\Sigma}$ part. That is, the variance ratio values corresponding to the $\pmb{\mu}$ part are lower and closer to $1/n$. This might be because the true variance value for each dimension becomes large, which tend to expand the range of generated values and affects the estimation for the covariance part. Then the fluctuation of corresponding diagonal entries (for the covariance part) in the estimated FIM are relatively obvious compared with that for the mean part.

Although Table 2 indicates that it may not be worth the time (compare the time cost ratio with part of the variance ratio values) to improve accuracy when $n=30$, the results for $n=100, 200$ from the independent method achieve a significant reduction on variance with acceptable and reasonable time cost, which means that our enhanced method is applicable in practice.

\subsection{Mixture Gaussian Distribution}

\quad Denote the generated vector $\pmb{Z}=[z_1,z_2,\ldots,z_n]^T$, and assume that $z_i \;(i=1,\ldots,n)$ are independently and identically distributed with probability density function:
$$f(z,\pmb{\theta})=\lambda\exp(-(z-\mu_1)^2/(2\sigma_1^2))/\sqrt{2\pi\sigma_1^2}+(1-\lambda)\exp(-(z-\mu_2)^2/(2\sigma_2^2))/\sqrt{2\pi\sigma_2^2},$$
where $\pmb{\theta}=[\lambda,\mu_1,\sigma_1^2,\mu_2,\sigma_2^2]^T$. There are $5(5+1)/2=15$ unique entries in the FIM. Furthermore, we use the average of the negative Hessian matrix based on the sample vector over $10^6$ independent replicates to approximate (since \cite{b19} showed that the Hessian matrix is attainable) the true FIM, which is not computable in this example.

Here consider the same setting as that in \cite{b11}, where $\pmb{\theta}=[0.2,0,4,1,9]^T$, $M=2$, $N=40000$, $c=0.0001$, and elements of the perturbation $\pmb{\Delta}_{k|i}$ follow Bernoulli $\pm1$ distribution for all $k\;(k=1,\ldots,N)$ and $i\;(i=1,\ldots,M)$. We estimate the true FIM $\pmb{F}(\pmb{\theta})$ given the gradient of the log-likelihood function, and measure the performance of methods by computing the sample mean of $\|\pmb{F}_{\text{est}}(\pmb{\theta})-\pmb{F}(\pmb{\theta})\|/\|\pmb{F}(\pmb{\theta})\|$ based on 50 independent replicates, where $\|\bullet\|$ represents the spectral norm.
Table 3 gives the typical estimated FIM based on the standard method and the independent method, respectively. By typical, we mean it is the estimated FIM corresponding to the 25th value of $\|\pmb{F}_{\text{est}}(\pmb{\theta})-\pmb{F}(\pmb{\theta})\|/\|\pmb{F}(\pmb{\theta})\|$ in descending order given 50 independent replicates. Note that the true FIM is approximated by sample averages over $10^6$ independent replicates. Table 4 summarizes the performances (accuracy and time cost) between the standard method and the independent approach.

\begin{table}[htbp]
	\caption{Simulation Results for Example 3.2}
	\begin{center}
		\begin{tabular}{|c|c|}
			\hline
			True FIM&$\left[
			\begin{array}{ccccc}796.8750&	0&	0&	0&	0\\
				0&	7.5000&	1.5000&	0&	0\\
				0&	1.5000&	3.1876&	0&	0\\
				0&	0&	0&	3.3333&	-0.0741\\
				0&	0&	0&	-0.0741&	0.1523\end{array}
			\right]
			$\\
			\hline
			Typical $\bar{\pmb{F}}(\pmb{\theta})_{\text{basic}}$&$\left[
			\begin{array}{ccccc}
				796.8752	&1.5414&	1.8820	&0.3499&	-0.7470\\
				1.5414&	7.5004&	1.4958&	0.01550&	-0.0075\\
				1.8820&	1.4958&	3.1868&	0.0157&	-0.02251\\
				0.3499	&0.0155	&0.0157&	3.3332&	-0.0863\\
				-0.7470&	-0.0075&	-0.02251&	-0.0863&	0.1521
			\end{array}
			\right]
			$\\
			\hline
			Typical $\bar{\pmb{F}}(\pmb{\theta})_{\text{indep}}$&$\left[
			\begin{array}{ccccc}
				796.8752&0.2536&	0.0278&	0.4273&	-0.0930\\
				0.2536	&7.4949&	1.4922&-0.0043&	-0.0009\\
				0.0278&1.4922&	3.1768	&-0.0019&	-0.0026\\
				0.4273	&-0.0043&	-0.0019&3.3323	&-0.0754\\
				-0.0930&	-0.0009&	-0.0026&	-0.0754&	0.1509
			\end{array}
			\right]$
			\\
			\hline
		\end{tabular}
		\label{tab1}
	\end{center}
\end{table}

\begin{table}[htbp]
	\caption{Comparison between Two Methods}
	\begin{center}
		\begin{tabular}{|c|c|c|c|}
			\hline
			&Independent &Standard&Ratio\\
			\hline
			Accuracy ($\|\pmb{F}_{\text{est}}(\pmb{\theta})-\pmb{F}(\pmb{\theta})\|/\|\pmb{F}(\pmb{\theta})\|$) &0.00063&	0.00330&	0.19\\
			\hline
			Time Cost/$seconds$&1760&	1212&	1.45\\
			\hline
		\end{tabular}
		\label{tab1}
	\end{center}
\end{table}
 Table 4 shows the obvious advantage of the enhanced method with independent perturbations while taking time cost into account. It is worth noting that our accuracy of the basic approach is higher than that shown in \cite{b11} because for each time, i.e., fixing $k$ and $i$, \cite{b11} only generates one sample following the mixture Gaussian distribution, but we generate $n$ samples, which might result in a more accurate estimation.

\section{Conclusion and Future Work}

\quad It is often hard to estimate the FIM in real-world models because of high dimension or nonlinear functions. This paper presents an enhanced Monte Carlo method with independent perturbations to estimate the Fisher information matrix in problems where it is not analytically available. The theoretical analysis and numerical examples show that this approach improves the estimation performance by reducing variance of elements compared to the basic Monte Carlo approach.

Moreover, the above results indicate that the performance might change even if we only change certain parameters in the model; for example, consider the time cost ratios in Example 1 and \cite{b10}. However, as $n$ increases, the advantage of the enhanced method is apparent with no doubt.

In future work, in addition to considering the performance of estimated FIM in other practical problems (for example, Bayesian statistics \cite{b20} and the EM algorithm \cite{b8}), it is important to account for the relative benefit and cost of reducing variance, especially when $n$ is very large, where we cannot ignore the cost of generating the additional perturbation vectors in the independent perturbation method.

\begin{singlespace}


\section{Introduction}
\IEEEPARstart{V}{isual} trackers based on Siamese networks \cite{li2019siamrpn++, wang2019fast, zhang2019deeper, siamban}  have achieved state-of-the-art results on various tracking benchmark datasets \cite{wu2015object, fan2019lasot, mueller2016benchmark, kiani2017need, muller2018trackingnet}. A key component of state-of-the-art trackers \cite{li2019siamrpn++, zhu2018distractor, li2018high} is based on anchor boxes, which are boxes of different sizes and aspect ratios that select as tracking candidates. Anchor boxes are widely used in visual tracking \cite{li2019siamrpn++, zhang2019deeper, zhu2018distractor, li2018high}, and can achieve a good trade-off in speed-accuracy. These methods \cite{li2019siamrpn++, zhang2019deeper, zhu2018distractor, li2018high} first generate anchor boxes densely over the image,
they then find the candidate regions of the target by learning a strong classifier, and they finally refine their coordinates through a regression network.

\begin{figure}[!t]
\centering
\setlength{\abovecaptionskip}{0.cm}
\includegraphics[width=8.8 cm, scale=.5]{heatmap.pdf}
\caption{Visualization results of feature maps (keep only the highest score for clarity) of the four boundaries (green box) predicted by Siamese network. The feature maps of the four boundaries are fed to a modified corner pooling layer to predict the corner heatmaps and offsets of the green box. Specifically, we combine the offset prediction and penalty strategies to make the green box close to the ground truth box (red box).}
\label{heatmap}
\end{figure}
\vspace{-0.1cm}

However, trackers based on anchor boxes have two drawbacks. First, these trackers need very large-scale anchor boxes, e.g. more than 2k in SiamRPN \cite{li2018high} and more than 3k in SiamRPN++ \cite{li2019siamrpn++}. It is because that the anchor-based tracker is trained to classify whether each anchor box contains a target region, and a large number of anchor boxes are required to ensure that they have sufficient overlapping regions with the ground truth boxes. In training, the positive and negative samples are defined by calculating the IoU between anchor boxes and ground truth boxes \cite{tian2019fcos}. They \cite{li2019siamrpn++, zhang2019deeper, zhu2018distractor, li2018high} simply think that the anchor boxes are positive samples if the IoU is larger than a hand-crafted threshold and are negative samples if the IoU is less than a hand-crafted threshold. However, recent research \cite{zhang2020bridging} indicates that this fixed hand-craft threshold will affect the performance of anchor-based algorithms. As a result, a tiny part of the anchor boxes can ensure overlap with ground truth boxes. In this design, this will produce a huge imbalance between positive samples and negative samples \cite{li2019siamrpn++} and slows down training \cite{law2018cornernet, lin2017focal}. Furthermore, the imbalance of positive and negative samples will have a huge impact on a model's performance \cite{zhang2020bridging, paa-eccv2020}. For example, the training for the classification model is inefficient as most locations are easy negative samples, which will lead to a degraded classification model \cite{lin2017focal}. If the classifier cannot produce a correct candidate region for the target location, then the regression network will not be accurate \cite{2017Faster, Ocean_2020}.

Second, trackers with hand-crafted anchor boxes will introduce many design options. This includes how to select the number, size, and aspect ratio of anchor boxes. This design is mostly inspired by ad-hoc heuristics and it can become more complex when a network predicts multiple-level feature maps, each of which uses different features and has specific anchor boxes \cite{li2019siamrpn++,law2018cornernet, lin2017focal,liu2016ssd}. These parameters in anchor-based trackers need to be tuned carefully. Even with careful design, it is difficult for trackers to deal with large changes of target shapes, especially for small targets, due to the fixed scales and aspect ratios of the anchor boxes.

In contrast, the location and boundary of the target that humans can identify are not obtained by predefined candidate boxes \cite{siamban}. Therefore, we can design a tracking algorithm without any predefined anchor boxes. Motivated by these analyses, we introduce a SiamCorners method for visual tracking, which makes our tracking algorithm more flexible and general by discarding anchor boxes. Our proposed method consists of four steps. First, we crop four boundary images containing target boundary information from an initial frame as the template images of the SiamCorners.  Specifically, the four boundary images regard the top, left, bottom and right boundary of the bounding box of the initial frame as the initial locations, respectively, and it then crops at a specific ratio of the width or height of the box. The SiamCorners then predicts the feature maps (see figure \ref{heatmap}) of the four boundaries of the bounding box based on the input of template images and search images. Second, inspired by the recently proposed CornerNet \cite{law2018cornernet, law2019cornernet}, which is used to predict a bounding box as a pair of corners in object detection. However, the original CornerNet is class-specific, and hence cannot be used directly for visual tracking. Consequently, we introduce a modified corner pooling layer to help the network integrate the target-specific feature in the Siamese framework. The entire corner pooling layer is used to predict the multiple top-left corners and bottom-right corners of the bounding box by the predicted feature maps of the four boundaries. Third, a decoding process is performed to predict the corner locations based on the heatmaps and offset maps. It is worth noting that we chose a Non-maximum-suppression (NMS) operation to remove some abnormal corners. Finally, multi-level features fusion and some penalty strategies are adopted to select the best tracking bounding box. Intuitively, features from earlier convolution outputs of deep networks will contain rich low-level feature information such as shape, color, which is important for tracking localization. Features from the latter convolution outputs of the deep networks will contain rich semantic information that is beneficial to deal with some challenging factors like huge deformation and background clutter. In our design, multi-layer features are fed into the multiple SiamCorners modules individually. We then use a new penalty function to select the best tracking box from the candidate corners generated by the layer-wise aggregation strategy. Specifically, we first employ the location and scale penalty function that considers the changes of width, height, and center of the tracking target compared to the previous frame. This penalty strategy is used to re-rank the corners' score to remove the outliers. We then utilize a penalty term to suppress large changes in the ratio and size of the corners. Subsequently, we select the corners with the largest score as a candidate tracking box. Finally, a linear interpolation is employed to make the tracking box changing smoothly. This design simplifies the output of the whole network and eliminates the complex design of the anchor boxes.

In summary, the proposed SiamCorners mainly includes the following contributions:

	\begin{itemize}
	    \item We propose an anchor-free proposal network (SiamCorners) that avoids many hand-crafted configurations during designing the anchor boxes and makes our method more flexible and general.
        \item We present a modified corner pooling layer to integrate target-specific information into the corner prediction.
	    \item We formulate a multi-layer features fusion strategy to help SiamCorners benefit from both low-level information and high-level information.
	    \item We introduce a novel penalty term that considers the changes of the center, width, and height of the target in the adjacent frames, which makes SiamCorners more suitable for tracking tasks.

	\end{itemize}

Based on these contributions, we demonstrate the effectiveness and efficiency of SiamCorners on five tracking benchmark datasets: OTB100 \cite{wu2015object}, UAV123 \cite{mueller2016benchmark}, LaSOT \cite{fan2019lasot}, NFS30 \cite{kiani2017need} and TrackingNet \cite{muller2018trackingnet}.  Our SiamCorners approach achieves state-of-the-art results on the five datasets. Moreover, the wide range of ablation experiments are performed on the UAV123 and the LaSOT datasets.

\section{RELATED WORKS}
Generic object tracking has made great progress in recent years following the development of several new methods \cite{li2019siamrpn++, wang2019fast, zhang2019deeper, siamban, ruan2018multi, liu2017robust, Yuankai2019Learning, liu2020multi-task}. In particular, methods based on Siamese networks play an important role in visual tracking.

\textbf{Siamese Networks:} Recently, trackers based on Siamese networks have received much attention because of their well-balanced in accuracy-efficiency. SiameseFC \cite{bertinetto2016fully} calculates the similarity prediction of template images and search images through a full-convolutional network, while running over 100FPS. CFNet \cite{valmadre2017end} introduces correlation filter layers into the Siamese tracking framework, which can achieve end-to-end representation learning of tracking targets. However, its performance is not significantly improved compared to SiameseFC. SiamRPN \cite{li2018high} develops the tracking task as a one-shot detection task by introducing region proposal subnetwork in the Siamese network, which is trained offline in an end-to-end manner on large-scale image pairs. The region proposal network contains two branches, one for classification branch and  another for regression branch. The classification branch predicts the target location based on the fixed anchor boxes, while the regression branch fine-tunes the predicted target coordinates. Therefore, the design of anchor boxes has a great impact on the final tracking results. DaSiamRPN \cite{zhu2018distractor} further improves the discriminative ability of the model by introducing an effective sampling strategy to control the imbalanced distribution of the training data. However, note that these trackers \cite{bertinetto2016fully, valmadre2017end, li2018high, zhu2018distractor} built their algorithm in shallow
networks similar to AlexNet \cite{krizhevsky2012imagenet}, the components they proposed did not benefit from low-level and high-level feature information due to relatively shallow networks. While our method uses multiple SiamCorners modules to predict multiple candidate corners according to
the layer-wise aggregation features of the deep networks. Similarly, SiamRPN++ \cite{li2019siamrpn++} enables SiamRPN tracker to take advantage of deep features by introducing deeper networks, such as ResNet-50 \cite{he2016deep} or deeper networks. Unlike anchor-based trackers \cite{li2019siamrpn++, zhu2018distractor, li2018high}, our SiamCorners tracker is more flexible and general since it eliminates the design of anchor boxes.
\begin{figure*}[!t]
\centering
\centerline{{\includegraphics[width=7.0in]{systemnetwork.pdf}}}
\caption{The overview of our SiamCorners architecture, which includes the Siamese feature extractor followed by the top-left corner and bottom-right corner branches in parallel. $\bigstar$ denotes depth-wise correlation operation. $z_{\theta}(\theta\in \{t,l,b,r\})$ denote the four boundary images cropped from the original image.
}
\label{sys}
\end{figure*}

\textbf{Anchor-free Mechanism:} Anchor-free trackers recently became popular in visual tracking tasks due to their relatively simple and superior performance. FCAF \cite{han2019fully} and Ocean \cite{Ocean_2020} respectively introduce a simple anchor-free tracker by using a variant of the FCOS \cite{tian2019fcos} detector.
They \cite{Ocean_2020, han2019fully} first generate some hand-crafted points instead of a large number of anchor boxes as training samples, and then use a classification network to find some points that may be the locations of a tracking target. Finally, they utilize a regression network to predict a 4D vector, which represents the distance from the target location to the four sides of the bounding box. Unlike these anchor-free trackers, our SiamCorners directly predicts the target box as a pair of corners, which does not include a classification network to determine the coarse location in some hand-crafted points. Therefore, SiamCorners maintains a high-efficiency running speed by eliminating the classification phase of the anchor-free trackers \cite{Ocean_2020, han2019fully}.

Recently, corner-based detectors have attracted widespread attention in the object detection community due to their relatively new paradigms. CornerNet \cite{law2018cornernet} was the first to propose a corner-based detector that converts the bounding box of the target into a pair of corner predictions. CornerNet uses an hourglass network \cite{newell2016stacked} to obtain the heatmaps including the top-left corners and the bottom-right corners of the bounding boxes. It then uses a series of feature embeddings to match these corners. To accurately predict the target, CornerNet introduces a corner pooling to localize the corner. CornerNet-Lite \cite{law2019cornernet} utilizes CornerNet as a baseline to introduce an attention mechanism so that the detector does not require detailed processing for each pixel in the input image, thereby improving the inference efficiency of the entire detector. CenterNet \cite{duan2019centernet:} regards each detection target as a triplet, which predicts the centers of the bounding boxes along with corners. MartixNets \cite{rashwan2020matrixnets} proposes an aspect ratio and scale aware architecture so that CornerNet can eliminate corner pooling operation. Instead of matching corners by feature embeddings in CornerNet, CentripetalNet \cite{CentripetalNet} matches the predicted corners by whether the centripetal shift results of the corners are aligned, where the centripetal shift denotes the spatial shift from the corner locations to the center of the bounding box.

There are two key differences between tracking and detection tasks. First, anchor-free detectors usually aim to detect different categories of targets in an input image, while in tracking, we aim to match specific targets in search images. Therefore, we follow a common Siamese architecture \cite{li2019siamrpn++} that can discriminate the target from the search images. Second, the detector based on anchor-free usually assumes that the target category is predefined. The common dataset MS COCO \cite{lin2014microsoft} of object detection only contains 91 target categories. Thus, independent corner pooling layers are trained in \cite{law2018cornernet, law2019cornernet, duan2019centernet:, rashwan2020matrixnets, CentripetalNet} for each object class. In contrast, the target class is unknown in tracking \cite{li2019siamrpn++, he2017robust, he2016connected}. Furthermore, unlike the object detection tasks, the target does not belong to any pre-defined classes or is represented by existing training datasets. Instead, target-specific corner pooling layers are required for test frames incorporating target information in Siamese framework.

\section{Proposed Method}
\subsection{Overview}

\begin{figure}[!t]
\setlength{\abovecaptionskip}{0.cm}
\setlength{\belowcaptionskip}{-0.cm}
\centering
\makeatletter\def\@captype{figure}\makeatother
{
   \label{crop}
   \includegraphics[width=8.1 cm]{Visio_croped.pdf}}

\caption{Examples of four boundary images fed to the template branch. The boundary images are cropped from the $t_{w,h}$ width or height of the red box of the input image.}
\label{crop}
\end{figure}

In this work, we propose an anchor-free method, SiamCorners, which detects a target as a pair of corners. A common Siamese network is used to integrate target-specific information into the tracking architecture.

Fig. \ref{sys} provides an overview of SiamCorners. As illustrated in Fig. \ref{sys}, the template branch has four inputs: i) the top-boundary image $z_{t}$, ii) the left-boundary image $z_{l}$, iii) the bottom-boundary image $z_{b}$, iv) the right-boundary image $z_{r}$ (see Fig. \ref{crop}). The target images of the top (left) boundary and the bottom (right) boundary respectively take the four boundaries of the target's bounding box as the initial location, and are cropped from $t_{w,h}$ (we set $t_{w,h}=0.5$ in all experiments) times the width (height) of the bounding box of the original image, respectively. Furthermore, the other  regions of the boundary images are constant padding. Then, we use the same data augmentation strategy as SiamRPN++ \cite{li2019siamrpn++} to resize the four boundary images to a fixed size. Afterwards, the four boundary images of the bounding box are fed to ResNet-50 \cite{he2016deep} network, respectively. The feature maps from different layers, namely $conv3$, $conv4$, and $conv5$ blocks of ResNet-50, are fed to the two branches. Specifically, the feature maps of the top and the left boundary images are fed to the top-left corner prediction branch, and the feature maps of the bottom and the right boundary images are fed to the bottom-right corner prediction branch. The features of different layers in ResNet-50 are followed by a convolution network. Depth-wise correlation operations are further performed to obtain correlation features. The correlation features are then fed to the corner pooling layers to get the corner heatmaps and offsets, respectively (see Section \ref{pooling}). After obtaining the predicted heatmaps and offsets, a simple decoding process is used to obtain the corner candidates (see Section \ref{decoding1}). Finally, the multi-layer feature fusion and the post-processing operation are used to get the final bounding box (see Section \ref{fusion}).

\subsection{Tracking Corners Loss}
We predict two sets of corner heatmaps that are used to calculate the top-left corners and the bottom-right corners of the target, respectively. Let $\tilde{x}\in\mathbb{R}^{\emph{H}\times\emph{W}}$ be a search image, where $\emph{H}$ and $\emph{W}$ are the height and width of the image, respectively. The corner heatmaps of the convolutional network are denoted as $\hat{X}\in \left [0,1\right]^{\frac{H}{s}\times\frac{W}{s}}$, where $\emph{s}$ is a scale factor mapped from the search image to the corner heatmaps.

We train the corner network following Law and Deng \cite{law2019cornernet}. For each heatmap, there is only one ground-truth label for a corner as a positive location and the other locations are negative, which will make the network difficult to learn. For efficient training, we reduce the penalty when the negative locations are still within the radius of the positive location. We determine the radius by the at least \emph{d} IoU (we set \emph{d} = 0.5 in all experiments) between a pair of corners and the ground-truth box. For each ground truth corner $c(c_{x},c_{y})\in\mathbb{R}^{2}$ in the search image, we compute the location of corner heatmap ${\tilde{c}(\tilde{c}_{x},\tilde{c}_{y})=\left \lfloor \frac{c(c_{x},c_{y})}{s} \right \rfloor}$. Given the radius and location ${\tilde{c}}$, the amount of penalty reduction is calculated by 2D Gaussian, $X_{xy}=exp\left(-\frac{(x-\tilde{c}_{x})^2+(y-\tilde{c}_{y})^2}{2\delta_{\tilde{c}}^2} \right)$, where $\delta_{\tilde{c}}$ is 1/3 of the radius. Again, $X_{xy}$ denotes the ``ground truth" heatmap designed with unnormalized Gaussians. Let $\hat{X}_{xy}$ be the prediction score of network at the location $(x,y)$. The training objective is a variant of focal loss \cite{lin2017focal}:
\begin{equation}
\mathcal{L}_{tr} =\frac{-1}{K} \sum_{xy}
\begin{cases}
(1-\hat{X}_{xy})^\alpha log(\hat{X}_{xy})& \mbox{if}X_{xy}=1 \\ \label{1}
(1-X_{xy})^\beta(\hat{X}_{xy})^\alpha log(1-\hat{X}_{xy}) & \mbox{otherwise}
\end{cases}
\end{equation}
where $\alpha$ and $\beta$ denote the weighting factors which control the contribution of each corner (we set $\alpha=2$ and $\beta=4$ in all experiments). $K$ is the number of candidate corners. In eq.(\ref{1}), the $(1-X_{xy})$ term reduces the penalty when the predicted corners are within the radius of the positive location.

However, the output size of the network is usually smaller than the original image as the network becomes deeper. There will be a discretization error when mapping a location $(m,n)$ of the original image to the location $({\left\lfloor\frac{m}{s}\right\rfloor,\left\lfloor\frac{n}{s} \right \rfloor})$ in the heatmaps. Moreover, some precision will be lost in tracking when remapping locations from heatmaps to the original image. To address the discretization error caused by network stride, we introduce an offset network that is used to slightly adjust the corner offset before remapping the heatmap to the original image. The corner offset on the heatmap is calculated as:

\begin{equation}
\emph{\textbf{o$_k$}}=\left(\frac{m_{k}}{s}-\left \lfloor \frac{m_{k}}{s}\right \rfloor,\frac{n_{k}}{s}-\left \lfloor \frac{n_{k}}{s}\right \rfloor \right)
\label{o}
\end{equation}
where $\emph{\textbf{o$_k$}}$ is the corner offset, $m_{k}$ and $n_{k}$ are the $k-th$ corner coordinates of the object in the original image. We employ the smooth L1 Loss \cite{girshick2015fast} at corner locations:
\begin{equation}
\mathcal{L}_{off} =\frac{1}{K}\sum_{k=1}^KSmoothL1Loss(\emph{\textbf{o$_k$}},\hat{\emph{\textbf{o}}}_{k})
\label{off}
\end{equation}
where $\emph{\textbf{o}}_{k}$ and $\hat{\emph{\textbf{o}}}_{k}$ are ``ground truth" offset and offset network prediction, respectively. Therefore, the whole method optimizes the loss function $\mathcal{L}$, which is expressed as:

\begin{equation}
\mathcal{L} = \mathcal{L}_{tr}+\lambda\mathcal{L}_{off} \label{L}
\end{equation}
where $\lambda$ (we set $\lambda=1$ in all experiments)  is a hyper-parameter that controls the balance of the two loss terms.

\subsection{Corner Pooling}
\label{pooling}

\begin{figure*}[!t]
\centering
\centerline{{\includegraphics[width=6.5 in]{cornerpoolingModule.pdf}}}
\caption{The corner pooling module starts with the two feature maps, which are modified from the residual block. The bottom pooling and right pooling use the max-pools methods of eq.\ref{5} and \ref{6}, respectively. The entire corner pooling layer has two branches that predict the target's heatmaps and offsets, respectively. In the figure, the $3\times3$ Conv-BN-Relu module consists of a $3\times3$ convolution layer, a Batch Normalization layer, and a ReLU layer. Similarly, the $3\times3$ Conv-BN-ReLU-Conv module consists of the Conv-BN-ReLU module and a $1\times1$ convolution layer.}
\label{poolingmodule}
\end{figure*}

In this section, we describe how to use a modified corner pooling to predict corner heatmaps and offsets in Siamese architecture. The existing corner pooling is originally proposed for object detection, we make modifications to adapt it to object tracking. Given the template images $z_{\theta},\theta \in \{t,l,b,r\}$ and the search image $\tilde{x}$, the feature maps are efficiently denoted by parameterizing ${f^{a}(z_{\theta})={\phi(z_{\theta})}}, a\in \{3,4,5\}$ and ${f^{a}(\tilde{x})={\phi(\tilde{x})}}$, where $\phi$ is a convolutional neural network (CNN). Here, $(f^{3}(z_{\theta}),f^{3}(\tilde{x}))$, $(f^{4}(z_{\theta}),f^{4}(\tilde{x}))$, and $(f^{5}(z_{\theta}),f^{5}(\tilde{x}))$ denote the output features of \emph{conv3}, \emph{conv4}, and \emph{conv5} branches in Fig. \ref{sys}, respectively. Afterward, following the Siamese architecture, we aim to predict the depth-wise correlation feature: ${{f_{\theta}^{a}}={f^{a}(z_{\theta})\bigstar f^{a}(\tilde{x})}}$. Therefore, the specific-target features are integrated into the SiamCorners network.

Here, we will describe in detail how to determine whether a pixel $(i,j)$ is a bottom-right corner in correlation maps. On a $\frac{H}{s}\times\frac{W}{s}$ correlation map, the corner pooling layer max-pools the all feature on $f_{b}^{a}$ to get a feature vector $b_{ij}^{a}$ from $(i,0)$ to $(i,\frac{W}{s})$ , and max-pools the feature map $f_{r}^{a}$ to get a feature vector $r_{ij}^{a}$ from $(0,j)$ to $(\frac{H}{s},j)$. This calculation stage can be expressed as:

\begin{equation}
{b_{ij}^{a}} =
\begin{cases}

max(f_{b_{ij}}^{a}, b_{i(j-1)}^{a}) &\mbox   {if } 0<j\le\frac{W}{s} \\ \label{5}
f_{b_{i0}}^{a}  & \mbox{otherwise}
\end{cases}
\end{equation}

\begin{equation}
{r_{ij}^{a}} =
\begin{cases}
max(f_{r_{ij}}^{a}, r_{(i-1)j}^{a}) &\mbox   {if } 0<i\le\frac{H}{s} \\ \label{6}
f_{r_{0j}}^{a}  & \mbox{otherwise}
\end{cases}
\end{equation}

We then add ${b_{ij}^{a}}$ and ${r_{ij}^{a}}$ to get the final corner pooling result. Similarly, the pooling results ${t_{ij}^{a}}$ and ${l_{ij}^{a}}$ are obtained by max-pooling the feature vectors horizontally from right to left and vertically from bottom to top on the correlation feature maps $f_{t_{ij}}^{a}$ and $f_{l_{ij}}^{a}$, respectively. Finally, the result of top-left corner pooling is obtained by the sum of the two pooling results.

The architecture of the whole corner pooling is shown in Fig.\ref{poolingmodule}.  The starting part of the corner pooling module corresponding to the two inputs ($f_{b}^{a}$ and $f_{r}^{a}$) is a modified variant of the residual block \cite{he2016deep}, respectively. After the design of the residual block, the corner pooling feature is followed by a $3\times3$ Conv-BN layer and its output is added to the two projection shortcuts. Finally, it predicts the corner heatmaps and offsets by feeding the modified residual feature through a $3\times3$ Conv-BN-ReLU  layer followed by a $3\times3$ Conv-BN-ReLU-Conv layer, respectively. It is worth noting that the modified corner pooling and the original corner pooling in CornerNet \cite{law2018cornernet} have different designs in the network structure and pooling mechanism. As shown in Fig. \ref{poolingmodule}, the modified corner pooling introduces bottom pooling and right pooling in the right and bottom boundary feature maps ($f_{r}^{a}$ and $f_{b}^{a}$) respectively. Both the bottom and right boundary features are followed by projection shortcuts to enhance the corner prediction accuracy. In contrast, the original corner pooling \cite{law2018cornernet}  for object detection task predicts the bottom-right corner by using the bottom-right corner pooling and a single projection shortcut in a unified feature map.

\subsection{Decoding for Corner Heatmaps and Offsets}
\label{decoding1}
In this section, we aim to predict the tracking bounding box by the given corner heatmaps and offsets. Following Section \ref{pooling}, the final outputs of corner pooling have four parts: i) the top-left corner heatmaps. ii) the bottom-right corner heatmaps. iii) the top-left corner offsets. iv) the bottom-right corner offsets. After corner pooling, the top-left and bottom-right corner heatmaps are first followed by a sigmoid function, respectively. Then, a NMS operation is adopted to remove redundant corner locations. We further obtain the corresponding top-left corner location $
(X{_{T_{l}}^{a}},Y{_{T_{l}}^{a}}) \in\mathbb{R}^{N\times 2}$ and bottom-right corner location $(X{_{B_{r}}^{a}},Y{_{B_{r}}^{a}}) \in\mathbb{R}^{N\times 2}$ by selecting the top $N$ (we set $N=15$ in all experiments) corner scores in top-left and bottom-right corner heatmaps, respctively. To overcome the impact of network stride, we denote $(\Delta X_{T_{l}}^{a}, \Delta Y_{T_{l}}^{a})$ and $(\Delta X_{B_{r}}^{a}, \Delta Y_{B_{r}}^{a})$ as the corresponding offsets of the top-left corners and the bottom-right corners, respectively. Afterward, the refined top $N$ corner sets $COR^{*}={(\hat X{_{T_{l}^{a}},\hat Y{_{T_{l}^{a}}},\hat X{_{B_{r}^{a}}},\hat Y{_{B_{r}^{a}}},S_{cor}^{a})}}$ can be predicted by the following formula:

\begin{equation}
\begin{split}
\hat X{_{T_{l}^{a}}}
=X{_{T_{l}}^{a}}+\Delta X_{T_{l}}^{a},\quad \hat Y{_{T_{l}^{a}}}
=Y{_{T_{l}}^{a}}+\Delta Y_{T_{l}}^{a}\\
\hat X{_{B_{r}^{a}}}
=X{_{B_{r}}^{a}}+\Delta X_{B_{r}}^{a},\quad \hat Y{_{B_{r}^{a}}}
=Y{_{B_{r}}^{a}}+\Delta Y_{B_{r}}^{a}\\
S_{cor}^{a}
=\frac{1}{2}(T_{l_{sco}}^{a}+B_{r_{sco}}^{a})
\end{split}
\end{equation}
where $S_{cor}^{a}$ is a pair of corner scores which are the average score of the top-left corner score $T_{l_{sco}}^{a}$ and bottom-right corner score $B_{r_{sco}}^{a}$. Here, $COR^{*}\in\mathbb{R}^{N\times 5}$. After the top $N$ corners $COR^{*}$ are obtained, we use multi-layer features fusion and some penalty strategies to select the optimal tracking corners.

\subsection{Multi-level Features Fusion and Corners Selection}
\label{fusion}
\textbf{Features Fusion:} Intuitively, visual tracking requires rich feature representations that include features from low-level to high-level, and resolutions from fine to coarse and scales from small to large \cite{li2019siamrpn++}. To this end, we introduce a multi-layer features fusion strategy to help SiamCorners benefit from low-level features and high-level features. Following Section \ref{decoding1}, we concatenate $\hat X{_{T_{l}^{a}}}$, $\hat Y{_{T_{l}^{a}}}$, $\hat X{_{B_{r}^{a}}}$ , $\hat Y{_{B_{r}^{a}}}$ and $S_{cor}^{a}$ and use a simple mapping to get the corner locations $X_{T_{l}}$, $Y_{T_{l}}$,  $X_{B_{r}}$, $Y_{B_{r}}$ and $S_{cor}$ corresponding to the original image, respectively. Therefore, the corner sets in the original image can be denoted as a $3N\times5$ dimension term $COR={(X{_{T_{l}},Y{_{T_{l}}},X{_{B_{r}}},Y{_{B_{r}}},S_{cor})}}$.

\textbf{Corners Selection:} After obtaining the corners $COR$, we introduce a post-process strategy to make the SiamCorners more suitable for tracking tasks. In visual tracking, the target in nearby frames almost does not have large variations in size and shape. We calculate the overall scales $c_{sc}$ of the corners by the following formula:

\begin{equation}
c_{sc}^2
=(w+p)(h+p) \label{14}
\end{equation}
where $w=X{_{B_{r}}}-X{_{T_{l}}}$ and $h=Y{_{B_{r}}}-Y{_{T_{l}}}$ denote the width and height of the target, and $p=\frac{1}{2}(w+h)$ denotes the padding. Similarly, the same principle is used to calculate the overall scales $\hat c_{sc}^{2}$ of corners at the last frame. We use a penalty term to suppress large changes in the size and ratio of the target:
\begin{equation}
C_{penalty}
=e^{\eta*max(\frac{c_{r}}{\hat c_{r}},\frac{\hat c_{r}}{c_{r}})*max(\frac{c_{sc}}{\hat c_{sc}},\frac{\hat c_{sc}}{c_{sc}})} \label{15}
\end{equation}
where $c_{r}$ and $\hat c_{r}$  denote the corners' ratio of height and width in the current frame and the last frame, respectively. Here, $\eta$ is a hyper-parameter. Therefore, the score $S_{pen}$ of each corner in the $COR$ is calculated by the weighted penalty as: $S_{pen}^u=C_{penalty}^u*S_{cor}^u$, $u=1,2,...,3N$. Unlike SiamRPN \cite{li2019siamrpn++,li2018high}, we further introduce a novel index term $U(u), u=1,2,...,3N$, which represents the index of the variation degree of each corner in $COR$ from large to small compared with the last frame. The metric of variation degree is determined by the sum of the absolute values of the difference between the center, width, and height of the target between two adjacent frames. For example, when $u=1$, $S_{pen}^{U(1)}$ means that the $U(1)-th$ corner in the corner sets $COR$ has the largest motion change relative to the last frame. Finally, the scores $S_{final}^{U(u)}$ of each corner can be expressed as:

\begin{equation}
S_{final}^{U(u)}
=S_{pen}^{U(u)}*(1-\gamma)+S_{han}^{u}*\gamma \label{16}
\end{equation}
where $\gamma$ is a hyper-parameter that controls the balance between the two score terms. $S_{han}^{u}$ is a variant of the Hanning window since we only take the monotonically increasing part of the Hanning window. After obtaining the final score $S_{final}$ of each corner, the expression of the corner set $COR$ is replaced by $COR_{final}={(X{_{T_{l}},Y{_{T_{l}}},X{_{B_{r}}},Y{_{B_{r}}},S_{final})}}$. Afterwards, we pick the corners with the highest score in the $COR_{final}$ as the final tracking bounding box. Furthermore, a linear interpolation operation is used to update the target size so that the shape can be changed smoothly.

Finally, we summarize the SiamCorners in algorithm {\ref{tracking_all}}.

\begin{algorithm}[!t]
  \caption{Tracking as one-shot}
  \label{tracking_all}
  \begin{algorithmic}[1]
    \Require
       initial template images: $z_{\theta}$;
       search image: $\tilde{x}$;
    \Ensure
     Tracking bounding box $x^{*}$
    \State ${f^{a}(z_{\theta})\gets{\phi(z_{\theta})}}$; \quad {\footnotesize \#  Feed the template images into CNN.}
    \Repeat
      \State ${f^{a}(\tilde{x})\gets{\phi(\tilde{x})}}$;
      \State ${{f_t^{a}},{f_l^{a},f_b^{a}},{f_r^{a}}\gets{f^{a}(z_{\theta})\bigstar f^{a}(\tilde{x})}}$ {\footnotesize \# Depth-wise correlation}
      \State $t^{a},l^{a},b^{a},r^{a}\gets{CornerPooling({{f_t^{a}},{f_l^{a},f_b^{a}},{f_r^{a}}})}$ {\footnotesize \# Section \ref{pooling}}
      \State $COR^{*}\gets {Decoding(t^{a},l^{a},b^{a},r^{a})}$  {\footnotesize \# Decoding process (Section \ref{decoding1})}
       \For{$a=3,4,5$}
            \State $COR\gets{{(\hat X{_{T_{l}^{a}},\hat Y{_{T_{l}^{a}}},\hat X{_{B_{r}^{a}}},\hat Y{_{B_{r}^{a}}},S_{cor}^{a})}}}$ {\footnotesize \# Concatenation}
        \EndFor
        \State $COR_{final}\gets{COR}$\quad {\footnotesize \# Update corner scores using eq.(\ref{14})-(\ref{16})}
        \State $x^{*}\gets{PostProcessing(COR_{final}})$ {\footnotesize {\# Selecting a pair of corners and using a linear interpolation operation} }

    \Until {end of sequence}
  \end{algorithmic}
\end{algorithm}

\section{Experimental Results}
Our method is implemented in Python using PyTorch. On a single NVIDIA RTX 2080 GPU, SiamCorners achieves a tracking speed of 42 FPS by employing ResNet-50 as backbone. To facilitate further research, the completed training and testing code will be released at https://github.com/yangkai12/SiamCorners.
\subsection{Training Details}
 We use ResNet-50 pre-trained on ImageNet \cite{russakovsky2015imagenet} as the backbone network. The training splits of the GOT-10k \cite{huang2019got}, LaSOT \cite{fan2019lasot}, COCO \cite{lin2014microsoft}, YouTube-BoundingBoxes \cite{real2017youtube}, ImageNet VID \cite{russakovsky2015imagenet}, and ImageNet DET datasets are used to train our network to learn how to measure the similarity of input image pairs in visual tracking. In both training and testing phases, we set the sizes of the template image and test image to 127 pixels and 255 pixels, respectively. To reduce overfitting, we implemented some data augmentation strategies in the input image pairs, such as random horizontal flipping, random color jittering, random scaling, random cropping, adjusting the brightness and contrast of an image. We use stochastic gradient descent (SGD) with a batchsize of 58 to optimize the overall training loss from eq. (\ref{L}). We set a learning rate of 0.001 for earlier 5 epochs to train the top-left branch and bottom-right branch. For the latter 15 epochs, the whole network is trained in an end-to-end manner, where the learning rate decays exponentially from 0.005 to 0.0005, while weight decay and momentum are set to 0.0001 and 0.9, respectively. Unlike SiamRPN++ \cite{li2019siamrpn++}, we do not have a finetune backbone operation in the last 10 epochs, which will increase training time exponentially.

\subsection{Testing Metrics}

\begin{figure}[!t]
\centering
\makeatletter\def\@captype{figure}\makeatother
\subfigure{
   \label{fig_Success_Plot}
   \includegraphics[width=4.1 cm]{OTB100_success.pdf}}
\hspace{0.005 in}
\subfigure{
   \label{fig_Precision_Plot}
   \includegraphics[width=4.1 cm]{OTB100_precision.pdf}}
\caption{Experimental results of the methods on the OTB100 dataset.}
\label{OTB}
\end{figure}

In this experiment, we evaluate our method on tracking benchmarks OTB100 \cite{wu2015object}, UAV123 \cite{mueller2016benchmark}, LaSOT \cite{fan2019lasot}, NFS30 \cite{kiani2017need}, and TrackingNet \cite{muller2018trackingnet}. We use success scores and precision scores to evaluate the performance of the trackers on OTB100, UAV123 and NFS30. The success scores are defined as the percentage of success frames when the overlap threshold between the predicted result and the ground truth varies from 0 to 1, where a success frame means that the tracking result is larger than a given threshold. For a fair comparison, all trackers are ranked by the area under curve (AUC) of the success
plot. The precision score is computed as the percentage of frames, which is defined that the Euclidean distance between the center of the prediction target and the ground truth is within a given distance threshold. Here, the different distance thresholds are used to evaluate all trackers. For the LaSOT and TrackingNet datasets, we use the success score, precision score, and normalized precision score to evaluate the performance of all trackers, where the definitions of success and precision score are the same as previously described. Meanwhile, the normalized precision is defined as the normalization of precision over the size of the ground truth bounding box. The trackers are then ranked by AUC score with normalized precision between 0 and 0.5. Generally, the AUC score of success plot is mainly used to rank all trackers.

\subsection{State-of-the-art Comparison}

\begin{figure}[!t]
\centering
\makeatletter\def\@captype{figure}\makeatother
\subfigure{
   \label{fig_Success_Plot}
   \includegraphics[width=4.1 cm]{UAV123_success.pdf}}
\hspace{0.005 in}
\subfigure{
   \label{fig_Precision_Plot}
   \includegraphics[width=4.1 cm]{UAV123_precision.pdf}}
\caption{Experimental results of the methods on the UAV123 dataset.}
\label{UAV}
\end{figure}

\begin{figure*}[!t]
\centering
\subfigure[scale variation]{\includegraphics[width=0.245\textwidth]{1.pdf}}
\subfigure[aspect ratio change]{\includegraphics[width=0.245\textwidth]{2.pdf}}
\subfigure[low resolution]{\includegraphics[width=0.245\textwidth]{3.pdf}}
\subfigure[fast motion]{\includegraphics[width=0.245\textwidth]{4.pdf}}
\subfigure[full occlusion]{\includegraphics[width=0.245\textwidth]{5.pdf}}
\subfigure[partial occlusion]{\includegraphics[width=0.245\textwidth]{6.pdf}}
\subfigure[out-of-view]{\includegraphics[width=0.245\textwidth]{7.pdf}}
\subfigure[background clutter]{\includegraphics[width=0.245\textwidth]{8.pdf}}
\subfigure[illumination variation]{\includegraphics[width=0.245\textwidth]{9.pdf}}
\subfigure[viewpoint change]{\includegraphics[width=0.245\textwidth]{10.pdf}}
\subfigure[camera motion]{\includegraphics[width=0.245\textwidth]{11.pdf}}
\subfigure[similar object]{\includegraphics[width=0.245\textwidth]{12.pdf}}
\caption{The success plots on UAV123 dataset for 12 challenging factors: scale variation,  aspect ratio change, low resolution, fast motion, full occlusion, partial occlusion, out-of-view, background clutter, illumination variation, viewpoint change, camera motion, and similar object. The legend shows the AUC score for each method.}
\label{uav123}
\end{figure*}

\begin{figure*}[!t]
\centering
\renewcommand\arraystretch{1.1}
\subfigure{\includegraphics[width=0.164\textwidth,height=1.635cm]{bicycle-9-79.jpg}}\hspace{-0.027in}
\subfigure{\includegraphics[width=0.164\textwidth,height=1.635cm]{bicycle-9-459.jpg}}\hspace{-0.027in}
\subfigure{\includegraphics[width=0.164\textwidth,height=1.635cm]{bicycle-9-1065.jpg}}\hspace{-0.027in}
\subfigure{\includegraphics[width=0.164\textwidth,height=1.635cm]{boat-4-491.jpg}}\hspace{-0.027in}
\subfigure{\includegraphics[width=0.164\textwidth,height=1.635cm]{boat-4-1070.jpg}}\hspace{-0.027in}
\subfigure{\includegraphics[width=0.164\textwidth,height=1.635cm]{boat-4-1899.jpg}}\hspace{-0.027in}\vspace{-0.117in}
\subfigure{\includegraphics[width=0.164\textwidth,height=1.635cm]{boat-17-711.jpg}}\hspace{-0.027in}
\subfigure{\includegraphics[width=0.164\textwidth,height=1.635cm]{boat-17-929.jpg}}\hspace{-0.027in}
\subfigure{\includegraphics[width=0.164\textwidth,height=1.635cm]{boat-17-1474.jpg}}\hspace{-0.027in}
\subfigure{\includegraphics[width=0.164\textwidth,height=1.635cm]{cat-3-133.jpg}}\hspace{-0.027in}
\subfigure{\includegraphics[width=0.164\textwidth,height=1.635cm]{cat-3-568.jpg}}\hspace{-0.027in}
\subfigure{\includegraphics[width=0.164\textwidth,height=1.635cm]{cat-3-1910.jpg}}\hspace{-0.027in}\vspace{-0.117in}
\subfigure{\includegraphics[width=0.164\textwidth,height=1.635cm]{crocodile-10-331.jpg}}\hspace{-0.027in}
\subfigure{\includegraphics[width=0.164\textwidth,height=1.635cm]{crocodile-10-1213.jpg}}\hspace{-0.027in}
\subfigure{\includegraphics[width=0.164\textwidth,height=1.635cm]{crocodile-10-1409.jpg}}\hspace{-0.027in}
\subfigure{\includegraphics[width=0.164\textwidth,height=1.635cm]{flag-5-361.jpg}}\hspace{-0.027in}
\subfigure{\includegraphics[width=0.164\textwidth,height=1.635cm]{flag-5-869.jpg}}\hspace{-0.027in}
\subfigure{\includegraphics[width=0.164\textwidth,height=1.635cm]{flag-5-1008.jpg}}\hspace{-0.027in}\vspace{-0.117in}
\subfigure{\includegraphics[width=0.164\textwidth,height=1.635cm]{giraffe-10-399.jpg}}\hspace{-0.027in}
\subfigure{\includegraphics[width=0.164\textwidth,height=1.635cm]{giraffe-10-544.jpg}}\hspace{-0.027in}
\subfigure{\includegraphics[width=0.164\textwidth,height=1.635cm]{giraffe-10-751.jpg}}\hspace{-0.027in}
\subfigure{\includegraphics[width=0.164\textwidth,height=1.635cm]{hand-16-1438.jpg}}\hspace{-0.027in}
\subfigure{\includegraphics[width=0.164\textwidth,height=1.635cm]{hand-16-1801.jpg}}\hspace{-0.027in}
\subfigure{\includegraphics[width=0.164\textwidth,height=1.635cm]{hand-16-1929.jpg}}\hspace{-0.027in}\vspace{-0.095in}
\subfigure{\includegraphics[width=0.9\textwidth,height=0.405cm]{label.pdf}}
\caption{Visualization results of different methods on LaSOT dataset (from left to right and top to bottom: \emph{bicycle-9}, \emph{boat-4}, \emph{boat-17}, \emph{cat-3}, \emph{crocodile-10}, \emph{flag-5}, \emph{giraffe-10} and \emph{hand-16}).}
\label{fig_visual}
\end{figure*}

\textbf{OTB100 \cite{wu2015object}:} OTB100 is one of the most widely used benchmark datasets in tracking, which includes 100 fully labeled video sequences. The Area Under Curve (AUC) of success plot is employed to compare SiamCorners and nine state-of-the-art trackers including TADT \cite{li2019target}, GCT \cite{gao2019graph}, CIResNet-22-RPN \cite{zhang2019deeper}, RT-MDNet \cite{jung2018real}, SiamRPN \cite{li2018high}, DaSiamRPN \cite{zhu2018distractor}, FCAF \cite{han2019fully}, Ocean \cite{Ocean_2020}, ATOM \cite{danelljan2019atom} and DIMP \cite{bhat2019learning}. As shown in Fig. \ref{OTB}, our method performs better than the top tracker, such as ATOM, and slightly lower than its optimized version, DIMP. Our method outperforms the anchor-based trackers CIResNet22-RPN and DaSiamRPN with relative gains of 0.6\% and 1.3\% in terms of success score, respectively. Further, our method achieves a gain of 2.2\% in AUC over FCAF, while on par with offline anchor-free tracker Ocean. This shows that SiamCorners can achieve the state-of-the-art results by eliminating the anchor boxes.

\newcommand{\tabincell}[2]{\begin{tabular}{@{}#1@{}}#2\end{tabular}}
\begin{table*}[!t]
\centering
\caption{A comparison between the SiamCorners and state-of-the-art trackers on TrackingNet dataset in terms of AUC score, precision score (P) and normalized precision score (P$_{norm}$) }.

\fontsize{9pt}{9pt}\selectfont
\begin{tabular}{ccccccccccc}
\hline
\toprule
 &\tabincell{c}{ECO \\ \cite{danelljan2017eco}}& \tabincell{c}{SiamFC \\ \cite{bertinetto2016fully}}&\tabincell{c}{CFNet \\ \cite{valmadre2017end}}& \tabincell{c}{MDNet \\ \cite{nam2016learning}}&\tabincell{c}{GFS-DCF \\ \cite{xu2019joint}}&\tabincell{c}{UPDT \\ \cite{bhat2018unveiling}}& \tabincell{c}{DaSiamRPN \\ \cite{zhu2018distractor}}&\tabincell{c}{UpdateNet \\ \cite{zhang2019learning}}&\tabincell{c}{ATOM \\ \cite{danelljan2019atom}}&\tabincell{c}{SiamCorners}\\
\midrule
AUC(\%)& 55.4  & 57.1  & 57.8  &60.6  &60.9  &61.1  &63.8  &67.7 &70.3  &69.5\\
P(\%)& 49.2& 53.3& 53.3&56.5  &56.6  &55.7&59.1 &62.5 &64.8&64.7\\
P$_{norm}$(\%)& 61.8&66.3& 65.4& 70.5 &71.8  &70.2&73.3 &75.2 &77.1&76.3\\
\bottomrule
\end{tabular}
\label{trackingnet}
\end{table*}

\begin{table*}[!t]
\centering
\caption{A comparison between the SiamCorners and state-of-the-art trackers on NFS dataset.}
\fontsize{9pt}{9pt}\selectfont
\begin{tabular}{cccccccccccc}
\hline
\toprule
&\tabincell{c}{BACF \\ \cite{kiani2017learning}}  &\tabincell{c}{SRDCF\\ \cite{lukevzivc2018discriminative}}  &\tabincell{c}{FCNT\\ \cite{wang2015visual}}  &\tabincell{c}{DaSiam-\\RPN \cite{Zhu_2018_ECCV}}  &\tabincell{c}{HDT\\ \cite{qi2016hedged}}  &\tabincell{c}{MDNet\\ \cite{nam2016learning}}  &\tabincell{c}{MKCF\\ \cite{tang2018high}}  &\tabincell{c}{ECO\\ \cite{danelljan2017eco}}  &\tabincell{c}{CCOT\\ \cite{danelljan2016beyond}} &\tabincell{c}{SiamRPN++\\ \cite{li2019siamrpn++}}  &SiamCorners\\
\midrule
AUC(\%)&34.2  & 35.3  & 39.3  &39.5  &40.0  &42.5  &45.5  &47.0  &49.2  &49.2  &53.7\\
FPS &35.3  & 13.0  &3.0  &160.0  &-  &1.0  &30.0  &8.0  &0.3 &35.0  &42.0\\

\bottomrule
\end{tabular}
\label{NFS}
\end{table*}

\begin{figure}[!t]
\centering
\makeatletter\def\@captype{figure}\makeatother
\subfigure{
   \label{fig_Success_Plot}
   \includegraphics[width=4.1 cm]{lasot_success.pdf}}
\hspace{0.005 in}
\subfigure{
   \label{fig_Precision_Plot}
   \includegraphics[width=4.1 cm]{lasot_Normalized.pdf}}
\caption{Experimental results of the methods on the LaSOT dataset.}
\label{LaSOT}
\end{figure}

\textbf{UAV123 \cite{mueller2016benchmark}:} We evaluate the proposed SiamCorners method on the UAV123 dataset in Fig. \ref{UAV}. The dataset is collected from a low-altitude UAV, which contains 123 video sequences and more than 110K frames in total. In this dataset,  we compare our method with recent trackers, \emph{i.e.} SiamRPN++ \cite{li2019siamrpn++}, DaSiamRPN \cite{zhu2018distractor}, SiamRPN \cite{li2018high}, RTMDNet \cite{jung2018real}, ECO \cite{danelljan2017eco}, GCT \cite{gao2019graph}, ASLA \cite{jia2012visual}, SAMF \cite{li2014scale}, MEEM \cite{zhang2014meem}, and DSST \cite{danelljan2014accurate}. In terms of
the success score and precision score, SiamRPN++ achieves scores of 61.2\% and 80.6\%, respectively. Our approach outperforms the top tracker SiamRPN++ with relative gains of 0.2\% and 1.3\% in terms of the success score and precision score, respectively. It shows that our method achieves superior accuracy when compared with the recent methods based on anchor boxes.

In Fig. \ref{uav123}, we analyze the performance of the trackers under 12 challenging factors in tracking scenarios, such as lighting changes, scale variation, and aspect ratio change. These attributes are used to evaluate the performance of the trackers against different aspects. For clarity, the legend only shows the top 10 AUC scores in each plot. We observe that the proposed SiamCorners tracker handles well in scale variation, aspect ratio change, low resolution, fast motion, full occlusion, partial occlusion, out-of-view, viewpoint change, camera motion, and similar object. In particular, the proposed SiamCorners tracker has a 7.8\% higher AUC score than DaSiamRPN when tracking targets suffer from fast motion. These results show that the introduced SiamCorners method has the potential to solve the fast motion problem.

\textbf{TrackingNet \cite{muller2018trackingnet}:} TrackingNet is a recently proposed large-scale dataset in the wild, which provides more than 30K videos and 14 million bounding box annotations. We evaluate SiamCorners on its test split with 511 videos. Table \ref{trackingnet} shows the results of several state-of-the-art trackers, which are obtained by an online evaluation server. Among the compared approaches, ATOM achieves the best tracking result with an AUC score of 70.3\%. Our SiamCorners achieves an AUC score of 69.5\%, competitive with the other state-of-the-art trackers.

\textbf{LaSOT \cite{fan2019lasot}:} To further validate the effectiveness of the SiamCorners tracker, we evaluate our tracker on LaSOT \cite{fan2019lasot} dataset. LaSOT is a recently proposed tracking benchmark dataset, which contains 1400 sequences with an average video length of 2500 frames. We evaluate SiamCorners on its test split containing 280 videos. Furthermore, we compare our SiamCorners method with some representative trackers, that is,  SiamRPN \cite{li2018high}, MDNet \cite{nam2016learning}, VITAL \cite{song2018vital}, GCT \cite{gao2019graph}, ASRCF \cite{dai2019visual}, SiamFC \cite{bertinetto2016fully},
StructSiam \cite{zhang2018structured}, DSiam \cite{guo2017learning},  SINT \cite{tao2016siamese1}, and ECO \cite{danelljan2017eco}. As shown in Fig. \ref{LaSOT}, SiamCorners achieves the best performance in both success plot and normalized precision plot. Our approach outperforms MDNet with relative gains of 8.0\% and 9.3\% in terms of success score and normalized precision score, respectively.

Figure \ref{fig_visual} shows the visualization results of the trackers on the LaSOT dataset. Although these methods generally track targets accurately, they do not perform well under some challenging factors. The ECO and ASRCF based on correlation filters do not perform well in testing sequences with aspect ratio change (bicycle-9, cat-3, and giraffe-10). The MDNet and VITAL trackers are all built on classification component of CNN, which are prone to drift off the target when tracking targets undergo illumination variation (bicycle-9), partial occlusion (giraffe-10), and aspect ratio change (crocodile-10).  The SiamRPN and StructSiam trackers that follow the common Siamese architecture are less effective in addressing sequences with deformation (hand-16) and scale variation (boat-4, cat-3, and giraffe-10) attributes. In particular, the tracker SiamRPN based on anchor boxes cannot handle the scale variation problem well since it is difficult to handle the scale variation of the target by setting fixed ratio anchor boxes. Overall, the tracking results of the proposed SiamCorners performs well on most sequences. This qualitatively demonstrates the superiority of SiamCorners compared with other state-of-the-art methods.

\textbf{NFS30 \cite{kiani2017need}:} The NFS30 dataset contains 100 videos and a length of 380K frames in total, which are mainly captured by high frame rate cameras in real-world scenarios. We evaluate the SiamCorners tracker on the 30 FPS version of this dataset, which contains some challenging factors for fast-moving targets. From Table \ref{NFS}, we can observe that SiamCorners achieves $4.5\%$ higher than SiamRPN++ in terms of the AUC score. Furthermore, our method achieves a speed gain of 7 FPS compared with SiamRPN++, which verifies that the speed gain can be obtained by discarding anchors to reduce network parameters.

\subsection{Ablation Study}
\begin{table*}
\renewcommand\arraystretch{1.5}

	\centering
	\caption{Ablation study of the SiamCorners tracker on UAV123 and LaSOT.}
	\fontsize{9.7pt}{9.7pt}\selectfont

	\begin{tabular}{c|ccc|c|cc|cc|cc|c|c}
	    \hline
        Backbone & L3 & L4 & L5 & Offset & \tabincell{c}{MCP} & OCP & PPF &OPF &FRT &OUT &UAV123 & LaSOT \\

		\hline

		\cline{2-13}

		\multirow{6}{*}&\checkmark &  &  &\checkmark  &\checkmark  &&\checkmark&&\checkmark&&0.527  &0.380 \\

		~ &  &\checkmark  &  &\checkmark  &\checkmark  &&\checkmark&&\checkmark&&0.587 &0.432\\

		~ &  &  &\checkmark &\checkmark  &\checkmark   &&\checkmark&&\checkmark&&0.523  &0.373\\

		~ &\checkmark &\checkmark  &  &\checkmark  &\checkmark  &&\checkmark&&\checkmark&&0.556  &0.401\\

		~ &\checkmark &   &\checkmark &\checkmark  &\checkmark   &&\checkmark&&\checkmark&&0.539  &0.391\\

		~ {ResNet-50} &  &\checkmark &\checkmark &\checkmark  &\checkmark   &&\checkmark&&\checkmark& &0.601  &0.447\\

		\cline{2-13}

		\multirow{3}{*}{} &\checkmark &\checkmark   &\checkmark  &\checkmark  &  &\checkmark &\checkmark&&\checkmark& &0.593\  &0.455\\

		~ &\checkmark  &\checkmark &\checkmark &\checkmark  &\checkmark   &&&\checkmark&\checkmark& &0.600 &0.463\\

		~ &\checkmark  &\checkmark &\checkmark &\checkmark  & &\checkmark  &\checkmark&&&\checkmark &0.551  &0.411\\

		\cline{2-13}
		\multirow{2}{*}{} &\checkmark &\checkmark   &\checkmark  &\  &\checkmark  &&\checkmark&&\checkmark& &0.599\  &0.441\\

		~ &\checkmark  &\checkmark &\checkmark &\checkmark  &\checkmark   &&\checkmark&&\checkmark& &0.614  &0.480\\

		\hline
	\end{tabular}
\newline
\newline

\begin{tablenotes}

\item[*] Note: L3, L4, and L5 represent the \emph{conv3}, \emph{conv4}, and \emph{conv5} outputs of ResNet-50, respectively. Offset represents w-

hether to add an offset network to the convolutional network to predict the corner offsets.
\end{tablenotes}

\label{ablation}
\end{table*}

To verify the effectiveness of each component in the proposed approach, we perform an ablation study on the UAV123 and LaSOT datasets in terms of AUC score of success plot. To be fair, we use the same training datasets and parameter settings in the training step for different components.

\textbf{Backbone Architecture.} Different convolutional layers of the backbone network produce different feature representations and parameters, which will affect the speed, memory usage, and performance of the tracker. As shown in Table \ref{ablation}, we empirically observe that \emph{con4} branch achieves the best AUC score in a single-level feature, while shallower or deeper layers obtain performance drops. Furthermore, the performance will be further improved by combining \emph{conv4} branch and \emph{conv5} branch, while the performance is not improved when combining the other two branches. Finally, SiamCorners achieves the best result that produces a 0.614 AUC score on UAV123 through aggregating all three layers, which is 2.7\% higher than the single-layer baseline.

\textbf{Offset.} We investigate the impact of the offset network by excluding it from our SiamCorners framework. As shown in Table \ref{ablation}, the proposed offset operation gains 1.5\% improvement on UAV123 and 3.9\% improvement on LaSOT, which demonstrates the importance of offset network. This happens because the introduced offset network can alleviate the impact of the network stride and thus improve the performance.

\textbf{Corner Pooling.} We investigate the impact of the original corner pooling (OCP) and the modified corner pooling (MCP) on the final tracking results. Note that the original corner pooling operation only performs corner prediction on a single feature map. We feed the feature maps of the top and left boundary images to the top-left corner pooling of OCP and the feature maps of the bottom and right boundary images to the bottom-right corner pooling of OCP in the ablation experiment. Table \ref{ablation} shows that introducing a MCP gains 2.1\% improvement on UAV123 and 2.5\% improvement on LaSOT. This shows that the MCP operation is a key component in the SiamCorners framework. We further investigate the proposed penalty function (PPF) and the original penalty function (OPF) on the tracking performance. Using a PPF to select the final tracking box yields 1.4\% improvement on UAV123 and 1.7\% improvement on LaSOT. Finally, we analyze the impact of the four respective templates (FRT) and one unified template (OUT) on the tracking performance. We use the OCP operation since OUT only generates a single feature map as input for corner pooling. Further, FRT obtains a 6.3\% improvement on UAV123 and a 6.9\% improvement on LaSOT. This suggests the FRT adds some performance improvements to our tracker.

\textbf{Component-wise Time.} We investigate the component-wise computation time of the SiamCorners tracker on the UAV123 dataset, shown in Table \ref{ablation2}. The last residual block \emph{conv5} of ResNet-50 takes the most computational time. Since the number of parameters on the feature extractor will be increased by using a deeper network, which increases the running time. It is worth noting that the calculation time of corner pooling ranks third. In the following research, we can design a more efficient pooling operation to reduce the running time.
\begin{table}[!t]

\centering
\caption{Component-wise computation time of the SiamCorners tracker on UAV123 dataset}

\fontsize{9pt}{9pt}\selectfont
\begin{tabular}{ccccccccc}
\hline
\toprule
 Component
 &\tabincell{c}{C3\\}
 &\tabincell{c}{C4\\}
 &\tabincell{c}{C5\\}
 &\tabincell{c}{Off\\}
 &\tabincell{c}{Cor\\}
 &\tabincell{c}{Dec\\}
 &\tabincell{c}{Pos\\}\\
\midrule
Time (ms)
&\footnotesize 2.30
&\footnotesize 4.43
&\footnotesize 5.69
&\footnotesize 0.63
&\footnotesize 4.38
&\footnotesize 3.29
&\footnotesize 0.82 \\
\bottomrule
\end{tabular}
\label{ablation2}
\newline
\newline
\footnotesize{\leftline{Note: C3, C4, and C5 represent  calculation time that the testing images are}}

\footnotesize{\leftline{fed to \emph{conv3}, \emph{conv4}, and \emph{conv5} of ResNet-50, respectively. Off and Cor rep-}}

\footnotesize{\leftline{resent the offset predictions of the corners and corner pooling operations, r-}}

\footnotesize{\leftline{espectively. Dec and Pos represent decoding operations and post-processing}}

\footnotesize{\leftline{operations that correspond to Section \ref{decoding1} and Section \ref{fusion}, respectively.}}

\end{table}

\section{Conclusions}
In this paper, we propose a novel tracking architecture, SiamCorners, which does not require multi-scale testing or pre-defined anchor boxes with different aspect ratios. SiamCorners directly tracks the targets as a pair of corners. Our architecture integrates layer-wise features in deep networks. By exploring a new penalty function, we can constrain the variation of center, width, and height of the bounding boxes in adjacent frames. We evaluate SiamCorners on OTB100, UAV123, NFS30, LaSOT, and TrackingNet datasets, which demonstrates its effectiveness and efficiency. We hope that our work will lead to further exploration for anchor-free methods in visual tracking.
\section*{Acknowledgment}

This research was supported in part by Special Research project on COVID-19 Prevention and Control of Guangdong Province (Grant No. 2020KZDZDX1227), in part by the Shenzhen Research Council (Grant No. JCYJ20170413104556946 and No. JCYJ20170413105929681), in part by the Natural Science Foundation of China under Grant No. 62006060, Grant No. U2013210,  Grant No. 62002241 and Grant No. 61972112, in part by the Guangdong Basic and Applied Basic Research Foundation under Grant No. 2021B1515020088. This work was also supported in part by funding (Grant No. 2019-INT021) from Shenzhen Institute of Artificial Intelligence and Robotics for Society (AIRS).

\bibliographystyle{IEEEtran}