\section{Introduction}

\IEEEPARstart{S}{ynthetic} Aperture Radar (SAR) has been widely used in Earth observation applications due to its capability to work under all weather and daylight conditions. The semantic segmentation of Polarimetric SAR (PolSAR) images, namely the pixel-wise classification of PolSAR images according to ground surface types, is beneficial to a large number of remote sensing applications (e.g., urban area management, disaster monitoring and land-cover mapping).

In recent years, with the rise of convolutional neural networks (CNNs) many methods have been developed for the semantic segmentation of natural images \cite{zhao2017pspnet}\cite{chen2018deeplabv3+} and remote sensing images \cite{ding2020lanet}\cite{mou2019relation}. However, a limited number of studies has been conducted on the semantic segmentation of PolSAR images based on deep CNNs \cite{geng2020multi}. There are two major barriers: i) PolSAR images contain intense speckle noise due to the coherent imaging mechanism of PolSAR systems. This speckle noise is a severe challenge for the automatic segmentation algorithms. ii) Large amount of training data is required to train an effective deep CNN. The manual annotation of SAR data is not only labor-intensive but also difficult, since the ground objects in SAR images can hardly be recognized by human observation without assisting data and expert knowledge \cite{duan2018multi}.

Most of the literature works on the semantic segmentation of SAR images are designed under the assumption that there are not enough labeled samples to train a CNN. In \cite{li2018novel}, a relatively small portion of pixels of the image to be classified are used for the training phase. The Fully Connected Network (FCN) has been used in a pipeline that contains sliding-crop and sparse coding operations. In \cite{wang2018hierarchical}, the FCN is combined with the sparse and low-rank subspace representations to alleviate the problem of insufficient training data. In \cite{duan2018multi}, a multi-scale CNN has been proposed for the semantic segmentation of SAR images. To solve the problem of lacking ground truth segmentation maps, it uses image scene labels for the training. In \cite{cao2019pixel}, the FCN is used in the complex domain of SAR data to include the phase information. In \cite{wu2019polsar}, extra datasets are used to pre-train CNN models for the semantic segmentation of SAR images. Since all these studies have been conducted on small datasets (most datasets contain only a single image), their objectives are mainly to reduce the dependence on training data keeping high generalization capabilities.

\begin{figure*}[htbp]
    \centering
    \includegraphics[width=1\linewidth]{pics/FlowChart.png}
    \caption{Atchitecture of the proposed Multi-path Residual Network (MP-ResNet).}
    \label{FigOverview}
\end{figure*}

\begin{figure}[htbp]
    \centering
    \includegraphics[width=8cm]{pics/Deconv.png}
    \caption{Feature deconvolution block of the proposed MP-ResNet.}
    \label{FigDeconv}
\end{figure}

Recently, with the emergence of several benchmark datasets, new developments on the CNN-based semantic segmentation of SAR images have been proposed. In \cite{mohammadimanesh2019new}, a 'encoder-decoder' CNN network with inception modules and skip connections is introduced for the semantic segmentation of wetland PolSAR images. In \cite{yue2020MANet}, a multi-scale attention based FCN (MANet) is presented combining multi-scale feature extraction and the attention mechanism. In \cite{wang2020hrsarnet}, a small yet efficient network (HR-SARNet) is proposed for the semantic segmentation of high-resolution SAR images. However, most of these literature studies are based on shallow CNNs to avoid over-fitting problems. Since shallow CNNs are not powerful enough to extract high-level semantic information, their accuracy is limited.
Differently from previous studies, this work aims to improve the semantic segmentation of high-resolution PolSAR images by proposing a novel network architecture, the Multi-path Residual Network (MP-ResNet). This network enables a multi-scale modeling of high-level semantic features through its parallel branches, which strengthens the learning of local discriminative features and reduces the effects of speckle noises. This network is based on a large benchmark dataset on PolSAR images made available by the Gaofen committee. Ablation studies show that the proposed approach achieves an large increase of 0.36\% and 0.64\% in terms of average OA and fwIoU, respectively, compared to the baseline FCN. It also surpasses several literature works in the comparative experiments.

\section{Methodology}

The key to improve the semantic segmentation of PolSAR images is to learn discriminative features from a larger image context, so that the effects of speckle noises can be mitigated. Based on this idea, we propose the Multi-path Residual Network (MP-ResNet) shown in Fig.\ref{FigOverview}. This network exploits a baseline fully convolutional network (FCN). It is inspired by the multi-path convolution design in high-resolution network (HRNet) \cite{wang2020hrnet} and the multi-scale feature fusion design in LinkNet \cite{chaurasia2017linknet}. In this section we introduce in details the design of the proposed MP-ResNet.

\subsection{Quick Embedding of Shallow Features}
PolSAR data generally contain 4 channels which are related to the 4 combination of linear polarizations (i.e. the HH and VV co-polarized channels and the HV and VH cross-polarized channels). The corresponding 4 images are stacked and given as input to a CNN. Since CNNs are capable of extracting semantic features from raw input data, no extra filtering operations have been applied to the input images. Instead, we merely apply maximum-suppression and normalization operations to stabilize and squeeze the value range of input data.

It is well-known that the signal at pixel level in PolSAR images is generally very noisy. The most discriminative features in these images are the signal values and texture patterns in local areas.  This makes it necessary to down-sample the features in early layers of a CNN and embed the features from a wider range. Considering this, we adopt the ResNet as the backbone feature extraction network. The first layers of ResNet are strided convolutions followed by a pooling layer, which quickly decreases the scaling rate of features to 1/4. Additionally, while 3$\times$3 convolutional kernels are widely used in CNNs to avoid over-fitting problems, the kernel size of the first convolutional layer in the ResNet is 7$\times$7. This helps the network to aggregate features from a larger pixel neighbourhood and alleviate the effect of speckle noises. After the first two residual modules of ResNet, the scaling rate of features is quickly decreased to 1/8. Considering that most of PolSAR images are not high-resolution and thus do not contain rich spatial details, this feature size is large enough for modeling the regional land-cover information. Therefore, 1/8 is set as the fundamental scaling rate for the aggregation of context information in the proposed network.

\subsection{Multi-path Semantic Information Embedding}
The size of valid receptive fields (VRFs) is known to be crucial to the embedding of context information in CNNs \cite{zhang2018context}. In the semantic segmentation of remote sensing data, the size of VRFs determines the spatial range from which the CNN can exploit discriminative features, which is related to the granularity of semantic segmentation results \cite{ding2020twostage}. Although the serial connection of pooling and strided convolutions can greatly enlarge the VRFs, it also brings the problem of losing spatial information. Therefore, how to simultaneously enlarge VRFs and preserve spatial information is one of the most crucial bottleneck problems in semantic segmentation tasks. PSPNet \cite{zhao2017pspnet} and DeepLab \cite{chen2018deeplabv3+} managed to enlarge VRFs without severe loss of spatial information through the use of additional context aggregation modules in the late layers of the CNNs. However, the context information is aggregated through pooling and dilated convolutions in these designs, which are less effective than stacked convolutional layers. The dilated convolutions may also cause gridding effects and enlarge computational costs.

Alternatively, HRNet presents a multi-path architecture that organizes multi-scale convolutional layers in a parallel manner. The highest scaling rate of HRNet is 1/32, which ensures a large VRF of the network. However, the serial convolutions in its parallel branches greatly increases the computational costs of this network. They may also cause over-fitting problems on small datasets. In addition, in HRNet the multi-scale feature branches are concatenated together to generate the semantic segmentation results, thus the features are not fully fused and utilized.

Inspired by the parallel feature embedding design of the HRNet, The proposed MP-ResNet uses ResNet as the feature extraction network but has 2 additional encoding branches after the second convolutional block. Therefore, the features are encoded both forwardly (size remain the same) and downstream-wise (size reduced). The same residual blocks are duplicated into the parallel branches so that each encoding branch has the same amount of convolutional layers. In this way, each encoding branch contains rich semantic information.

Differently from the HRNet, the parallel  branches in MP-ResNet focus on the embedding of high-level features (the lowest scaling rate is 1/8). Although the parallel embedding of larger feature maps are potentially feasible in other applications, for the semantic segmentation of PolSAR images our objective is to exploit discriminative features from a wider range. In this way, the segmentation results become less sensitive to pixel noise. The computational costs of MP-ResNet is only slightly higher than ResNet and far less than the HRNet (see discussion in Section \ref{sec.results}). Another significant difference of the MP-ResNet is that a decoder network is employed to fuse the features learned from its parallel branches.

\subsection{Fusion of Multi-scale Features}
Decoder networks are commonly used in semantic segmentation to recover the spatial details of encoded features. A common design is to concatenate or add the encoded features with the features from early layers of the encoder networks (e.g. UNet\cite{ronneberger2015unet} and SegNet\cite{badrinarayanan2017segnet}). Although this design can aggregate spatial information from low-level features, it also introduces
redundant information (minor details and noises). For the segmentation of high-resolution PolSAR images this problem can be critical. However, in the proposed MP-ResNet, the fused multi-scale features are the ones encoded by the parallel branches of the encoder. These features contain rich semantic information, thus their fusion does not lead to noise problems. A feature deconvolution module \cite{chaurasia2017linknet} is introduced to enlarge the spatial size of features from higher branches before the fusion. Fig.\ref{FigDeconv} shows this spatial deconvolution block. It adopts a channel-wise 'Bottleneck' design to reduce the computation and refine the crucial semantic information. In this way, the multi-branch features are fused in a coarse-to-fine manner in the decoder.

\section{Experimental Results}\label{sec.results}
\subsection{Dataset Descriptions and Evaluation Metrics}
The experiments dataset of this study are developed on the Gaofen dataset provided by the '2020 Gaofen contest on automated high-resolution earth observation image interpretation'. The PolSAR images are collected from the Gaofen-3 satellite. Their ground sampling distance is between 1m and 3m. The ground truth maps are annotated according to 5 land-cover types: background, built-up area, vegetation, water and bare soil. The accessible training data are 500 pairs of PolSAR images and label maps each with 512$\times$512 pixels. The testing data are not visible to users, but a scoring system is provided to evaluate the uploaded algorithms.

We adopts 3 metrics for the evaluation of semantic segmentation results. They are overall accuracy (OA), F1 score and frequency weighted intersection over union (fwIoU). OA is the ratio of the number of correctly classified pixels among all pixel numbers. F1 is calculated as:
    \begin{equation}
    \rm F1=2 \cdot {\frac{precision \cdot recall}{precision+recall}}
    \label{F1}
    \end{equation}
FwIoU is the evaluation metric suggested by the contest organizer. It is calculated as:
    \begin{equation}
        fwIoU = \frac{1}{\sum_{i=0}^{N}\sum_{j=0}^{N}}\sum_{i=0}^{N}\frac{\sum_{j=0}^{N}s_{ij}s_{ii}}{\sum_{j=0}^{N}s_{ij}+\sum_{j=0}^{N}s_{ji}-s_{ii}}
    \end{equation}
where N is the number of total classes, $s_{ij}$ denotes the number of i-th class pixels that are classified into the j-th class.

\subsection{Multi-fold Ablation Study}\label{sec.ablation}
To quantitatively evaluate the improvement of the proposed MP-ResNet over the baseline method (FCN), we conduct a multi-fold ablation study on the Gaofen dataset. The training and validation sets are randomly divided from all available training data (500 image pairs) with a numeric ratio of 9:1. In this way, a total of 10 training and validation sets are obtained. To reduce the effects of random factors, the ablation study has been conducted on all the 10 training-validation sets. The results are reported in table \ref{Table.Ablation}. Compared to the baseline method (FCN), the proposed MP-ResNet shows average improvements of 0.36\%, 0.54\% and 0.64\% in OA, mean F1 and fwIoU, respectively. The improvements have also been verified on the contest test set which is not directly available. According to the fwIoU scores obtained from the contest, the proposed MP-ResNet has an increase of 1.23\% in fwIoU compared to the FCN.

\begin{table}[htbp]
    \centering
    \caption{Results of the multi-fold ablation study.}
    \resizebox{1\linewidth}{!}{%
        \begin{tabular}{c|ccc|ccc}
        \toprule
        \multirow{2}*{Datasets} & \multicolumn{3}{c|}{FCN (Baseline)} & \multicolumn{3}{c}{MP-ResNet (Proposed)}\\
        \cline{2-7}
         & OA(\%) & mF1(\%) & fwIoU(\%) & OA(\%) & mF1(\%) & fwIoU(\%)\\
        \hline
         Val 1 & 93.08 & 91.33 & 88.14 & 93.95(\textcolor{green}{+0.87}) & 92.47(\textcolor{green}{+1.14}) & 89.63(\textcolor{green}{+1.49}) \\
         Val 2 & 93.41 & 91.81 & 89.04 & 93.88(\textcolor{green}{+0.47}) & 92.58(\textcolor{green}{+0.77}) & 89.88(\textcolor{green}{+0.84}) \\
         Val 3 & 91.70 & 88.20 & 86.35 & 91.86(\textcolor{green}{+0.16}) & 88.24(\textcolor{green}{+0.04}) & 86.71(\textcolor{green}{+0.36}) \\
         Val 4 & 90.08 & 85.52 & 84.65 & 90.22(\textcolor{green}{+0.14}) & 86.11(\textcolor{green}{+0.59}) & 84.58(\textcolor{green}{+0.24}) \\
         Val 5 & 93.06 & 89.64 & 88.34 & 93.45(\textcolor{green}{+0.39}) & 90.46(\textcolor{green}{+0.82}) & 88.97(\textcolor{green}{+0.63}) \\
         Val 6 & 91.96 & 89.28 & 86.37 & 92.85(\textcolor{green}{+0.89}) & 90.76(\textcolor{green}{+1.48}) & 87.81(\textcolor{green}{+1.44}) \\
         Val 7 & 92.97 & 89.22 & 88.25 & 93.33(\textcolor{green}{+0.36}) & 89.20(\textcolor{red}{-0.02}) & 88.90(\textcolor{green}{+0.65}) \\
         Val 8 & 91.43 & 87.21 & 86.09 & 91.60(\textcolor{green}{+0.17}) & 87.58(\textcolor{green}{+0.37}) & 86.52(\textcolor{green}{+0.43}) \\
         Val 9 & 93.11 & 90.66 & 88.66 & 93.24(\textcolor{green}{+0.13}) & 90.72(\textcolor{green}{+0.06}) & 88.91(\textcolor{green}{+0.25}) \\
         Val 10 & 93.40 & 91.59 & 89.50 & 93.62(\textcolor{green}{+0.22}) & 91.76(\textcolor{green}{+0.17}) & 89.85(\textcolor{green}{+0.35}) \\
         Test set & - & - & 69.42 & - & - & 70.65 (\textcolor{green}{+1.23}) \\
        \hline
         \multicolumn{4}{c|}{Average Improvements} & +0.36 & +0.54 & +0.64 \\

        \hline
        \end{tabular}}
   \label{Table.Ablation}
\end{table}

\subsection{Comparative Experiments}
To further assess the improvements brought by the proposed MP-ResNet, we compared it with several literature works. Apart from the baseline method FCN, several well-established methods in the field of semantic segmentation have been considered, including the SegNet \cite{badrinarayanan2017segnet}, the UNet \cite{ronneberger2015unet}, the PSPNet \cite{zhao2017pspnet} and the DeepLabv3+ \cite{chen2018deeplabv3+}. Since the proposed MP-ResNet is inspired by the multi-path architecture of HRNet\cite{wang2020hrnet} and the feature fusion design of LinkNet\cite{chaurasia2017linknet}, we also included these networks in comparison. Moreover, several literature methods presented for the semantic segmentation of SAR images have also been tested, including the mult-scale FCN (MS-FCN) in \cite{wu2019polsar}, the Inception FCN (Inc-FCN) in \cite{mohammadimanesh2019new} and the HR-SARNet \cite{wang2020hrsarnet}. The training and validation sets used in the comparison are these related to the first row of table \ref{Table.Ablation}. For fairness, the same parameter settings have been applied during the training process (e.g. training epochs, batch size, learning rate).

Table \ref{Table.Compare} presents the quantitative results obtained on the Gaofen dataset. Due to the effects of speckle noise, the performances of shallow networks (SegNet, HR-SARNet, Inc-FCN and UNet), which rely heavily on low-level features, is unsatisfactory. Although there is also a multi-scale design in the MS-FCN, there is no enhancement in its branches, thus its accuracy is lower than that of the FCN. The simple FCN without any sophisticated design ranks at the 3rd place. Although HRNet shows better performance in the semantic segmentation of VHR optical remote sensing images \cite{zhang2020multi}, its size is too large to be fully trained on the Gaofen dataset, thus its accuracy is far lower than that of the FCN. LinkNet also has a deep encoding network (ResNet34) but it has skip connections with low-level features, which introduce noise and degrades its accuracy. The PSPNet has an additional multi-scale average pooling head compared with the FCN, which does not prove to be effective on the PolSAR dataset. The multi-scale atrous convolution head in DeepLabV3+ increases the mean F1, OA and fwIoU by 0.57\%, 0.30\% and 0.36\%, respectively, compared to the FCN. The proposed MP-ResNet achieves the best accuracy in nearly every metrics except the F1 of the 'others' class. Its improvements over the DeepLabV3+ are 0.57\%, 0.57\% and 1.13\% in mean F1, OA and fwIoU, respectively. This proves the effectiveness of the proposed network with which we have won the 2nd place in the 'Gaofen Challenge' contest.

\begin{table*}[htbp]
    \centering
    \caption{Quantitative Results of the comparative study.}
    \resizebox{1\linewidth}{!}{%
        \begin{tabular}{c|cccccc|ccc}
        \toprule
        \multirow{2}*{Methods} & \multicolumn{6}{c|}{F1(\%)} & \multirow{2}*{mean F1(\%)} & \multirow{2}*{OA(\%)} & \multirow{2}*{fwIoU(\%)$\blacktriangledown$} \\
        \cline{2-7}
         & Water & Built-up area & Industrial area & Grassland & Barren & Others & & & \\
        \hline
        SegNet\cite{badrinarayanan2017segnet} & 45.78 & 24.13 & 36.42 & 65.80 & 0.00 & 0.00 & 28.69 & 45.79 & 36.25\\
        HR-SARNet\cite{wang2020hrsarnet} & 73.10 & 76.29 & 57.31 & 78.67 & 0.00 & 0.10 & 47.58 & 74.01 & 62.87\\
        Inc-FCN\cite{mohammadimanesh2019new} & 83.62 & 74.90 & 52.92 & 78.58 & 4.41 & 15.05 & 51.58 & 73.74 & 62.97\\
        MS-FCN\cite{wu2019polsar} & 69.68 & 77.14 & 60.02 & 79.13 & 0.00 & 0.01 & 47.66 & 74.63 & 63.09\\
        UNet\cite{ronneberger2015unet} & 83.58 & 77.54 & 59.63 & 79.67 & 0.02 & 22.36 & 53.80 & 76.14 & 65.43\\
        HRNet\cite{wang2020hrnet} & 91.30 & 90.19 & 85.09 & 90.34 & 73.96 & 92.98 & 87.31 & 90.02 & 83.41\\
        LinkNet\cite{chaurasia2017linknet} & 90.18 & 91.63 & 88.43 & 91.11 & 74.95 & 92.90 & 88.20 & 91.17 & 85.19\\
        PSPNet\cite{zhao2017pspnet} & 92.71 & 93.36 & 89.71 & 93.12 & 84.03 & 94.91 & 91.31 & 92.98 & 88.05\\
        FCN\cite{long2015fcn} & 91.89 & 93.53 & 90.26 & 93.20 & 84.79 & 94.28 & 91.33 & 93.08 & 88.14\\
        DeepLabV3+\cite{chen2018deeplabv3+} & 93.44 & 93.67 & 90.37 & 93.41 & 85.15 & \textbf{95.32} & 91.90 & 93.38 & 88.50\\
        MP-ResNet (proposed) & \textbf{94.85} & \textbf{94.29} & \textbf{90.78} & \textbf{93.93} & \textbf{86.03} & 94.96 & \textbf{92.47} & \textbf{93.95} & \textbf{89.63}\\
        \hline
        \end{tabular}}
   \label{Table.Compare}
\end{table*}

Fig.\ref{Fig.Compare} shows a comparison of the segmentation results on several sample areas. Due to its multi-path modeling and multi-scale feature fusion design, the proposed MP-ResNet is capable of modeling context information from a wider image range. Therefore, some critical areas for other networks are correctly segmented and the detected object boundaries are more continuous.

\begin{figure}[thpb]
\centering
    {\includegraphics[height=0.5cm]{Compare/ColorBar.png}}\\
    \setlength{\tabcolsep}{1.0pt}
    \begin{tabular}{>{\centering\arraybackslash}m{1.7cm}>{\centering\arraybackslash}m{1.7cm}>{\centering\arraybackslash}m{1.7cm}>{\centering\arraybackslash}m{1.7cm}>{\centering\arraybackslash}m{1.7cm}}
        \includegraphics[width=1.7cm]{Compare/55_img.png} &
        \includegraphics[width=1.7cm]{Compare/55_gt_a.png} &
        \includegraphics[width=1.7cm]{Compare/55_Res34_a.png} &
        \includegraphics[width=1.7cm]{Compare/55_Deeplab_a.png} &
        \includegraphics[width=1.7cm]{Compare/55_MPResNet_a.png}\\
        \includegraphics[width=1.7cm]{Compare/142_rgb.png} &
        \includegraphics[width=1.7cm]{Compare/142_gt_a.png} &
        \includegraphics[width=1.7cm]{Compare/142_Res34_a.png} &
        \includegraphics[width=1.7cm]{Compare/142_Deeplab_a.png} &
        \includegraphics[width=1.7cm]{Compare/142_MPResNet_a.png}\\
        \includegraphics[width=1.7cm]{Compare/291_rgb.png} &
        \includegraphics[width=1.7cm]{Compare/291_gt_a.png} &
        \includegraphics[width=1.7cm]{Compare/291_Res34_a.png} &
        \includegraphics[width=1.7cm]{Compare/291_Deeplab_a.png} &
        \includegraphics[width=1.7cm]{Compare/291_MPResNet_a.png}\\
        \includegraphics[width=1.7cm]{Compare/133_rgb.png} &
        \includegraphics[width=1.7cm]{Compare/133_gt_a.png} &
        \includegraphics[width=1.7cm]{Compare/133_Res34_a.png} &
        \includegraphics[width=1.7cm]{Compare/133_Deeplab_a.png} &
        \includegraphics[width=1.7cm]{Compare/133_MPResNet_a.png}\\
        PolSAR images (pseudo color) & Ground truth & FCN & Deeplabv3+ & MP-ResNet (Proposed)\\
    \end{tabular}
    \caption{Comparison of segmented maps obtained by different methods on sample testing areas.}
    \label{Fig.Compare}
\end{figure}

\begin{table}[htbp]
\centering
    \caption{Comparison of model size and computational cost expressed in terms of params (Mb) and FLOPS (Gbps), respectively.}
    \begin{tabular}{l|cc}
        \toprule
        Methods & Params (Mb) & FLOPS (Gbps) $\blacktriangledown$\\
        \hline
        HR-SARNet\cite{wang2020hrsarnet} & 0.06 & 6.12 \\
        MS-FCN\cite{wu2019polsar} & 23.54 & 21.73 \\
        LinkNet \cite{chaurasia2017linknet} & 21.65 & 27.76\\
        FCN \cite{long2015fcn} & 21.35 & 90.97 \\
        HRNet\cite{wang2020hrnet} & 65.85 & 93.64 \\
        DeeplabV3+ \cite{chen2018deeplabv3+} & 22.29 & 94.84  \\
        PSPNet\cite{zhao2017pspnet} & 22.73 & 95.55\\
        MP-ResNet (proposed) & 54.97 & 115.93 \\
        UNet \cite{ronneberger2015unet} & 9.16 & 221.68 \\
        Inc-FCN\cite{mohammadimanesh2019new} & 6.14 & 329.37\\
        SegNet\cite{badrinarayanan2017segnet} & 3.51 & 341.06 \\
        \bottomrule
    \end{tabular}
    \label{Table.Calculations}
\end{table}

Table \ref{Table.Calculations} presents the sizes and computational costs of the compared methods. The floating point operations per second (FLOPS) are calculated based on the input size [4, 512, 512] of the images in the Gaofen dataset. The literature methods for the semantic segmentation of SAR images (HR-SARNet and Inc-FCN) are generally shallow, thus and their sizes are relatively small. UNet, Inc-FCN and SegNet need the most FLOPS since they apply many convolution operations on the early-layer features. Although the parameter size of the MP-ResNet is much larger than the FCN, its FLOPS do not increase significantly.

\section{Conclusion}
The semantic segmentation of PolSAR images is challenging due to the intense speckle noise and the lack of large training datasets. Taking advantage of the open dataset from the Gaofen contest, we propose a Multi-Path Residual Network (MP-ResNet) for the semantic segmentation of high-resolution PolSAR images. Compared to the baseline FCN, the MP-ResNet has three parallel semantic embedding branches to strengthen the aggregation of context information. It also adopts a multi-scale feature fusion design in its decoder to take advantage from each encoding branch. As a result, the VRF of the MP-ResNet is significantly enlarged, thus allowing the aggregation of discriminative features from a large range and alleviating the impact of noise. The multi-fold ablation study conducted on the Gaofen dataset has proved the effectiveness of our designs. The comparative experiments with several state-of-the-art methods show that the proposed method has a significant improvement in all accuracy metrics.

Since the objectiveness of this work is to propose a general architecture for the semantic segmentation of PolSAR images, we did not add sophisticated designs. However, in future works we plan to combine additional context aggregation modules to improve the performance (e.g., dilated convolutions, attention mechanisms).

\ifCLASSOPTIONcaptionsoff
  \newpage
\fi

\bibliographystyle{IEEEtran}


\section{Introduction}
Sequence labeling is the problem of labeling each token of an input sequence to produce an output sequence of labels. It is a fundamental problem in many fields such as machine learning, natural language processing, computer vision, and speech recognition.

To solve sequence labeling, the linear chain conditional random field (CRF) \cite{lafferty2001conditional} is developed which assumes dependency relations between neighboring labels conditioned on the input sequence. Using hand-crafted feature representations of the input sequence, linear chain CRFs have been an important benchmark for many real world sequence labeling tasks \cite{lafferty2001conditional,sha2003shallow,gimpel2011part}.  The recent renaissance in neural networks has led to significant progress on automatic feature extraction with deep neural networks. By modeling the potential functions in CRFs with deep neural networks, neural CRF models achieve significant performance boost over traditional feature-based CRFs on many tasks \cite{collobert2011natural,chen2015long,ma2016end,lample2016neural}.

In many applications, however, there exist long-range dependency relations between output labels. To model such long-range dependencies, the potential function of a CRF can be factorized by cliques containing more than two neighboring labels, resulting in a high order CRF. A major disadvantage of high order CRFs is that the computational complexity of training and decoding is exponential in the order.

In this paper, we present a novel latent variable model called the Neural Latent Dependency Model (NLDM) that solves the long-range dependency problem by modeling the dependency relations between output labels with a latent tree structure. Similar to neural CRFs, our model employs a deep neural network (a Bi-LSTM) for automatic feature extraction from the input sequence.
Unlike high order CRFs, our model does not place an upper bound to the dependency length between labels. By not allowing edge crossing in the latent tree structure, our model can be trained in an end-to-end way in polynomial time with respect to the sequence length and the label number. Although exact decoding (predicting the output labels given an input sequence) may not be tractable in our model, we propose an efficient approximate decoding algorithm.

To empirically verify the benefits of capturing long-range dependencies, we generated synthetic datasets from a neural \emph{infinite order} generative model and show that our model consistently outperforms the neural CRF on these datasets. We then empirically show that our model also outperforms the neural CRF and other baselines on nine real datasets of sequence labeling.

\section{Background}
\label{sec:bg}
\subsection{Sequence Labeling with Linear Chain CRFs} \label{sec:2.1}
Given an input sequence $\mathbf{x}=x_1, x_2, ..., x_n$ of length $n$ and an output label sequence $\mathbf{y}=y_1, y_2, ..., y_n$ where $y_i$ is the label of input token $x_i$, a conditional random field model \cite{lafferty2001conditional} is used to assign a score $\phi(\mathbf{x}, \mathbf{y})$ to the pair $(\mathbf{x}, \mathbf{y})$. Linear chain CRFs make the so-called first-order assumption that given the input sequence $\mathbf{x}$ and a label $y_i$, the labels to the left of $y_i$ are independent of the labels to the right of $y_i$. Therefore, the score $\phi (\mathbf{x}, \mathbf{y})$ can be factorized as follows,

\[\phi(\mathbf{x}, \mathbf{y})=\sum_{i=1}^{n} \left(\phi (\mathbf{x}, y_i, i) + \psi(y_{i-1}, y_{i})  \right)\]
In the equation, $\phi(\mathbf{x}, y_i, i)$ is the emission score, which is the score of a particular label $y_i$ given sequence $\mathbf{x}$ at position $i$.
$\psi(y_{i-1}, y_{i})$ is the transition score, which is the score of the transition from label $y_{i-1}$ to label $y_i$. With $m$ possible labels for each token, the transition score is parameterized with matrix $\mathbf{t} \in \mathbb{R}^{m\times m}$.

To train the model, we typically aim to minimize the regularized negative conditional log-likelihood of the training data, which is,
\begin{align*}
& J(\boldsymbol{\Theta}) = -\sum_{\alpha=1}^{N} \log P(\mathbf{y}^{\alpha} | \mathbf{x}^{\alpha}) + \Omega ||\mathbf{\Theta}||_2^2 \\
&= -\sum_{\alpha=1}^{N} \left( \Phi(\mathbf{x}^{\alpha}, \mathbf{y}^{\alpha})-\log \mathbf{Z} (\mathbf{x}^{\alpha}) \right) + \Omega ||\mathbf{\Theta}||_2^2
\end{align*}
where $\Omega$ is a hyper-parameter, and $\mathbf{Z}(\mathbf{x}^{\alpha})$ is the partition function, which can be efficiently calculated using the forward algorithm with time complexity $O(n \times m^2)$.

\subsection{Sequence Labeling with Deep Neural Networks} \label{sec:2.2}
Linear chain CRFs have been extremely popular in sequence labeling for their simplicity and efficiency.
However, hand-crafted feature engineering in CRFs requires domain-specific knowledge and task-specific resources.

Recent progress in deep learning enables effective automatic feature extraction for many sequence labeling tasks in natural language processing, including word segmentation \cite{chen2015long}, named entity recognition \cite{ma2016end,lample2016neural}, POS tagging \cite{ma2016end} and so on. All these approaches employ a non-linear deep neural network that takes as input pre-trained word embeddings \cite{mikolov2013efficient,pennington2014glove} and outputs a continuous representation for each word in the input sequence that captures not only the word information but also information of its context.

Since the neural network already captures sufficient context information of the input tokens, it can be used with or without the original CRF model for sequence labeling. Mathematically, after the hidden representation $\mathbf{h}_i \in \mathbb{R}^{D_{h}}$ is produced by the deep neural network for each token $x_i$ at position $i$, the emission score can be computed using a bilinear function,
\[\phi(\mathbf{x}, y_i, i) = \mathbf{h}_{i}^{T} \mathbf{W}_{e} \mathbf{t}_{y_{i}}\]
where $\mathbf{t}_{y_{i}} \in \mathbb{R}^{D_l}$ is the embedding for label $y_i$, $\mathbf{W}_{e} \in \mathbb{R}^{D_h\times D_l}$ is a weight matrix
The label for each token can then be predicted independently by maximizing the following conditional probability
\[P(\mathbf{y}|\mathbf{x}) = \prod_{i=1}^{n} P(y_i|\mathbf{x})=\prod_{i=1}^{n} \frac{e^{\phi(\mathbf{x}, y_i, i)}}{\sum_{y'} e^{\phi(\mathbf{x}, y', i)}}\]

We call this approach neural softmax. On the other hand, we can also explicitly model label correlation by using a CRF on top of the output representation of the neural network, resulting in the so-called neural CRF model with the following conditional probability:
\[P(\mathbf{y}|\mathbf{x})=\frac{e^{\sum_{i=1}^{n} \big( \phi(\mathbf{x}, y_i, i) + \psi (y_{i-1}, y_{i})\big)}}{\sum_{\mathbf{y'}}e^{\sum_{i=1}^{n} \big( \phi(\mathbf{x}, y'_i, i) + \psi (y'_{i-1}, y'_i)\big)}}\]

In the literature, for the neural softmax model, \newcite{collobert2011natural} proposed a simple yet effective feed-forward neural network for predicting the label for each token by encoding its context information within a fixed window size. For the neural CRF model, \newcite{huang2015bidirectional} was the first to show the promise of applying a (bi-directional) Long-Short-Term-Memory network (LSTM) together with a CRF layer in sequence labeling tasks. \newcite{lample2016neural} enhanced word embeddings with character-based word representations produced by a Bi-LSTM network and showed its success in NER tasks for many languages. \newcite{ma2016end} employed a convolutional neural network (CNN) to represent the morphological information of a word from character embeddings in POS tagging and NER tasks.

\subsection{High Order CRFs}
Although linear chain CRFs have been shown to be effective in sequence labeling, they only model dependencies between neighboring labels and therefore have limited capacity in capturing long-range dependencies. High order CRFs extend linear chain CRFs by explicitly modeling all the label dependencies no longer than their order. However, the computational complexity of training and inference in high order CRFs grows exponentially in the order \cite{ye2009conditional}. For example, the time complexity of exact inference in a $k$-th order CRF is $O(n\times m^{k+1})$. When $k$ is large, the computational cost is not affordable.

Several approaches have been proposed to alleviate this problem under different conditions. \newcite{ye2009conditional} showed that by restricting the number of distinct label sequences used in the features, the time complexity for decoding can be reduced to polynomial. \newcite{qian2009sparse} proposed a sparse higher order model in which learning and inference with high order features can be efficient under the pattern sparsity assumption.

Compared with these approaches,
our model achieves efficient learning and inference without any assumption on the features or the properties of the output labels.

\subsection{Dependency Grammars and Its Variants}
\label{sec:dg}
Our approach is closely related to dependency grammars.
A dependency grammar is a simple yet effective way of covering syntactic relations between words in a sentence in the natural language processing field \cite{nivre2005dependency}. The syntactic relations form a dependency parse tree of the sentence. Given a sentence, discovering its parse tree is called dependency parsing. There are two types of dependency grammars that differ in whether the dependencies in the parse tree are allowed to cross or not. Crossing is allowed by non-projective dependency grammars and is not allowed by projective dependency grammars.
Note that a dummy root node is typically added at the beginning of the sentence which connects to the actual root word of the parse tree. If multiple edges are allowed to emit from the dummy root node, then we have a dependency forest over the sentence.

The quality of a dependency parse tree $\mathbf{z}$ (represented by a set of directed dependency edges) over a sequence of tokens $\mathbf{x}$ is measured by a score function $S(\mathbf{x}, \mathbf{z})$. One simple way of formulating the score function is to decompose the score by dependency edges in the tree, which is called first-order dependency factorization:
\[S(\mathbf{x}, \mathbf{z})=\sum_{(i,j) \in \mathbf{z}} s(\mathbf{x}, i, j)\]
where $s(\mathbf{x}, i, j)$ is defined as the score of the edge $(i,j)$.
For first order projective dependency parsing, \newcite{eisner-1996-coling} proposed an efficient dynamic programming algorithm for supervised parsing and \newcite{paskin2001cubic} proposed an inside-outside algorithm for unsupervised parsing. \newcite{mcdonald2005non} proposed a first order non-projective dependency parsing algorithm based on maximum spanning trees (MST) in directed graphs.

Our model draws inspiration from these approaches. However, instead of modeling the dependency relations between tokens of the input sequence, we model the relations between labels in the output label sequence.

\subsection{Latent Tree Models for Downstream Tasks}
Besides sequence labeling, latent tree models are useful in many other applications. \newcite{yogatama2016learning} proposed to learn task-specific constituency tree structures for natural language sentences with the REINFORCE algorithm \cite{williams1992simple}. Similarly, \newcite{choi2017unsupervised} proposed to learn task-specific constituency trees by employing the Gumbel-softmax trick \cite{jang2016categorical}. \newcite{williams2017learning} compared the two approaches empirically and found that these training methods lead to inconsistent trees of natural language sentences and the Gumbel-softmax trick leads to better classification accuracy than the REINFORCE algorithm.  There also exist previous work on latent dependency tree models that utilizes structured attention mechanisms \cite{kim2017structured,liu2017learning} for sentence and document classification problems.

\section{Neural Latent Dependency Model}
The neural latent dependency model can be seen as an extension of linear chain neural CRFs in which the dependency relations between the labels form a projective dependency tree instead of a linear chain.
The dependency tree has a dummy root node and we allow multiple edges connected to the dummy root. Given a label sequence with multiple tokens, there is more than one possible dependency tree and we regard the tree structure as a latent variable that shall be marginalized during inference. Figure \ref{fig:network} shows our model architecture with an example input sequence.

\begin{figure}[t]
	\centering
	\includegraphics[scale=0.2]{pic/network.pdf}
	\caption{\textbf{Left}: our model architecture for the input sequence \text{``OMG I love"} with the output label sequence \text{``!, O, V"}. \textbf{Right}: the set of 16 possible dependency trees, among which 12 are projective trees that are considered by our model.}
	\label{fig:network}
\end{figure}

\subsection{Model} \label{sec:model}
\paragraph{Input Representation} Each token in the input sequence is represented by a continuous embedding $\mathbf{x}_i \in \mathbb{R}^{D_x}$. A bi-directional recurrent neural network (RNN) \cite{rumelhart1986learning,hochreiter1997long} is applied on top of the embeddings to produce a vector for each token that captures context information of the token in the sequence. The bi-directional RNN produces two hidden vectors $\mathbf{h}_{lt}$ and $\mathbf{h}_{rt}$ for the left-to-right RNN and the right-to-left RNN seperately.

We concatenate the two hidden vectors from the two RNNs as our final representation for the $t$-th token $\mathbf{h}_t = [\mathbf{h}_{lt}; \mathbf{h}_{rt}] \in \mathbb{R}^{2D_{h}}$.

\paragraph{Edge Scoring} Given two positions $i$ and $j$ in a sequence $\mathbf{x}$ and its label sequence $\mathbf{y}$, we can calculate the score of a dependency edge from label $y_i$ to label $y_j$ as,
\begin{align*}
s(\mathbf{x}, \mathbf{y}, i, j) &= \phi (\mathbf{x}, y_j, j) + \psi(y_{i}, y_{j}) \\
&= \mathbf{h}_{j}^{T} \mathbf{W}_{e} \mathbf{t}_{y_{j}} + \bm{\psi}^{dir}_{y_i y_j}
\end{align*}
where $\phi (\mathbf{x}, y_j, j)$ is the emission score at position $j$, $\psi (y_{i}, y_{j})$ is the transition score of the dependency edge, $\mathbf{t}_{y_{j}}$ is the embedding of label $y_{j}$, $\mathbf{W}_{e}$ is a weight matrix, and $\bm{\psi}^{dir}$ is a direction-specific transition matrix chosen based on the direction of the dependency edge. Note that the dummy root node is at position 0 and has a fixed dummy label.

Instead of using a transition matrix $\bm{\psi}^{dir}$, we may also model the transition score with a bilinear function over label embeddings. Further, we may merge the emission and transition scores and compute the edge score with a trilinear function over the embeddings of the token and two labels, which is more expressive than the original score function,
\begin{align*}
    s(\mathbf{x}, \mathbf{y}, i, j) &= \mathbf{h}_{j}^{T} \mathbf{U} \mathbf{t}_{y_{i}} \mathbf{t}_{y_{j}}
\end{align*}
where $\mathbf{U}$ is a order-3 tensor. To reduce computational complexity, we assume that $\mathbf{U}$ has rank of $D_r$ and can be represented as the product of three matrices $\mathbf{U}_1\in \mathbb{R}^{D_r \times D_h}$, $\mathbf{U}_2$ and $\mathbf{U}_3$ $\in \mathbb{R}^{D_r \times D_l}$,
\begin{align*}
    & s(\mathbf{x}, \mathbf{y}, i, j) = \sum_{k=1}^{D_r} (\mathbf{v}_1 \circ \mathbf{v}_2 \circ \mathbf{v}_3)_k \\
    & \mathbf{v}_1 = \mathbf{U}_{1}\mathbf{h}_i;\quad\mathbf{v}_2 =  \mathbf{U}_{2}\mathbf{t}_{y_{i}};\quad\mathbf{v}_3 =  \mathbf{U}_{3}\mathbf{t}_{y_{j}};\quad
\end{align*}
where $\circ$ is the element-wise product.

\paragraph{Tree Scoring} Given an input sequence $\mathbf{x}$ and an output label sequence $\mathbf{y}$, the score of a dependency tree $\mathbf{z}$ over the labels (represented by a set of directed dependency edges) is factorized by edges:
\[s(\mathbf{x}, \mathbf{y}, \mathbf{z}) = \sum_{(i, j) \in \mathbf{z}} s(\mathbf{x}, \mathbf{y}, i, j)\]
This factorization is motivated by first order dependency parsing in natural language processing \cite{eisner-1996-coling,paskin2001cubic,mcdonald2005non}. We require that the dependency tree is projective, which means the edges in the tree cannot cross each other.

\paragraph{Label Sequence Scoring} In most applications, the dependency relations between output labels are not visible and therefore the dependency tree $\mathbf{z}$ can be seen as a latent variable. For a given input sequence $\mathbf{x}$ and an output label sequence $\mathbf{y}$, their score can be computed by marginalizing the latent variable $\mathbf{z}$:
\[S(\mathbf{x}, \mathbf{y}) = \sum_{\mathbf{z}} e^{s(\mathbf{x}, \mathbf{y}, \mathbf{z})} = \sum_{\mathbf{z}} e^{\sum_{(i,j) \in \mathbf{z}} s(\mathbf{x}, \mathbf{y}, i, j)}\]
or can be approximated by the maximum score over the latent variable $\mathbf{z}$.
\[S(\mathbf{x}, \mathbf{y}) \approx \max_{\mathbf{z}} e^{s(\mathbf{x}, \mathbf{y}, \mathbf{z})} = \max_{\mathbf{z}} e^{\sum_{(i,j) \in \mathbf{z}} s(\mathbf{x}, \mathbf{y}, i, j)}\]

The summation in the first equation can be efficiently computed using the inside algorithm \cite{baker1979trainable} and the maximization in the second equation can be efficiently computed using the CYK algorithm \cite{eisner-1996-coling}. Both algorithms have $O(n^3)$ time complexity and their running time can be significantly improved with vectorized operations in GPUs. Note that the original versions of the two algorithms output a dependency tree over the sequence of tokens, but here we run them to produce a tree over the sequence of labels.

\paragraph{Partition Function} We can also perform marginalization over $\mathbf{y}$ to compute the partition function,
 \[\mathbf{Z}(\mathbf{x}) = \sum_{\mathbf{y}} S(\mathbf{x}, \mathbf{y})=\sum_{\mathbf{y}} \sum_{\mathbf{z}} e^{s(\mathbf{x}, \mathbf{y}, \mathbf{z})}\]
We can efficiently compute the partition function using a modified inside algorithm. The basic idea behind the algorithm is that we regard our model as a weighted context-free grammar (WCFG) that contains two layers: the first layer is a projective dependency grammar (a subclass of WCFG) over the label sequence and the second layer is a lexicon that contains unary production rules in the form of a label producing a token. The inside algorithm on this WCFG computes the partition function of our model. The time complexity is $O(m^2 \times n^3)$ ($m$ is the number of possible labels for each token) and the running time can be improved with vectorized operations in GPUs.

There are two main reasons as to why we choose to model the dependencies between the output labels with a projective tree rather than a non-projective tree. One reason is that a projective tree contains more short dependencies which are desired in many sequence labeling problems. The other reason is that computing the partition function for non-projective trees is $\#$P-hard.

\subsection{Inference}
At test time, given an input token sequence $\mathbf{x}$, we aim to discover the output label sequence with the highest score,
\[\hat{\mathbf{y}}=\arg\max_{\mathbf{y}} S(\mathbf{x}, \mathbf{y})=\arg\max_{\mathbf{y}} \sum_{\mathbf{z}} e^{s(\mathbf{x}, \mathbf{y}, \mathbf{z})}\]
This problem is known as the Maximum A Posteriori problem \cite{murphy2012machine}, which is very difficult in general. Instead of trying to find an exact solution, we propose to solve this problem approximately,
\begin{align*}
\hat{\mathbf{y}}=\arg\max_{\mathbf{y}} S(\mathbf{x}, \mathbf{y})  &\approx \arg\max_{\mathbf{y}} \max_{\mathbf{z}} e^{s(\mathbf{x}, \mathbf{y}, \mathbf{z})} \\
&=\arg\max_{\mathbf{y}} \max_{\mathbf{z}} s(\mathbf{x}, \mathbf{y}, \mathbf{z})
\end{align*}
The last argmax in the equation can be exactly solved in polynomial time by simply replacing the \textbf{sum} operation with the \textbf{max} operation in our modified inside algorithm.

\subsection{Learning}
Given a training set of sequences $\mathbf{x}^1, \mathbf{x}^2, ..., \mathbf{x}^N$ and their label sequences $\mathbf{y}^1, \mathbf{y}^2, ..., \mathbf{y}^N$, the learning objective function is the regularized conditional log likelihood:
\begin{align*}
& J(\mathbf{\Theta}) = -\frac{1}{N}\sum_{\alpha = 1}^{N} \log P(\mathbf{y}^{\alpha} | \mathbf{x}^{\alpha}) + \Omega ||\mathbf{\Theta}||_2^2 \\
& =-\frac{1}{N}\sum_{\alpha = 1}^{N} \Big( \log S(\mathbf{x}^{\alpha}, \mathbf{y}^{\alpha}) - \log \mathbf{Z}(\mathbf{x}^{\alpha}) \Big)  + \Omega ||\mathbf{\Theta}||_2^2
\end{align*}
where $\mathbf{\Theta}$ contains the parameters of the LSTM and the edge score function, and $\Omega$ is a hyperparameter.
The gradient of $\log S(\mathbf{x}^{\alpha}, \mathbf{y}^{\alpha})$ can be computed via the higher-dimensional chain rule,
\begin{align*}
& \frac{\partial \log S(\mathbf{x}^{\alpha}, \mathbf{y}^{\alpha})}{\partial \mathbf{\Theta}} \\
& = \sum_{i} \sum_{j} \underbrace{\frac{\partial \log S(\mathbf{x}^{\alpha}, \mathbf{y}^{\alpha})}{\partial s(\mathbf{x}^{\alpha}, \mathbf{y}^{\alpha}, i, j)}}_{\text{term 1}} \underbrace{\frac{\partial s(\mathbf{x}^{\alpha}, \mathbf{y}^{\alpha}, i, j)}{\partial \mathbf{\Theta}}}_{\text{term 2}}
\end{align*}

It can be shown that term 1 is equal to ,
\begin{align*}
& \frac{\partial \log S(\mathbf{x}^{\alpha}, \mathbf{y}^{\alpha})}{\partial s(\mathbf{x}^{\alpha}, \mathbf{y}^{\alpha}, i, j)} \\
&= \frac{1}{S(\mathbf{x}^{\alpha}, \mathbf{y}^{\alpha})} \frac{\partial S(\mathbf{x}^{\alpha}, \mathbf{y}^{\alpha})}{\partial s(\mathbf{x}^{\alpha}, \mathbf{y}^{\alpha}, i, j)} \\
&=\frac{1}{S(\mathbf{x}^{\alpha}, \mathbf{y}^{\alpha})} \sum_{\mathbf{z}} e^{s(\mathbf{x}, \mathbf{y}, \mathbf{z})} \frac{\partial s(\mathbf{x}, \mathbf{y}, \mathbf{z}) }{\partial s(\mathbf{x}^{\alpha}, \mathbf{y}^{\alpha}, i, j)}\\
&= \sum_{\mathbf{z}} \frac{e^{s(\mathbf{x}, \mathbf{y}, \mathbf{z})}}{S(\mathbf{x}^{\alpha}, \mathbf{y}^{\alpha})}  \mathbbm{1}[(i, j) \in \mathbf{z} ]\\
&=  \sum_{\mathbf{z}} P(\mathbf{z}|\mathbf{x},\mathbf{y})  \mathbbm{1}[(i, j) \in \mathbf{z} ]\\
&=  \mathbbm{E}_{\mathbf{z} \sim P(\mathbf{z}|\mathbf{x},\mathbf{y})} \mathbbm{1}[(i, j) \in \mathbf{z} ]\\
\end{align*}

which can be exactly computed by the inside-outside algorithm \cite{baker1979trainable,eisner-1996-coling,paskin2001cubic}.
Term 2 can be computed using any automatic differentiation toolkit, such as those in TensorFlow and PyTorch.

In a similar manner, the gradient of the partition function can also be exactly computed using a modified inside-outside algorithm (by running the inside-outside algorithm on the WCFG introduced in section \ref{sec:model}) and automatic differentiation. Therefore, we can perform end-to-end parameter learning with gradient descent.

\subsection{Limiting Dependency Length}
The NLDM can capture dependency relations between labels of arbitrary distance. In most cases, however, very long dependencies are rare and therefore we can set an upper-bound  $k$ of the dependency length. The time complexity of both the inside-outside algorithm and the CYK algorithm can then be reduced to $O(m^2 \times n \times k^2)$.

\section{Relations to Previous Work on Sequence Labeling}

We show that the neural softmax model and the linear chain neural CRF can be seen as special cases of NLDMs, which is illustrated in Figure \ref{fig:comps}.

\begin{figure}[tb]
	\centering
	\includegraphics[scale=0.2]{pic/comp.pdf}
	\caption{\textbf{Top}: The neural softmax model shown as a special case of the NLDM (left) and the Markov network (right). \textbf{Bottom}: The linear chain CRF shown as a special case of the NLDM (left) and the Markov network (right).}
	\label{fig:comps}
\end{figure}

\paragraph{Relation to Neural Softmax}

In a NLDM, if we only allow dependencies from the root label $y_0$ to the other labels, then there is only one possible tree. Based on the transition-emission edge score formulation, the score of the $(\mathbf{x}, \mathbf{y})$ pair becomes,
\begin{align*}
S(\mathbf{x}, \mathbf{y}) &= \sum_{\mathbf{z}} e^{s(\mathbf{x}, \mathbf{y}, \mathbf{z})}
= e^{\sum_{i=1}^n s(\mathbf{x}, \mathbf{y}, 0, i)}
\\ &= e^{\sum_{i=1}^n \phi (\mathbf{x}, y_i, i) + \psi(y_{0}, y_{i})}
\end{align*}
If we further assume that $\psi(y_{0}, y_{i})$ is always 0, then we have
\begin{align*}
S(\mathbf{x}, \mathbf{y}) = e^{\sum_{i=1}^n \phi (\mathbf{x}, y_i, i)} = e^{\sum_{i=1}^n \mathbf{h}_{i}^{T} \mathbf{W}_{e} \mathbf{t}_{y_{i}}}
\end{align*}

Then the conditional probability of the labels given the input sequence becomes,
\begin{align*}
P(\mathbf{y}|\mathbf{x}) &= \frac{e^{\sum_{i=1}^n \phi(\mathbf{x}, y_i, i)}}{\sum_{\mathbf{y'}}e^{\sum_{i=1}^n \phi(\mathbf{x}, y'_i, i)}} \\ &= \frac{\prod_{i=1}^{n} e^{\phi(\mathbf{x}, y_i, i)}}{\prod_{i=1}^{n} \sum_{y'} e^{ \phi(\mathbf{x}, y', i)}} \\
&= \prod_{i=1}^{n} \frac{e^{\phi(\mathbf{x}, y_i, i)}}{\sum_{y'} e^{ \phi(\mathbf{x}, y', i)}}
\end{align*}
We can see that in this case the NLDM becomes equivalent to a neural softmax model described in Section \ref{sec:2.2}.

\paragraph{Relation to Linear Chain Neural CRFs}

In a NLDM, if we only allow dependencies from every label to its next label, then again there exists only one possible tree. Based on the transition-emission edge score formulation, the score of the $(\mathbf{x}, \mathbf{y})$ pair becomes,
\begin{align*}
    &S(\mathbf{x}, \mathbf{y})= \sum_{\mathbf{z}} e^{s(\mathbf{x}, \mathbf{y}, \mathbf{z})
    = e^{\sum_{i=1}^n s(\mathbf{x}, \mathbf{y}, i-1, i)} \\
    &= e^{\sum_{i=1}^n \left( \phi (\mathbf{x}, y_i, i) + \psi(y_{i-1}, y_{i}) \right)}

\end{align*}
which is exactly the linear chain neural CRF model illustrated in Section \ref{sec:2.1}.

\section{Experiments}
We evaluate our model on both synthetic and real data.
For comparison, we choose the following baseline models for sequence labeling.

\paragraph{Neural Softmax} Neural Softmax is a simple and strong baseline in sequence labeling. It considers labeling sequence as several independent supervised classification problems. In our experiments, we used the same feature extraction layer for comparison.

\paragraph{Neural CRFs} Neural CRFs are important baselines since our model can be seen as their extension. We use 1/2-order CRFs as our two CRFs baselines. In our experiments, we used the same feature extraction layer for neural CRFs and NLDMs in order to make a fair comparison. We also perform the trilinear function for neural CRFs models.

We use a Bi-LSTM as the neural encoder in all the models. The code of our approach and the baseline systems, the hyper-parameters, and the trained models will be released at \url{https://anonymous.link}.

\subsection{Hyper-parameters}
For all the models, we tune the learning rate in the range of \{0.03, 0.1, 0.3\} and the LSTM hidden size $D_h$ in the range of \{100, 200, 400\}. For CRF and NLDM models, we select the label embedding size $D_t$ from \{10, 50, 100\}, tensor rank $D_r$ from \{100, 400, 600\}, and the dependency length limit $k$ from \{2, 3, \dots, 10\}.

\subsection{Synthetic data}
We compared the NLDM and the neural CRF on a synthetic dataset generated from a novel infinite order hidden Markov model in which, at each step, the transition distribution is computed by an LSTM that takes the existing label sequence as input, and the emission distribution is computed by another LSTM that takes the existing word sequence as input.
To generate the dataset, we set the label number to 5, word dictionary size to 1000, and maximum sequence length to 10. For the LSTMs, the size of both the hidden state and the output state is set to 50 while the parameters are randomly generated. The embeddings of the labels and words are also randomly generated.  We generated 1k, 3k, 10k, 30k, and 100k samples and randomly selected 80\% for training, 10\% for validation and 10\% for testing.
For all the experiments, we did not tune the hyper-parameters and kept them the same in the neural CRF and NLDM.

\begin{figure}[t]
\centering
\includegraphics[scale=0.5]{pic/syn-results.pdf}
\caption{Results on the synthetic data. The neural CRF and the NLDM utilize the Bi-LSTMs as the neural encoder.}
\label{fig:syn-results}
\end{figure}

The experimental results are shown in Figure \ref{fig:syn-results}. As can be seen, with more training samples, the performances of both the neural CRF and the NLDM increase. When the sample size is small, the NLDM has significant advantage over the neural CRF model. For larger sample sizes, the NLDM still achieves better accuracies although with smaller margins.
\begin{table*}
\setlength\tabcolsep{4.3pt}
\renewcommand\arraystretch{1.16}
	\centering

	\begin{tabular}{cccccccccc}
		\hline
		{\bf } & {\bf\textsc{Tw}} &{\bf\textsc{En}}& {\bf\textsc{Fr}}& {\bf\textsc{It}}& {\bf\textsc{Te}}& {\bf\textsc{Pl}}& {\bf\textsc{Mt}}& {\bf\textsc{Ru}}& {\bf\textsc{Sv}} \\\hline
        \# Vocab  &4806 & 94217& 35373 & 124445&6465&130967&44162&98000&90960\\
		\# Labels  & 25 &  17&  15 & 16&14&15&17&16&17\\
		\# Training & 1327 & 3176 &  1153 & 5368&1051&13774&1123&3850&2738\\
        \# Dev & 500 & 1032& 907& 671&131&1745&433&579&912\\
        \# Test & 547 & 1035&726& 674&146&1727&518&601&914\\ \hline
	\end{tabular}
	\caption{Statistics of the nine sequence labeling datasets.}
    \label{table:pos-data}
\end{table*}

\begin{table*}
    \centering
    \setlength\tabcolsep{5.3pt}
\renewcommand\arraystretch{1.16}

    \begin{tabular}{ccccccccccc}
    \hline
              & {\bf\textsc{Tw}} & {\bf\textsc{En}}& {\bf\textsc{Fr}}& {\bf\textsc{It}}& {\bf\textsc{Te}}& {\bf\textsc{Pl}}& {\bf\textsc{Mt}}& {\bf\textsc{Ru}}& {\bf\textsc{Sv}}& {\bf\textsc{Avg.}} \\ \hline
    Neural Softmax &  88.73&90.22 & 89.24 & 89.70 & 82.79 & 89.02 & \bf89.16 & 82.78 & 89.50&87.90\\
    Neural 1-order CRF & 89.09 & 90.22 & 89.24 & 89.71 & 83.36 & 89.01 & 88.89 & 82.79 & 89.39&87.97\\
    Neural 2-order CRF &   89.29   &90.18 & 89.45 & 89.71 & 85.35 & \bf89.22 & 89.13 & 82.95 & 89.48   &   88.31   \\
    NLDM & {\bf 89.37} &{\bf90.41} & \bf89.64 & \bf90.58 & \bf86.45 & 89.18 & 89.11 & \bf83.23 & \bf89.77   &  \bf88.64    \\ \hline
    \end{tabular}
\caption{POS tagging accuracy.}
\label{table:twitter-result}
\end{table*}

\subsection{Real data}\label{sec:real-data}
We performed experiments on the POS tagging task. We used the Twitter datasets from \newcite{gimpel2011part} and \newcite{owoputi2013improved}. We merged the 1000-tweet OCT27-TRAIN set and the 327-tweet OCT27-DEV set as our training dataset, used the 500-tweet OCT27-TEST set for validation, and used the 547-tweet DAILY547 test set for testing. We also selected eight languages from Universal Dependencies tree banks (UD). The statistics of the datasets are showed in Table \ref{table:pos-data}.
We report the experimental results in Table \ref{table:twitter-result}. It can be seen that our model achieves the best results in most of the experiments.

\begin{table}
    \centering

\begin{tabular}{ccc}
    \hline {\bf\textsc{Dependency Len}} & {\bf\textsc{Freq.}} & {\bf\textsc{Freq. (gold)}}\\ \hline
    1 & 	65.3 & 65.4\\
    2-10 & 	20.7 & 34.3\\
    $>$ 10 & 	14.0 & 0.3\\ \hline
\end{tabular}
\captionof{table}{Frequency (\%) of different dependency lengths in parses decoded by our model and in gold parses on the Twitter dataset.}
\label{table:freq-edge-tree}
\end{table}

\section{Analysis}
\subsection{Impact of Dependency Length Limit}
We analyze the performance of our model on the Twitter POS tagging dataset with different values of the dependency length limit $k$. We repeat the experiments in Section \ref{sec:real-data} with $k\in\{1,2,5,10,15\}$ and report the results in Figure \ref{fig:tree-acc}.
Note that $k=1$ corresponds to linear chain neural 1-order CRFs. It can be seen that allowing long-range dependencies is indeed beneficial compare with the linear chain baseline. However, setting the $k$ value beyond 5 seems to decrease the accuracy, possibly because dependencies longer than 5 are relatively rare on the dataset.

\subsection{Analysis of Induced Tree Structures}
We use the learned NLDM from Section \ref{sec:real-data} as a dependency parser by solving $\mathbf{z}^*=\arg\max_{\mathbf{z}} \max_{\mathbf{y}} s(\mathbf{x}, \mathbf{y}, \mathbf{z})$ and compare the dependency trees decoded by the NLDM with the gold linguistic trees on the Twitter dataset.

\paragraph{Comparison on Dependency Length}
We first aim to analyze whether the NLDM indeed learns to capture long-range dependency relations. We count the frequencies of different dependency lengths in the decoded dependency parses and report the results in Table \ref{table:freq-edge-tree}. From the table we can see that over 30 percent of the dependencies are longer than 1, suggesting that our model is able to capture long-range dependencies. In addition, the NLDM can be seen to induce longer dependencies compared with the gold parses.

\paragraph{Comparison with Gold Trees}
We compare the dependency trees decoded from our model with the gold linguistic dependency trees from the Tweebank \cite{kong-EtAl:2014:EMNLP2014}. The unlabeled dependency accuracy of our model is 17.9\%, which suggests that our model induces dependency trees that are significantly different from the gold trees, echoing the findings from previous work \cite{yogatama2016learning,williams2017learning}.

\begin{figure}
\centering
\includegraphics[scale=0.45]{pic/analysis.pdf}
\captionof{figure}{Accuracy on the Twitter test set with different dependency length limits.}
\label{fig:tree-acc}
\end{figure}

\section{Conclusion}
We propose a novel latent variable probabilistic model for sequence labeling called the Neural Latent Dependency Model (NLDM). A NLDM allows dependencies between labels of arbitrary length and models the set of dependencies with a latent projective tree structure. We prove that the neural softmax model and the linear chain neural CRF model are two special cases of our model. We present an end-to-end training algorithm and an efficient approximate decoding algorithm. We tested our model on both synthetic and real datasets and show that our model achieves superior performance over baselines. For future work, we plan to design better decoding algorithms and apply our model to other more complicated sequence labeling tasks such as semantic role labeling.

\section*{Acknowledgement}
This work was supported by Alibaba Group through Alibaba Innovative Research Program. This work was also supported by the National Natural Science Foundation of China (61976139).


\section{Introduction}
\label{sec:intro}
The past few years have seen the advent of high-energy neutrino astronomy, with the identification and characterization of a diffuse astrophysical flux of neutrinos at TeV energies and above~\citep{Aartsen:2016xlq, Aartsen:2013jdh,Aartsen:2014gkd}, as well as the first observation of an extragalactic source of such high energy neutrinos~\citep{IceCube:2018dnn, IceCube:2018cha}. Neutrinos nevertheless remain some of the most underutilized of the astronomical particle messengers. Their production is theorized in a wide range of astrophysical scenarios, and may be the only observable signatures of processes in dense regions of the Universe. The IceCube Neutrino Observatory (IceCube) is composed of a surface array as well as an in-ice array consisting of $\sim$1 cubic kilometer of instrumented ice located at the South Pole~\citep{Aartsen:2016nxy}. Most IceCube searches for astrophysical neutrino emission are optimized for TeV-PeV energies, relying on the enhanced effective volume (due to the increased cross sections as well as the longer propagation distances of the secondary particles) of the detector at these high energies and the steeply falling spectrum of atmospheric neutrino and muon backgrounds. In the context of considering a `stacking' analysis of sources, a 10-year data study by IceCube finds evidence for neutrino emission in a set of 4 objects, including TXS0506+056~\citep{Aartsen:2019fau}.

Searches for astrophysical sources of neutrinos at sub-TeV energies have been carried out with Super-Kamiokande~\citep{Thrane:2009tw, Abe:2018mic}, as well as with a single year of IceCube data~\citep{Aartsen:2015eai}.  The analyses have been predominantly focused on using only $\nu_\mu$, and the untriggered `tabula rasa’ search in Super-Kamiokande tested for power-law spectra. All searches suffer from significantly larger background event rates compared to the high-energy searches due to the steeper power laws that characterize atmospheric background spectra. However, they remain a promising avenue for the discovery of transient emitters, which are in general not expected to have power law spectra~(\cite{longair_1994}, see chapter 21). Scenarios that predict flaring neutrino emission which peak in the 10-100\,GeV range along with a suppression of emission at higher energies are of particular interest. This may occur due to the subphotospheric emission of neutrinos from gamma-ray bursts (GRBs), wherein the jet is dissipated by inelastic collisions and quasithermal neutrinos are produced via hadronuclear reactions~\citep{Murase:2013hh}. Detecting sub-TeV neutrinos from GRBs can thus delineate the subphotospheric scenario from the more classical case, in which gamma-rays are attributed to synchrotron emission from nonthermal electrons accelerated at internal shocks~\citep{Toma:2010xw}. Neutrino emission that peaks at $\sim$GeV energies is also naturally expected in core collapse supernovae in which the jets are not energetic enough to break through the stellar envelope and gamma-ray emission is thus `choked' off~\citep{Ando:2005xi}. The observed correlation between long duration GRBs and a small fraction of core-collapse supernovae~\citep{Modjaz:2011bm} hints that a larger number of core collapse supernovae may develop such choked jets. It has been suggested~\citep{Taboada:2010xf} that neutrino emission from such events occurring  within $\sim$10\,Mpc should be detectable by the IceCube-DeepCore sub-array.

This work presents the search for transient emitters of neutrinos at $\sim$GeV energies with three years of IceCube data and marks the first use of an all-flavor sample of neutrinos at these low energies for such purposes. While the inferior angular resolution of, and limited effective volume for, $\nu_e$ (+ $\bar{\nu}_e$) and $\nu_{\tau}$ (+ $\bar{\nu}_{\tau}$) events at high energies (TeV) discourage their use in high-energy searches~\citep{Aartsen:2019fau}, in low-energy searches targeting transient emission, they form an additional signal component that enhances the potential for discovery.

In this work, the analysis has been performed on archival data, however, the possibility of employing this as an online real-time search (alert system) is being investigated.

\section{Detector and event sample}
\label{sec:DetectorEventSample}
The IceCube Neutrino Observatory consists of roughly 1 cubic kilometer of instrumented glacial ice located at the South Pole~\citep{Aartsen:2016nxy}. The detector contains 86 vertical strings each composed of 60 digital optical modules (DOMs) including a downward facing 10-inch photomultiplier tube. The DOMs are deployed from 1.45\,km to 2.45\,km below the surface of the ice mostly in a hexagonal pattern. The largest part of the detector, the in-ice array, is optimized for neutrinos at $\sim$ TeV -- PeV energies and has some sensitivity to the detection and analysis of neutrinos in the $\mathcal{O}(100)$\,GeV region. IceCube-DeepCore is the more densely instrumented region located approximately in the middle of the detector~\citep{Collaboration:2011ym}. It consists of the most central standard IceCube strings as well as 8 specialized low-energy strings whose DOMs have a higher quantum efficiency than the regular IceCube DOMs. Hence, this part of the detector is sensitive to neutrinos at energies down to $\sim$5\,GeV.

\begin{table}[b]
\centering
\caption{Relative fractions of event types for the GRECO dataset for each neutrino flavor (assuming an atmospheric flux) interacting through charged current (CC), all-flavor neutrinos interacting through neutral current (NC), atmospheric muons, and accidental events caused purely by detector noise \citep{Aartsen:2019tjl}. This is the event distribution used for background in our search.}
\label{tab:GRECO_ratio}
\begin{tabular}{ c c }
\hline \hline
 \bf{Event Type} & \bf{Relative Fraction (\%)}  \\  \hline
 $\nu_\mu + \bar{\nu}_\mu$ (CC) &  57.40  \\ \hline
 $\nu_e + \bar{\nu}_e$ (CC) &  21.64 \\ \hline
 $\nu_\tau + \bar{\nu}_\tau$ (CC) & 2.90  \\ \hline
 $\nu$ (NC) & 9.83  \\ \hline
 $\mu$ & 8.07  \\ \hline
 pure-noise & 0.01\\ \hline \hline
\end{tabular}
\end{table}

In this work, we use the GRECO (GeV Reconstructed Events with Containment for Oscillation) event selection. This selection was originally developed for the tau neutrino appearance analysis which is labelled as `analysis $\mathcal{A}$' in previous work~\citep{Aartsen:2019tjl}. The three years of data are obtained from April 2012 through May 2015 and have a livetime of $\sim$2.756 years (1006 days) and an average rate of 0.87\,mHz. The effective area for the event sample is shown in Figure~\ref{fig: Effective_area}, where also the individual contributions for the different neutrino flavors are indicated. This is compared to the effective area for the previous low-energy transient point-source search performed with IceCube data \citep{Aartsen:2015eai}. The background to any potential transient neutrino signal is a combination of atmospheric neutrinos, atmospheric muons, and PMT/DOM glass/electronics noise events that manage to trigger the array and pass the event selection. For the GRECO data the relative contributions of each of the major background types are shown in Table~\ref{tab:GRECO_ratio}.

\begin{figure}
	\center
	\includegraphics[width=0.7\linewidth]{effective_area_paper_largertext.pdf}
	\caption{All sky (4$\pi$\,sr) average neutrino effective area for the GRECO event selection. The effective area of the previous low-energy transient analysis from IceCube \citep{Aartsen:2015eai}, which was only $\nu_\mu$, is indicated by the purple dashed line.}
	\label{fig: Effective_area}
\end{figure}

\section{Analysis method}
\label{sec:AnalysisMetod}
This analysis consists of two parts: a kernel density estimation (KDE) used on the event times with the purpose of selecting time windows with event densities that are larger than a threshold defined to exclude background. This is followed by a maximum likelihood analysis that tests for spatial clustering of the events observed close together in time.

While establishing the performance of the analysis, time scrambled data is used for the background estimation. This means that we use experimental data obtained with IceCube-DeepCore but assign random observation times within the run-time of the detector for the individual events. Using scrambled data ensures that the algorithms and parameters values are developed and selected without knowledge of the final analysis results.

\subsection{Kernel density estimation}
\label{sec:KDE}

\begin{figure}
	\center
	\includegraphics[width=0.7\linewidth]{KDEvsMS100siginj1day.pdf}
	\caption{Example of time KDE. The x-axis is unit fraction of one day. Gaussian kernels for individual events are summed in order to compute the density of events with respect to time. No normalization is performed as this would change the scale of the y-axis for different numbers of events. }
	\label{fig: time KDE}
\end{figure}

The purpose of the time KDE is to identify the regions in time with large densities of events and thereby select the time windows of interest. Transient sources are expected to emit neutrinos primarily within a relatively short amount of time, and therefore the selected time windows with larger event densities are more likely to contain astrophysical transient emission. This way, the time KDE functions as a preselection of the most promising events. Thereby the maximum likelihood analysis is only performed on a small sub-set of the data sample, which will decrease the computing time.

When evaluating the KDE, a Gaussian kernel defined by
\begin{linenomath}
\begin{equation} \label{eq: kernel}
    K(t,t_n,\sigma) = \frac{1}{\sqrt{2\pi}\sigma} \exp \left(-\frac{|t-t_n|^2}{2\sigma^2} \right),
\end{equation}
\end{linenomath}
is created around the time of each event with a bandwidth $\sigma$ that should correspond to the neutrino emission time of the astrophysical transient source class under consideration. The value of $t-t_n$ is the difference between the kernel center and a given event. Since the width of the time distribution of neutrinos emitted from gamma-ray bursts has yet to be determined, we assume a bandwidth of 100\,s corresponding approximately to the gamma-ray emission time for long duration gamma-ray bursts~\citep{Harmon:2004wm}.

\begin{figure*}[t]
	\center
	\includegraphics[width=0.49\linewidth]{Compilation_track_final.pdf}
	\includegraphics[width=0.49\linewidth]{Compilation_cascade_final.pdf}
	\caption{Median angular resolutions, i.e., opening angles between true and reconstructed event directions, used for event angular resolution splines for three example declinations. Left panel shows the result for tracks, while the right panel shows the result for cascades. The error bars indicate the 68\% central intervals and are located at the energy bin centers. The asymmetry in the error bars illustrates that there is a `longer tail' of larger angular resolutions than smaller resolutions. Orange and purple error bars have been shifted slightly towards smaller and larger values, respectively, in order to avoid overlapping. }
	\label{fig: spline}
\end{figure*}

The sum of the individual contributions from the kernels is evaluated by
\begin{linenomath}
\begin{equation} \label{eq: R_KDE}
    R_{KDE}(t,\sigma) = \sum _{n=1}^{N} K(t,t_n,\sigma),
\end{equation}
\end{linenomath}
where $N$ is the total number of events (kernels). All maxima in time, denoted as {\sl cluster centers}, are found, and the maxima above a chosen threshold on the value of $R_{KDE}$ are selected. A time window with a potential transient (TWPT) of $\pm$3 bandwidths (600\,s) is created around each cluster center passing the threshold. The threshold is chosen to allow a given average number of cluster centers per year to pass. One time window is created around each, resulting in an expected average number of time windows per year (this parameter will be denoted as $N_{TWPT}$). In this case, an $R_{KDE}$-threshold of 905.80 is chosen corresponding to $N_{TWPT}=$100. A relatively small number of clusters passing the threshold is chosen because a larger number provides no appreciable improvements in the sensitivity. Figure~\ref{fig: time KDE} shows an example for one day of scrambled data.

\subsection{Likelihood Maximization}
\label{sec:MaximumLikelihood}

A maximum likelihood analysis is performed on the events within each selected time window created based on the time KDE. The likelihood is given by
\begin{linenomath}
\begin{equation}
    \mathcal{L}(n_s) = \frac{(n_s + n_b)^N}{N!}e^{-(n_s + n_b)} \prod_{i=1}^N \left( \frac{n_s \mathcal{S}_i}{n_s + n_b} + \frac{n_b \mathcal{B}_i}{n_s + n_b} \right),
\end{equation}
\end{linenomath}
where the first two factors describe the Poisson distributed probability of observing $N$ events from the source under consideration.
The log-likelihood ratio is then given by
\begin{linenomath}
\begin{equation} \label{eq: log-likelihood ratio}
    \ln \left( \frac{\mathcal{L}(n_s)}{\mathcal{L}(0)} \right) = -n_s + \sum_{i=1}^{N} \ln \left( \frac{n_s \mathcal{S}_i}{\left< n_b \right> \mathcal{B}_i} + 1 \right).
\end{equation}
\end{linenomath}
The parameter $n_s$ and $n_b$ are the numbers of signal and background events within the given time window, while $\left< n_b \right>$ is the expected average based on the average background rate and the size of the time window. The parameter $n_s$ refers to the number of signal events being fit during the maximization. The variables $\mathcal{S}_i$ and $\mathcal{B}_i$ are the signal and background probability distribution functions (PDFs), i.e, the probabilities that the event $i$ within a time window originates from the source and the background, respectively.

The background PDF $\mathcal{B}_i$ is assumed to be constant in azimuth, but changes as a function of the zenith angle of the event $i$. The value is calculated by
\begin{linenomath}
\begin{equation}
    \mathcal{B}_i = \mathcal{P}_{i, zenith} \cdot \mathcal{P}_{i, azimuth} = \mathcal{P}_{i, zenith} \cdot \frac{1}{2\pi} ,
\end{equation}
\end{linenomath}
where $\mathcal{P}_{i, zenith}$ is found using a spline with 25 bins in the cos(zenith) distribution of events, i.e., the background PDF roughly follows the zenith distribution of atmospheric neutrinos with a small contribution from atmospheric muons, and is based on experimental data.

The signal PDF $\mathcal{S}_i$ contains spatial information about the potential sources as well as the events. The expected angular resolutions of the events are estimated from their reconstructed energies and declinations using a spline of the median angular resolution found from simulated GRECO events. In this case, the resolution refers to the opening angle between the true and reconstructed event directions, and the data have been weighted according to the expected energy spectrum for atmospheric neutrinos described by~\citep{Honda:2015fha}. Separate splines are employed for each event type (tracks and cascades) in order to account for the better angular resolutions of the tracks. In this case, tracks are defined as events with reconstructed track lengths of at least 50\,m. The splines are functions of energy and declination for the events. Bins in declination are distributed linearly in the range -85$^\circ$ to 85$^\circ$, with bin centers 5$^\circ$ apart. Bins in energy are logarithmically distributed with an additional bin at higher and/or lower energies (see Figure~\ref{fig: spline}).

The signal PDF $\mathcal{S}_i$ is defined as a Kent distribution, i.e., a probability distribution analogous to the bivariate normal distribution but normalized correctly on a sphere. Assuming circular angular errors on the events, i.e, angular resolutions in azimuth and zenith are equal, the Kent distribution takes the form
\begin{linenomath}
\begin{equation} \label{eq: signal PDF}
    \mathcal{S}_i = \frac{\kappa_i}{4\pi \sinh(\kappa_i)} \exp \left( \kappa_i \cos(|x_{\mathrm{source}}-x_i|) \right).
\end{equation}
\end{linenomath}
Here, $|x_{\mathrm{source}}-x_i|$, is the angular distance between the source and the event $i$. The parameter $\kappa_i$ is the concentration parameter for the distribution, related to the angular resolution of the events by
\begin{linenomath}
\begin{equation} \label{eq: kappa}
    \kappa_i \approx \frac{1}{\sigma_i^2},
\end{equation}
\end{linenomath}
where the median of the differences between true and reconstructed event directions are used as the angular resolution. This will be discussed more thoroughly below. The angular extent of the source is omitted since transient emitters are expected to be point sources.

In addition to the $R_{KDE}$-threshold, a threshold is applied to the value of signal events $n_s$ in order to require a minimal amount of spatial clustering (spatial correlation) between the events within a time window. In this work, the $n_s$-threshold is set to 2.
The log-likelihood ratio defined in Equation~\ref{eq: log-likelihood ratio} is maximized with respect to the value of $n_s$ and the source direction, $x_{source}$, included in the Kent distribution. Several distinct directional maxima may exist, and therefore we perform the maximization once for each event within the time window, while using the best-fit event position as the initial guess for the position of a possible transient point source. This ensures that all relevant maxima are found, and subsequently, the largest maximum passing the $n_s$-threshold is selected.

\subsection{Performance}
\label{sec:Performance}

\begin{figure*}
	\center
	\includegraphics[width=0.49\linewidth]{Source_uncertainty_SPS.pdf}
	\includegraphics[width=0.49\linewidth]{Source_uncertainty_dec.pdf}
	\caption{Uncertainty of a simulated source position as function of the number of injected neutrinos, $n_s$. Left panel: The uncertainty decreases with the neutrino mean energy for a fixed declination of 16 degrees. Right panel: Same plot but varying the declination for a fixed mean energy of 100\,GeV.}
	\label{fig: source_uncertainty}
\end{figure*}

\iffalse
\begin{figure*}
	\center
	\includegraphics[width=0.48\linewidth]{Tw_SPS20.pdf}
	\includegraphics[width=0.48\linewidth]{Sen_ns_threshold.pdf}
	\caption{The 90\% sensitivity for a source emitting neutrinos with a mean energy of $E_{\mathrm{mean}}=20$\,GeV located at declination of 16 degrees. The number of signal events on the y-axis refers to the average number of injected signal events. Left: The $n_s$-threshold is fixed at 2, while the KDE-threshold/average number of time windows (N$_{TWPT}$) is varied. Right: Average number of time windows is fixed at 100 per year and the $n_s$-threshold is varied. }
	\label{fig: sen_thresholds}
\end{figure*}
\fi

The performance of the analysis is estimated with simulations by injecting simulated signal events into the scrambled background data. Thereby the uncertainty on the likely source position is found when maximizing the natural log of the likelihood ratio as well as the sensitivity and discovery potential for the analysis.

For the injection of signal events into the background data, the signal neutrinos are weighted according to a Dirac spectrum that mimics the shape expected from an astrophysical transient~\citep{Murase:2013hh}. The Dirac spectrum is parameterized as
\begin{linenomath}
\begin{equation}
    \Phi(E) = N \cdot \frac{\left( \frac{E}{kT} \right)^2}{e^\frac{E}{kT} + 1},
\end{equation}\label{eq: Dirac_fluence}
\end{linenomath}
where $kT$ depends on the the mean energy of the neutrinos
\begin{linenomath}
\begin{equation}
    kT = \frac{E_{\mathrm{mean}}}{3.15}.
\end{equation}
\end{linenomath}
In this work, we use the mean energies of ${E_{\mathrm{mean}}=\mathrm{20}}$\,GeV and $E_{\mathrm{mean}}=$100\,GeV as benchmarks.

\begin{figure*}
	\center
	\includegraphics[width=0.49\linewidth]{Even_sen_dis.pdf}
	\includegraphics[width=0.49\linewidth]{fluence_sen_dis.pdf}
	\caption{Left panel: Event sensitivity and discovery potential for the two different neutrino emission energies. The number of signal events on the y-axis refers to the average number of injected signal events. Right panel: The event sensitivity and discovery potential converted to fluences using the effective area for the GRECO event selection. The peak of the fluence is used as reference point for the energy (See also Section \ref{sec:UnblindResults}).}
	\label{fig: Sen_dis}
\end{figure*}

The uncertainty on the source position is defined as the difference between the true injected source position and the likely position found when performing the likelihood maximization. In this case, only the time windows actually containing the source are considered. The uncertainty is estimated and the result is shown in Figure~\ref{fig: source_uncertainty}.
The left panel shows that the uncertainty decreases with the number of signal events $n_s$ and with the neutrino energy. The source position is generally estimated more accurately when more neutrinos from the source are detected. Higher energy neutrinos are reconstructed better, leading to the events being observed closer to the true source position. This pattern is also seen in Figure~\ref{fig: spline} for the angular resolution splines. The right panel of Figure~\ref{fig: source_uncertainty} shows that the uncertainty on the source position decreases with larger (more positive) values of the source declination, which is also a result of the event angular resolution being better at larger declinations. The angular resolutions strongly affects the signal PDF since it is being squared in Equation~\ref{eq: kappa} and because $\kappa$ is included in the exponent in equation~\ref{eq: signal PDF}. The central median values of the angular and energy resolution are used in this transient analysis for each candidate neutrino within the examined time window(s). While the 68\% confidence intervals for the resolutions are not explicitly used, they reflect the spread of neutrino angular resolutions that are used for the pseudo-trials necessary to establish the sensitivity to an astrophysical transient source.

In this work, we define the test statistic (TS) as the largest value of the maximized natural log of the likelihood ratios for either observed data or a background-only hypothesis. Owing to the fact that Wilks' theorem is not satisfied for this analysis we are unable to define an analytical expression for the TS distributions. Therefore, the background TS distribution used to establish the sensitivity and discovery potential for this analysis is generated from scrambled trials, where a scrambled trial corresponds to a year of time scrambled data. This motivates an estimation of a discovery potential at a 3$\sigma$ significance rather than at the more traditional 5$\sigma$ due to computational challenges. Additionally, a 3$\sigma$ discovery potential is further motivated by the possibility of using this analysis for sending transient astrophysical alerts~\citep{Smith:2012eu,Rutledge:1998hc,Aartsen:2016lmt}, in which case 3$\sigma$ is a reasonable threshold.

Figure~\ref{fig: Sen_dis} shows the 90\% sensitivity and 3$\sigma$ discovery potential for the analysis in terms of the number of observed signal events as well as fluence. The 90\% sensitivity is defined as the signal strength that results in 90\% of the scrambled trials containing signal having a value of the TS equal to or larger than the median background TS. The 3$\sigma$ discovery potential requires the median of the scrambled trials containing signal to coincide with the 3$\sigma$ background TS value. The sensitivity in terms of the number of signal events is slightly worse around the horizon, since the background rate at those declinations is higher. The sensitivity in terms of the fluence is better for larger declinations due to the effective area of the detector being better in the Northern Sky.

For a source emitting neutrinos with a mean energy of $E_{mean}$ \nolinebreak= 20 GeV located at a declination of 16$^{\circ}$, it was found that increasing the number of time windows considered (when $n_s$ threshold is fixed to 2) does not significantly improve sensitivity (which plateaus at $\sim$8 signal events). However, the computation time is significantly increased, motivating us to consider not more than 100 windows on average per trial. The sensitivity when the average number of time windows considered is fixed to 100, is also not significantly affected by the $n_s$ threshold (the sensitivity improves from 7.8 to 7.6 signal events when the $n_s$ threshold is varied from 1.5 to 4.0), which is therefore kept fixed at 2 in order to require a minimum of spatial clustering between the events.

The background rate is not constant as a function of the azimuth angle due to the detector geometry having slightly different atmospheric muon rejection efficiency and neutrino acceptance, and is also not constant in time as a result of seasonal temperature variations in the atmosphere. In order to account for the fluctuations in time and azimuth, the discovery potential is estimated with an increased background rate of 15\%. The result of an increase in the background is shown in Figure~\ref{fig: Sen_dis_15}, where the discovery potential is worsened by less than 5\%.

\begin{figure}[t]
	\center
	\includegraphics[width=0.7\linewidth]{Sen_dis_15.pdf}
	\caption{Comparison between the discovery potential using the actual event rate and using a rate increased by 15\%. The number of signal events on the y-axis refers to the average number of injected signal events. The red dots in the ratio plot show the ratio between the discovery potential for a source with a mean energy of 100\,GeV with a normal background rate and with an increased rate. Similarly, the blue points in the ratio plot corresponds to the ratio between the discovery potentials with normal rate and increased rate for a source with a mean energy of 20\,GeV.}
	\label{fig: Sen_dis_15}
\end{figure}

\section{Results}
\label{sec:ResultsConlusion}
We used blind analysis techniques for GRECO data collected between April 2012 and May 2015. The unblinded results are presented in this section.

\subsection{Unscrambled results}
\label{sec:UnblindResults}

When unscrambling the three years of GRECO data we find 300 cluster centers above the $R_{KDE}$-threshold resulting in 300 time windows. After applying the $n_s$-threshold, 267 time windows remain.
The comparison between the TS for the best-fit signal and source position and the TS distribution for three years of scrambled background data is shown in Figure~\ref{fig: Result}. The 3$\sigma$ TS value for the background distribution is indicated by the black dashed line. The best-fit TS value is smaller than the 3$\sigma$ TS for background, and hence, this  analysis reveals no transient neutrino emission from point sources in the data.

\begin{figure}[t]
	\center
	\includegraphics[width=0.7\linewidth]{Unblinding_result.pdf}
	\caption{TS for the best-fit signal and source position in the three years of GRECO data (2012-2015) compared to the background TS distribution based on one year of data. }
	\label{fig: Result}
\end{figure}

Since no transient neutrino emission from point sources was found in the data, we place upper limits on the volumetric rate ($\dot{\rho}$) of the transient phenomena in the local universe.
The discovery potential and sensitivity are given in terms of peak fluence, which is given by
\begin{linenomath}
\begin{equation} \label{eq: Peak fluence}
E^2 \Phi_\nu = \int_0^\infty \mathrm{d}z \frac{\mathcal{E}}{4\pi d_L^2(z)} \frac{\mathrm{d}N}{\mathrm{d}z} e^{-N(z)},
\end{equation}
\end{linenomath}
where the parameter $\mathcal{E}$ is the bolometric energy of the source and $N(z)$ describes the number of expected sources within redshift $z$ and is parameterized as
\begin{linenomath}
\begin{equation} \label{eq: Number_src}
N(z) = \dot{\rho} T f_{sky} \frac{4\pi}{3} d^3_c(z),
\end{equation}
\end{linenomath}
with $T$ being the livetime of the data used and $f_{sky}$ denotes the fraction of the sky considered. Since this work is a full sky search,  $f_{sky}=$1.
The factor
$\frac{\mathrm{d}N}{\mathrm{d}z} e^{-N(z)}$
in Equation~\ref{eq: Peak fluence} is a measure of the probability of having the closest transient neutrino source at redshift $z$.

The relationship between the luminosity distance $d_L$ and the comoving distance $d_c$ is
\begin{linenomath}
\begin{equation} \label{eq: Luminosity distance}
d_L(z) = (1+z)d_c(z),
\end{equation}
\end{linenomath}
where the comoving distance~\citep{Hogg:1999ad} can be expressed as
\begin{linenomath}
\begin{equation} \label{eq: comoving distance}
d_c(z) = \int_0^z \frac{c}{H(z')} \mathrm{d}z',
\end{equation}
\end{linenomath}
with $c$ being the speed of light and $H(z)$ is the Hubble parameter that takes the form
\begin{linenomath}
\begin{equation} \label{eq: Hubble parameter}
H(z)^2 = H_0^2 \left( (1+z)^3 \Omega_m + \Omega_\Lambda \right),
\end{equation}
\end{linenomath}
by assuming a matter and dark energy dominated universe without curvature.
In this work, we adopt the value \mbox{74.03 $\pm$ 1.42\,$\text{km\,s}^{-1}\,\text{Mpc}^{-1}$} for the Hubble expansion rate based on the recent local universe measurements~\citep{Riess:2019cxk}.
The cosmological parameters are taken to be $\Omega_M =$ 0.3 and $\Omega_\Lambda =$ 0.7, respectively.

The probability distribution describing the closest transient neutrino source with respect to redshift is then found by combining Equations~\ref{eq: Number_src} and \ref{eq: comoving distance} and differentiating with respect to the redshift. This yields
\begin{linenomath}
\begin{equation} \label{eq: dNdz}
\frac{\mathrm{d}N}{\mathrm{d}z} = \dot{\rho} T 4\pi \frac{c}{H(z)}  d_c^2.
\end{equation}
\end{linenomath}
By combining Equations~ \ref{eq: Peak fluence}, \ref{eq: Luminosity distance}, and \ref{eq: dNdz} we obtain
\begin{linenomath}
\begin{equation} \label{eq: Peak fluence 2}
\begin{split}
E^2 \Phi_\nu &= \int_0^\infty \mathrm{d}z \frac{\mathcal{E}}{4\pi (1+z)^2 d_c^2(z)} \dot{\rho} T 4\pi \frac{c}{H(z)} d_c^2(z) e^{-N(z)}, \\
&= \dot{\rho} T c \mathcal{E} \int_0^\infty \mathrm{d}z \frac{1}{(1+z)^2 H(z)} e^{-N(z)}.
\end{split}
\end{equation}
\end{linenomath}
Finally, by using Equations~\ref{eq: Number_src}, \ref{eq: comoving distance}, and \ref{eq: Peak fluence 2} we get
\begin{linenomath}
\begin{equation} \label{eq: limit equation}
E^2 \Phi_\nu = \dot{\rho} T c \mathcal{E} \int_0^\infty \mathrm{d}z \frac{1}{(1+z)^2 H(z)} \exp \left( -\dot{\rho} T \frac{4\pi}{3} \int_0^z \frac{c}{H(z')} \mathrm{d}z'^3 \right).
\end{equation}
\end{linenomath}

From Equation~\ref{eq: limit equation} correlated values of the volumetric rate and the bolometric energy are found. The result is shown in Figure~\ref{fig: Limits}, where the 3$\sigma$ discovery potentials for the data with a 15\% increase of rate have been used in order to account for the systematic uncertainties arising when assuming a constant background rate in time and azimuth. The discovery potential in terms of fluence varies with the declination of the source as seen in Figure~\ref{fig: Sen_dis}, and consequently the volumetric rate is declination dependent as well.

\begin{figure*}
	\center
	\includegraphics[width=0.49\linewidth]{Limit_SPS20.pdf}
	\includegraphics[width=0.49\linewidth]{Limit_SPS100.pdf}
	\caption{Upper bounds (3$\sigma$ C.L.) on the volumetric rate of transient neutrino point sources as a function of their bolometric neutrino energy. This is compared to models for high and low luminosity gamma-ray bursts, see~\citep{Murase:2013hh,Liang:2006ci}. The light blue band shows the declination dependence of the discovery potential. Left panel shows result based on sources with a mean energy of 20\,GeV, while right panel is based on a mean energy of 100\,GeV.}
	\label{fig: Limits}
\end{figure*}

\subsection{Systematic uncertainties}
\label{sec:Uncertainties}

We have estimated the systematic uncertainties on the sensitivities by considering simulated data based upon different ice models~\citep{Aartsen:2013rt} and DOM efficiency~\citep{Aartsen:2016nxy}, i.e., the optical efficiency for a DOM to detect a photon. In these calculations, the most conservative ice models, i.e., the models resulting in the largest possible uncertainties, have been used.

The systematic uncertainty on the sensitivity for a source following a subphotospheric energy spectrum with a mean energy of 100\,GeV is just above 25\% for downgoing neutrinos and around 17\% for upgoing events. The largest contribution arises due to the uncertainty on the model describing the bulk ice, e.g., different assumptions for the photon scattering and absorption within the bulk ice.

The systematic uncertainties for the spectrum with mean energy 20\,GeV is slightly smaller, $\sim$20\% for downgoing and $\sim$15\% for upgoing neutrino events. In this case, the uncertainty on the DOM efficiency is the dominant factor.

\section{Conclusion and Outlook}
\label{sec:Outlook}

This analysis has been performed on three years of GRECO data as an offline astrophysical search. This is the first all-flavor low-energy transient source search from IceCube and no astrophysical point sources were found. Upper limits of $\sim 705-2301\, \text{Gpc}^{-3}\, \text{yr}^{-1}$ have been placed on the volumetric rate of the transient neutrino sources, assuming neutrino spectra consistent with that from subphotospheric emission with flare times of up to approximately 600\,s, a mean neutrino energy of 100\,GeV, and a bolometric energy of $10^{52}$\,ergs. Effectively, this applies to gamma-ray bursts following the criteria described. Systematic uncertainties have been included, by using discovery potentials calculated based on varying parameters as described above, ultimately making the limit more conservative.

The analysis presented here is a potentially powerful tool for a real-time alert system. The ability to complement/extend the existing high-energy real-time alerts from IceCube, down to $\sim 10$\,GeV with an all-flavor neutrino selection opens up a new and relatively unexplored energy regime for the contribution of neutrinos to multimessenger astronomy. Benefiting future online, and offline, analyses will be improvements in the angular reconstruction algorithms for low-energy neutrino events as well as improved multi-PMT photon sensors which will be deployed within IceCube as part of the IceCube Upgrade~\citep{Ishihara:2019aao} construction in the early to mid 2020s. A multimessenger campaign centered around an alert system with such improvements may render neutrino emission from high luminosity GRBs identifiable (see Figure \ref{fig: Limits}) in the near future.

\acknowledgments
The IceCube collaboration acknowledges the support from the following: USA {\textendash} U.S. National Science Foundation-Office of Polar Programs,
U.S. National Science Foundation-Physics Division,
Wisconsin Alumni Research Foundation,
Center for High Throughput Computing (CHTC) at the University of Wisconsin{\textendash}Madison,
Open Science Grid (OSG),
Extreme Science and Engineering Discovery Environment (XSEDE),
Frontera computing project at the Texas Advanced Computing Center,
U.S. Department of Energy-National Energy Research Scientific Computing Center,
Particle astrophysics research computing center at the University of Maryland,
Institute for Cyber-Enabled Research at Michigan State University,
and Astroparticle physics computational facility at Marquette University;
Belgium {\textendash} Funds for Scientific Research (FRS-FNRS and FWO),
FWO Odysseus and Big Science programmes,
and Belgian Federal Science Policy Office (Belspo);
Germany {\textendash} Bundesministerium f{\"u}r Bildung und Forschung (BMBF),
Deutsche Forschungsgemeinschaft (DFG),
Helmholtz Alliance for Astroparticle Physics (HAP),
Initiative and Networking Fund of the Helmholtz Association,
Deutsches Elektronen Synchrotron (DESY),
and High Performance Computing cluster of the RWTH Aachen;
Sweden {\textendash} Swedish Research Council,
Swedish Polar Research Secretariat,
Swedish National Infrastructure for Computing (SNIC),
and Knut and Alice Wallenberg Foundation;
Australia {\textendash} Australian Research Council;
Canada {\textendash} Natural Sciences and Engineering Research Council of Canada,
Calcul Qu{\'e}bec, Compute Ontario, Canada Foundation for Innovation, WestGrid, and Compute Canada;
Denmark {\textendash} Villum Fonden, Carlsberg Foundation;
New Zealand {\textendash} Marsden Fund;
Japan {\textendash} Japan Society for Promotion of Science (JSPS)
and Institute for Global Prominent Research (IGPR) of Chiba University;
Korea {\textendash} National Research Foundation of Korea (NRF);
Switzerland {\textendash} Swiss National Science Foundation (SNSF);
United Kingdom {\textendash} Department of Physics, University of Oxford.

\bibliographystyle{jhep}