\section{Introduction}
\label{intro}

The proliferation of the new coronavirus, from now COVID-19, is the current threat to humanity and has spread rapidly around the world starting from January 2020. The 30th of January 2020 is a reference for history because it has been declared by the World Health Organization (WHO) as the official start of international public health emergency, better known as a pandemic. Currently, no countries is immune to the virus and, clearly, the situation appears to be critical. The virus manifests itself after 5-6 days of the onset of the disease with specific and non-specific symptoms. The first are fever, dry cough, sore throat, loss of taste or smell, while the second are fatigue, headache and breathlessness. Unfortunately, COVID-19 has also occurred in animals through transmission between them. Looking back, viruses with similar behavior are Severe Acute Respiratory Syndrome Coronavirus (SARS-CoV) and the Middle East Respiratory Syndrome Coronavirus (MERS Coronavirus) with related major respiratory problems. Currently, the medical protocol takes more than 24 hours to detect the virus in the human body. It is important to detect the disease during the starting phase in order to isolate the infected person because there is no effective cure. Diagnosis can be made through real-time polymerase chain reaction (RT-PCR). RT-PCR is not very reliable due to the high false negative rates and finalization time. Otherwise, COVID-19 can be detected in healthy people due to a false positive. It is clear that the low sensitivity of RT-PCR test is not satisfactory in the current pandemic situation. In some cases, the infected are not recognized in time and do not receive adequate care. As an alternative to RT-PCR, the Thorax Computer Tomography (CT) is probably a more reliable, effective, and faster approach for virus detection and treatment. In almost all hospitals CT image screening is available and can be adopted for a first analysis of the virus. Unfortunately, the Thorax CT images requires a radiologist and a lot precious time is lost. Therefore, the automated analysis of Thorax CT images can speed up the diagnosis in order to help specialist medical staff, above all, and not only, to avoid delays in the start of treatment. In the last few years, deep learning has proven effective for the management, analysis, representation and classification of medical images. In particular, the success of deep neural networks, applied to the image classification task, is connected to different interesting aspects such as the spread of software in terms of open source license, the constant growth of hardware power and the availability of large datasets. Specifically, for the treatment of COVID-19, deep neural networks are adopted both in segmentation and detection phases. However, uncertainty in COVID-19 diagnosis and data imbalance have a decisive impact on performance, hampering model generalization. In order to provide a solution to the above issues, we introduce a framework based on transfer deep learning and ensemble classification for COVID-19 diagnosis. It works based on three integrated stages. A first, which performs image preprocessing operations such as image resize and augmentation. A second, which redesigns and retrains multiple deep neural networks. A third, which combines different predictions provided by deep neural networks with the aim of making the best decision (COVID-19/not-COVID-19). The framework provides the following main contributions:
\begin{itemize}
\item A deep and ensemble learning based framework, to simultaneously address variation between classes and class imbalance for COVID-19 diagnosis task.
\item A framework that provides multiple classification models, based on deep transfer learning.
\item The demonstration that choosing multiple models, suitably combined, is better than a single model and can strengthen the decision during the diagnosis by a specialist doctor.
\item Some experimental greater improvements over existing methods on recent state of art dataset about COVID-19 detection task.
\end{itemize}

The paper is structured as follows. Section \ref{relatedwork} provides an overview of state of art about COVID-19 classification approaches. Section \ref{MM} describes in detail proposed framework. Section \ref{res} provides a wide experimental phase, while section \ref{conc} concludes the paper.

\section{Related work}
\label{relatedwork}
In this section, we briefly analyze the most important approaches working on COVID-19 diagnosis now existing in literature. In this field are included numerous works that address the task according to different aspects. Some offer important contributions about image representation, by implementing segmentation algorithms or new descriptors. Instead, others implement complex mechanisms of learning and classification.

In \cite{loey2020deep} authors propose an architecture in order to improve performance in recognizing COVID-19 from chest radiograph images. It consists of two main components: image augmentation and transfer learning. This combination improves performance measurements such as accuracy, sensitivity, specificity, precision, accuracy, and F1 score.

Authors in \cite{amyar2020multi} present a multitask deep learning model to jointly identify COVID-19 patient and segment COVID-19 lesion from chest computed tomography images. The proposed architecture includes three phases: COVID-19 vs normal vs other infections classification, COVID-19 lesion segmentation, image reconstruction. Furthermore, algorithm flow is based on a common created encoder for the three tasks. It takes a CT scan as input, and its output is then adopted for image reconstruction via a first decoder, to the segmentation via a second decoder, and to the classification of COVID-19 vs normal vs other infections via multilayer perceptron.

In \cite{he2020sample} authors build a public available dataset containing hundreds of CT scans COVID-19 positive and implement a sample-efficient deep learning approach that can obtain high diagnosis accuracy on a limited training set CT images. The approach integrates contrastive self-supervised learning with transfer learning layer, to learn powerful and unbiased features representation. Aim to reducing the risk of overfitting, a large and consistent dictionary on-the-fly based on the contrastive loss to fulfill this auxiliary task is built.

In \cite{el2020novel} authors propose a features selection and voting classifier framework for COVID-19 CT image classification. Firstly, the features are extracted using convolutional neural network (AlexNet). Secondly, a proposed features selection algorithm, Guided Whale Optimization Algorithm (SFS-Guided WOA) based on Stochastic Fractal Search (SFS), is then applied followed by balancing algorithm. Finally, a voting approach, Guided WOA based on Particle Swarm Optimization (PSO), that aggregates different classifier such as Support Vector Machine (SVM), neural networks, K-Nearest Neighbor (KNN) and decision trees predictions, to choose the most voted class in an ensemble learning way, is adopted.

In \cite{shah2020diagnosis} authors design a neural architecture, called CTnet-10, for the COVID-19 diagnosis from CT images. It is formed by a max-pooling layer of dimension $62\times62\times32$ followed by 2 convolutional layers of dimensions $60\times60\times32$, $58\times58\times32$ respectively and a pooling layer of dimension $29\times29\times32$. The last levels are: a flattened layer which is connected out to an fully connected layer of 4096 neurons, in which the dropout layer was used in each of these. Last layer, a single neuron sigmoid and linear, classifies CT scan images as COVID-19 positive or negative. Tests results are compared with known neural networks (DenseNet-169, VGG-16, ResNet-50, InceptionV3, and VGG-19).

In \cite{zhao2020relation} authors build an open-source COVID-19 CT image dataset and a diagnosis method based on multi-task learning and self-supervised learning. To address the overfitting issue, they study two strategies: one is to add additional information including segmentation masks of lung regions and fed them into the feature extraction network.

In \cite{bernheim2020chest} authors does a retrospective study on chest CT scans images with purpose to find the relationship to the time between symptom onset and the initial CT scan. The hallmarks of COVID-19 infection on images are bilateral and peripheral ground-glass and consolidative pulmonary opacities. With a longer time after the onset of symptoms, CT findings are more frequent, including consolidation, bilateral and peripheral disease, greater total lung involvement, linear opacities, crazy-paving pattern, and the reverse halo sign.

In \cite{gozes2020rapid} authors develop an AI-based automated CT image analysis tools for detection, quantification, and tracking of COVID-19. The system utilizes robust 2D and 3D deep learning models, modifying and adapting existing AI models and combining them with clinical understanding. The first step is the lung crop stage, in which the lung region of interest is extracted using a lung segmentation module. The following step detects COVID-19 related abnormalities using deep convolutional neural network architecture. To overcome the limited amount of images dataset, data augmentation techniques (image rotations, horizontal flips and cropping) are applied.

In \cite{zheng2020deep} author propose a 3D deep convolutional neural network, named DeCoVNet, to detect COVID-19 from CT volumes. DeCoVNet is composed of three blocks. First, called network stem, which consists in a vanilla 3D convolution, a batchnorm layer and a pooling layer. Second, composed of two 3D residual blocks (ResBlocks). In each ResBlock, a 3D feature map is passed into both a 3D convolution with a batchnorm layer and a shortcut connection containing a 3D
convolution. The output feature maps are added in an element-wise manner. Third, a progressive classifier (ProClf), which is composed of three 3D convolution layers and a fully-connected (FC) layer. A softmax activation function progressively abstracts the information in the CT volumes by 3D max-pooling and finally directly outputs the probabilities of being COVID-19 or not.

In \cite{ai2020correlation} authors investigate the diagnostic value and consistency of chest CT as comparison with RT-PCR test. For patients with multiple RT-PCR, the dynamic conversion of RT-PCR results (negative to positive, positive to negative) is analyzed as compared with serial chest CT scans for those with a time interval between RT-PCR tests of 4 days or more. Chest CT has a high sensitivity for diagnosis of COVID-19 and may be considered as a primary tool for the current detection in epidemic areas.

In \cite{fang2020sensitivity} authors examine the sensitivity, specificity, and feasibility of chest CT in detecting COVID-19 compared with RT-PCR test. Sensitivity and specificity of chest CT in their various steps are compared using RT-PCR as a gold standard. A reverse calculation approach is applied to chest CT as a hypothetical gold standard and compared to RT-PCR to it point out the flaw of the standard approach. The study want to prove that the sensitivity and specificity of the chest CT in COVID-19 diagnosis and the radiation exposure have to be taken into account together.

\section{Materials and Methods}
\label{MM}
In this section we introduce the proposed framework which is composed of two well known methodologies: deep neural networks and ensemble learning. The main idea is to combine several deep neural networks with purpose to classify images. The result is a set of competitive models providing a range of confidential decisions useful for making choices during classification. The framework is structured into three level. A first, which performs preprocessing in terms of image resize and augmentation. A second, which learns different deep neural networks, previously redesigned for the specific task. A third, in which different models provided by deep neural networks are combined, through ensemble rules, for classification purpose. Finally, the framework iterates through a predetermined number of times in a supervised learning context.

\subsection{Image augmentation}
\label{imgaug}

Many approaches have been developed to address the complications associated with the limited amount of data in machine learning. Image augmentation \cite{shorten2019survey} is a functional technique for increasing and/or changing the size of the training set without acquiring additional images. The concept is basic and consists of duplicating and/or modifying the images with some kind of variation so that more samples can help to train the model. As general idea, the image is augmented in a way that preserves key features for making predictions, but reworked so that the pixels present in some form of noise. The augmentation will be harmful if produces images that are very dissimilar to those used to test the model, so it is clear that this process must be organized in detail. In the proposed framework, we have adopted random reflection, translation, and scaling in order to enhance and augment the image content. As described in the experimental section, this step turns out to be fundamental in improving the performance of the proposed approach.

\subsection{Image resize}
\label{resize}

One of the defects of neural networks concerns the fixed size of the images to be processed. To this end, a resize step is performed based on the input layer dimension claimed by the deep neural networks (details can be found in table \ref{nets} column 5). Most of the networks need this trick but it does not alter the image information content in any way. The size normalization is essential because images of different or large dimensions cannot be processed for the network training and classification stages.

\subsection{Network design and transfer learning}
\label{netmod}

The transfer learning approach has been selected for classification purpose. The basic idea is to transfer the knowledge extracted from a source domain to a destination one, in our case the COVID-19 diagnosis. Generally, a pretrained network is chosen as starting point in order to learn a new task. It is the easiest and fastest solution to adopt the representational power of pretrained deep networks. Clearly, it turns out to be much faster and easier to tune a network with transfer learning than training a new network from scratch with randomly initialized weights. For COVID-19 diagnosis, deep learning architectures are selected based on their structure and performance skills. The goal is to train networks on images by redesign their structures in the final layer according to the needs of the addressed task (two outgoing classes: COVID-19 and not-COVID-19). Table \ref{nets} supports the provided below description of adopted networks.

Resnet18 \cite{Resnet} is inspired by pyramidal cells contained in the cerebral cortex. It uses particular skip connections or shortcuts to jump over some layers. It is composed of 18 layers deep, which with the help of a technique known as skip connection has paved the way for residual networks.

Densenet201 \cite{huang2017densely} is a convolutional neural network with 201 layers deep. Unlike standard convolutional networks composed of $L$ layers with $L$ one-to-one connections between the current layers and the nexts, it contains $\frac{L(L+1)}{2}$ direct connections. Specifically, each layer adopts the feature-maps of all preceding layers and its own feature-maps into all subsequent layers as inputs.

Mobilenetv2 \cite{sandler2018mobilenetv2} is a convolutional neural network with 53 layers deep. It was built based on inverted residual structure with shortcut connections between the thin bottleneck layers. The intermediate expansion layer uses lightweight depthwise convolutions to filter features as a source of non-linearity. Furthermore, non-linearities in the narrow layers are removed with purpose to maintain representational power.

Shufflenet \cite{zhang2018shufflenet} is a convolutional neural network with 173 layers deep designed for mobile devices with very limited computing power. The peculiarity of this architecture concerns the introduction of pointwise group convolution and channel shuffle operations in order reduce computation cost and maintain accuracy.

\begin{table}[!ht]
\centering \caption{Description of adopted pretrained network.}
\footnotesize
\begin{tabular}{l l l l l}
\hline
  Network & Depth & Size (MB) & Parameters (Millions) & Input Size\\ \hline
  Resnet18 & 18 & 44 & 11.7 & 224 $\times$ 224\\ \hline
  Densenet201 & 201 & 77 & 20.0 & 224 $\times$ 224\\ \hline
  Mobilenetv2 & 53 & 13 & 3.5 & 224 $\times$ 224\\ \hline
  Shufflenet & 50 & 6.3 & 1.4 & 224 $\times$ 224\\ \hline
\end{tabular}
\label{nets}
\end{table}

Deep neural networks have been adapted to COVID-19 classification problem. Originally, the main training phase is performed on the Imagenet dataset \cite{deng2009imagenet}, which include a million images divided into 1000 classes. The result consists in a rich features representation for a wide range of images. The network processes an image and provides a prediction about a classes to which it could belong with an attached probability. Commonly, the first layer of the network is the image input layer. The input requires images with 3 color channels. Immediately after is followed by convolutional layers, which work with purpose to extract image features. Particularly, last learnable layer and the final classification layer are adopted to classify the input image. In order to make the pretrained network compliant to classify new images, the two last layers with new layers are replaced. Frequently, the last layer, with related learnable weights, is fully connected. This layer is removed and replaced by a new which is fully connected with outputs related to number of classes of new data (COVID-19 and not-COVID-19). In addition, learning phase of the new layer, compared to the transferred layers, can be speeded up by increasing the rate factors. Optionally, the weights of previous levels can be left unchanged by setting their learning rate to zero. This variation prevents the updating of the weights, during training, and a consequent lowering of the execution time as the gradients of the relative layers do not have to be calculated. This aspect has a strong impact in the case of small datasets in order to avoid overfitting.

\subsection{Ensemble Learning}
\label{ense}

The contribution of different deep neural networks can be mixed in an ensemble context. Considering the set, with cardinality $k$, of images belonging to $x$ classes, to be classified

\begin{equation}
Imgs=\{i_{1},i_{2},\ldots,i_{k}\}
\end{equation}

each element of the set will be treated with the procedure below. Let's consider the set $C$ composed of $n$ deep neural networks

\begin{equation}
C=\{\beta_{1},\beta_{2},\ldots,\beta_{n}\}
\end{equation}

which are combined to classify the images through the set $CN$

\begin{equation}

CN = \begin{bmatrix}
    \beta_{1}i_{1} & \dots & \beta_{1}i_{k} \\
    \vdots & \ddots & \\
    \beta_{n}i_{1} &        & \beta_{n}i_{k}
    \end{bmatrix}
\end{equation}
each $\beta_{n}$ provides a decision $d \in I \{-1,1\}$, where $1$ stands for not-COVID and $-1$ for COVID, with reference to $i_k \in Imgs$. The set of decisions $D$ can be defined as follows

\begin{equation}
D = \begin{bmatrix}
    d_{\beta_{1}i_{1}} & \dots & d_{\beta_{1}i_{k}} \\
    \vdots & \ddots & \\
    d_{\beta_{n}i_{1}} &        & d_{\beta_{n}i_{k}}
    \end{bmatrix}
\end{equation}

it should be noted that each element of the matrix $D$ corresponds to the result of the deep neural network and image combination of the matrix $CN$ in terms of position, such as $\beta_{n}i_{k}\to d_{\beta_{n}i_{k}}$. Furthermore, a score value $s$, $s \in S \{0,\dots,1\}$, is associated with each decision $d$ and represents the posterior probability $P(i|x)$ that an image $i$ could belong to class $x$.  In addition, the set of scores $S$ can be defined as follows

\begin{equation}
S = \begin{bmatrix}
    P(i_{1}|x)_{d_{\beta_{1}i_{1}}} & \dots & P(i_{k}|x)_{d_{\beta_{1}i_{k}}} \\
    \vdots & \ddots & \\
    P(i_{1}|x)_{d_{\beta_{n}i_{1}}} &        & P(i_{k}|x)_{d_{\beta_{n}i_{k}}}
    \end{bmatrix}
\end{equation}

also in this case each element of the matrix $S$ corresponds to the result of the deep neural network and image combination of the matrix $CN$ with related posterior probability in terms of position, such as $\beta_{n}i_{k}\to d_{\beta_{n}i_{k}} \to P(i_{k}|x)_{d_{\beta_{n}i_{k}}}$. At this point, let's introduce the concept of mode, defined as the value which is repeatedly occurred in a given set

\begin{equation}
mode=l+\left(\frac{f_1-f_0}{2f_1-f_0-f_2}\right)
\times h
\end{equation}

where $l$ is the lower limit of the modal class, $h$ is the size of the class interval, $f_1$ is the frequency of the modal class, $f_0$ is the frequency of the class which precedes the modal class and $f_2$ is the frequency of the class which successes the modal class. The columns of matrix $D$ are analyzed by applying the mode, in order to obtain the values of the most frequent decisions. This step is performed in order to verify the best responses of the different deep neural networks, contained in the set $C$. Moreover, the meaning of mode is twofold. First, the most frequent value. Second, its occurrences in terms of indices. For each most frequent occurrence, modal value, the corresponding score from the matrix $S$ is extracted. In this regard, a new vector is generated

\begin{equation}
DS=\{ds_{P(i_{1}|x)_{d_{\beta_{1,\dots,n}i_{1}}}}, \ldots, ds_{P(i_{k}|x)_{d_{\beta_{1,\dots,n}i_{k}}}} \},\end{equation}

where each element $ds$ contains the average of the decision scores with higher frequency, extracted through the mode, in the related column of the matrix $D$. Also, the modal value of each column of the matrix $D$ is stored in the vector $DM$

\begin{equation}
DM=\{dm_{d_{\beta_{1,\dots,n}i_{1}}},\dots, dm_{d_{\beta_{1,\dots,n}i_{k}}}\},\end{equation}

each value $dm$ contains the modal value of the class to which image $i$ could belong with the average probability score $ds$. In essence, this is the class to which an image could belong based on the votes given by different deep neural networks.

\section{Experimental results}
\label{res}

This section describes the experiments performed on public dataset. In order to produce compliant performance, the settings reported in recent COVID-19 classification methods are adopted. The experimental phase is structured in two parts with purpose to address the COVID-19 detection task. The first concerns a comparison with deep neural networks, in order to prove how a multiple model can provide better guidance than a single. While, the second deals with recent methods, which adopt a different logic than the proposed approach.

\subsection{Dataset}
\label{dataset}

The COVID-19 CT adopted dataset is publicly available \footnote{https://github.com/UCSD-AI4H/COVID-CT} and the details are described in \cite{zhao2020covid}. It is composed of 746 Thorax Computer Tomography (CT) images, where 349 contain clinical findings of COVID-19 from 216 patients and 397 are obtained from non-COVID-19 patients. The CT images are a collection selected from COVID-19 related papers published in medRxiv, bioRxiv, NEJM, JAMA, Lancet, and others. The reliability of this dataset has been validated by a senior radiologist, of Tongji Hospital, Wuhan, China, that has worked on diagnosis and treatment of a large number of COVID-19 patients during the period of maximum emergency between January and April 2020.

\subsection{Settings}
\label{settings}

The framework consists in different modules written in Matlab language and extends the code available in \cite{pham2020comprehensive}. The pretrained networks were taken from ImageNet Large-Scale Visual Recognition Challenge (ILSVRC) \cite{russakovsky2015imagenet}. Various combinations of networks were selected, because the test phase was split by applying and not image augmentation. The choice is not random but was made based on single network performance and an in-depth analysis of their architecture (number of layers, applications in literature, etc). A different combination did not provide expected feedback. Densenet201, Mobilenetv2, Resnet18 were adopted with image augmentation, while Shufflenet, Resnet18 and Mobilenetv2 were adopted without image augmentation. Table \ref{nets} shows some important details related to the adopted networks. Among all the computational stages, the training process was certainly the most expensive. As is certainly known, the networks are composed of fully connected layers that make the structure extremely dense and complex. This aspect certainly increases the computational load. In order to compare the results with those obtained in \cite{pham2020comprehensive}, networks were trained by setting the mini batch size to $5$, the maximum epochs to $6$, the initial learning rate to $3 \cdot 10^{-4}$ (constant for the training stage), momentum value to $0.9$, gradient threshold method to L2 norm, factor for L2 regularization (weight decay) to $1 \cdot 10^{-4}$,  minimum batch size to $10$ and the optimizer is stochastic gradient descent with momentum (SGDM) algorithm. $80\%$ and $20\%$ of images are randomly included in train and test set respectively, for a number of iteration equal to $5$ with the aim of calculate the relevance feedback measures reported in table \ref{Metric}. For each training and before each network validation epoch the data were shuffled. The images were converted into RGB space and resized to align with the input format of each pretrained network.
The training was performed with and without image augmentation. Random reflection, translation and scaling were performed for the option of image augmentation. About random reflection, each image was reflected vertically with probability equal to $0.5$ in the top-bottom direction. Again, a horizontal and vertical translation to each image was applied. The translation distance was selected randomly from a continuous uniform distribution within the specified range $[-30,30]$. Similarly, images were scaled vertically and horizontally by selecting in random way the scale factor from a continuous uniform distribution within the specified range $[0.9, 1.1]$.

\subsection{Relevance feedback}
\label{RF}

Table \ref{Metric} summarizes the metrics adopted for the performance evaluation. The goal is to provide an uniform comparison with approaches working on the same task and to understand, from experimental phase, what information can be useful for COVID-19 diagnosis.

\begin{table}[!ht]
\centering
\caption{Evalutation metrics adopted during relevance feedback stage.}
\begin{tabular}{c c}
 &      \\

\hline \textbf{Metric}   & \textbf{Equation}  \\

\hline Sensitivity         & $\frac{TP}{TP + FN}$                     \\

\hline  Specificity      & $\frac{TN}{TN + FP}$                    \\

\hline Accuracy         & $\frac{TP+FN}{TP + FP + TN + FN}$                     \\

\hline $F_{1}$         &$\frac{2 \cdot TP}{2 \cdot TP+FP+FN}$\\

\hline
 &
\end{tabular}
\label{Metric}
\end{table}

The Sensitivity, also known as True Positive rate, concerns the portion of images containing COVID-19 disease elements that are correctly identified. The measure provides important information because highlights the skill to detect images containing the disease and contributes to increase the degree of robustness of result. At the same, it is possible to state of the Specificity, also known as True Negative rate, which instead measures the portion of negatives, images not containing COVID-19 disease elements, that have been correctly identified. Differently, accuracy, a well-known performance measure, is the proportion of true results among the total number of cases examined. In our case provides an overall analysis, certainly a rough measurement compared to the previous ones, about the skill of a classifier to distinguish an image of patient affected to COVID-19 from an image of patient not affected to COVID-19. Furthermore, $F_{1}$ is defined as the combination of precision and recall of the model in term of harmonic mean. In addition, the AUC was calculated using the trapezoidal integration to estimate the area under the ROC curve and represents the measure of performance of a classifier. ROC is a probability curve built by showing the True Positive rate against the False Positive rate with different threshold values. The AUC value is contained in the range between 0.5 and 1, where the value 0.5 represents the performance of a random classifier and the value 1 indicates a perfect one. A high AUC value provides positive classification indications.

\subsection{Discussion}
\label{discussion}

Tables \ref{reswithau} and \ref{reswithoutau} describe the comparison with different deep neural networks respectively applying and not image augmentation. The provided performance can be considered satisfactory compared to different neural architectures. In terms of accuracy, although it provides a rough measurement, we have provided the best result with and without image augmentation. Sensitivity, a measure that provides greater confidence about addressed problem, is very high for both cases. Otherwise, Specificity, which also provides a high degree of information related to the absence of COVID-19 within the image, is the best value for both cases. Regarding the remaining measures, F$_{1}$ score and AUC, considerable values were obtained. Table \ref{comparison} provides comparison results with existing COVID-19 classification methods in term of accuracy. As shown, the proposed approach is only surpassed by Ai et al. \cite{ai2020correlation} and Fang et al. \cite{fang2020sensitivity}. For the remaining methods the performance provided is better. The effectiveness of the results can be attributed to two main aspects: deep neural networks and competitive model for classification. First, the deep neural networks chosen for images learning and classification are the main strong point. Furthermore, the framework provides multiple learning models that certainly constitute a different starting point than a standard approach, in which a single model is provided. This aspect is relevant for improving performance. Second, the classification stage which provides multiple choices in decision making. In fact, at each iteration, the framework selects which networks are suitable for recognizing COVID-19 in the images on test set. Certainly, the computational load is greater but produces better results than a single classification approach. Not negligible issue concerns the image size normalization, with respect to the request of the first layer of the neural networks, before the leaning phase, which not produce a performance degradation. In other cases, degradation of image details, quality and content is due to normalization. Otherwise, the weak point is the computational load even if pretrained networks include layers with already tuned weights. Surely, the time required for training stage is long and computational resources are high but less than a network created from scratch. Moreover, the addressed binary classification has not been greatly disadvantaged by the problem of class imbalance because the relationship between the number of images per class is not very unbalanced. In many cases, data imbalance impacts a low prediction of accuracy for the minority class. It is open problem but the solutions are many such as undersampling of majority class or oversampling of minority classes using image augmentation, weighted loss method by updating the loss function to result in the same way for all classes. This behavior is often seen in medical data due to the limitations of patient samples and cost of acquiring annotated data. Furthermore, in the case of COVID-19 diagnosis it could be relevant as data relating to patients are not completely publics.

\begin{table}[!ht]
\centering
\caption{Classification results with data augmentation.}
\tiny
\begin{tabular}{cccccc}
\hline Network & Accuracy & Sensitivity & Specificity & F$_1$ score & AUC\\
\hline
AlexNet & 74.50 $\pm$ 4.40 & 70.46 $\pm$ 6.37 & 79.05 $\pm$ 8.61 & 0.75 $\pm$ 0.04 & 0.83 $\pm$ 0.04\\\hline
GoogLeNet & 78.97 $\pm$ 3.70 & 75.95 $\pm$ 13.69 & 82.38 $\pm$ 10.53 & 0.79 $\pm$ 0.06 & 0.91 $\pm$ 0.04\\\hline
SqueezeNet & 78.52 $\pm$ 7.56 & 91.56 $\pm$ 7.63 & 63.81 $\pm$ 23.79 & 0.82 $\pm$ 0.04 & 0.90 $\pm$ 0.01\\\hline
ShuffleNet & 86.13 $\pm$ 10.16 & 83.54 $\pm$ 19.89 & 89.05 $\pm$ 5.77 & 0.86 $\pm$ 0.12 & 0.93 $\pm$ 0.06 \\\hline
ResNet-18 & 90.16 $\pm$ 2.36 & 89.45 $\pm$ 7.31 & 90.95 $\pm$ 9.29 & 0.91 $\pm$ 0.02 & 0.96 $\pm$ 0.05 \\\hline
ResNet-50 & 92.62 $\pm$ 4.19 & 91.14 $\pm$ 3.35 & 94.29 $\pm$ 5.15 & 0.93 $\pm$ 0.04 & 0.98 $\pm$ 0.01 \\\hline
ResNet-101 & 89.71 $\pm$ 10.05 & 82.28 $\pm$ 20.09 & 98.10 $\pm$ 2.18 & 0.89 $\pm$ 0.12 & 0.97 $\pm$ 0.03 \\\hline
Xception & 85.68 $\pm$ 6.76 & 90.72 $\pm$ 4.79 & 80.00 $\pm$ 19.64 & 0.87 $\pm$ 0.05 & 0.94 $\pm$ 0.04 \\\hline
Inception-v3 & 91.28 $\pm$ 8.25 & 90.30 $\pm$ 5.12 & 92.38 $\pm$ 11.98 & 0.92 $\pm$ 0.08 & 0.97 $\pm$ 0.02 \\\hline
Inception-ResNet-v2 & 86.35 $\pm$ 5.71 & 88.19 $\pm$ 6.37 & 84.29 $\pm$ 14.50 & 0.87 $\pm$ 0.05 & 0.95 $\pm$ 0.05 \\\hline
VGG-16 & 78.52 $\pm$ 10.02 & 74.68 $\pm$ 30.14 & 82.86 $\pm$ 15.91 & 0.76 $\pm$ 0.17 & 0.91 $\pm$ 0.04 \\\hline
VGG-19 & 83.22 $\pm$ 5.85 & 90.72 $\pm$ 3.19 & 74.76 $\pm$ 12.96 & 0.85 $\pm$ 0.04 & 0.90 $\pm$ 0.05 \\\hline
DenseNet-201 & 91.72 $\pm$ 6.52 & 88.61 $\pm$ 8.86 & 95.24 $\pm$ 4.36 & 0.92 $\pm$ 0.07 & 0.97 $\pm$ 0.03 \\\hline
MobileNet-v2 & 87.25 $\pm$ 10.46 & 95.78 $\pm$ 2.64 & 77.62 $\pm$ 21.63 & 0.89 $\pm$ 0.08 & 0.95 $\pm$ 0.04 \\\hline
NasNet-Mobile & 83.45 $\pm$ 7.36 & 84.81 $\pm$ 2.19 & 81.90 $\pm$ 17.46 & 0.85 $\pm$ 0.05 & 0.94 $\pm$ 0.04 \\\hline
NasNet-Large & 85.23 $\pm$ 8.25 & 79.32 $\pm$ 16.28 & 91.90 $\pm$ 5.77 & 0.84 $\pm$ 0.10 & 0.93 $\pm$ 0.05 \\\hline
\textbf{Ensemble} & 96.38 $\pm$ 4.31 & 95.95 $\pm$ 4.51 & 96.86 $\pm$ 7.03 & 0.97 $\pm$ 0.04 & 0.98 $\pm$ 0.03 \\\hline
\end{tabular}
\label{reswithau}
\end{table}

\begin{table}[!ht]
\centering
\tiny
\caption{Classification results without data augmentation.}
\begin{tabular}{cccccc}
\hline Network & Accuracy & Sensitivity & Specificity & F$_1$ score & AUC\\
\hline
AlexNet  &  86.85 $\pm$ 13.66  &  80.25 $\pm$ 22.49 & 94.29 $\pm$ 4.84 & 0.85 $\pm$ 0.16 & 0.94 $\pm$ 0.04\\\hline
GoogLeNet & 93.83 $\pm$ 6.97 & 96.71 $\pm$ 4.06 & 90.57 $\pm$ 10.53 & 0.94 $\pm$ 0.06 & 0.96 $\pm$ 0.04\\\hline
SqueezeNet & 87.52 $\pm$ 6.45 & 86.84 $\pm$ 10.11 & 88.29 $\pm$ 12.01 & 0.88 $\pm$ 0.06 & 0.94 $\pm$ 0.06\\\hline
ShuffleNet & 95.97 $\pm$ 5.09 & 95.44 $\pm$ 7.47 & 96.57 $\pm$ 2.96 & 0.96 $\pm$ 0.05 & 0.97 $\pm$ 0.03\\\hline
ResNet-18 & 95.44 $\pm$ 8.02 & 98.99 $\pm$ 1.65 & 91.43 $\pm$ 15.25 & 0.96 $\pm$ 0.07 & 0.98 $\pm$ 0.03\\\hline
ResNet-50 & 93.62 $\pm$ 6.17 & 95.57 $\pm$ 6.27 & 91.43 $\pm$ 6.06 & 0.94 $\pm$ 0.06 & 0.98 $\pm$ 0.02\\\hline
ResNet-101 & 93.29 $\pm$ 5.69 & 96.20 $\pm$ 1.79 & 90.00 $\pm$ 10.10 & 0.94 $\pm$ 0.05 & 0.98 $\pm$ 0.02\\\hline
Xception & 91.11 $\pm$ 10.14 & 89.56 $\pm$ 12.55 & 92.86 $\pm$ 7.80 & 0.91 $\pm$ 0.10 & 0.96 $\pm$ 0.03\\\hline
Inception-v3 & 93.62 $\pm$ 5.22 & 96.20 $\pm$ 0.00 & 90.71 $\pm$ 11.11 & 0.94 $\pm$ 0.07 & 0.97 $\pm$ 0.04\\\hline
Inception-ResNet-v2 & 88.59 $\pm$ 7.59 & 89.24 $\pm$ 2.69 & 87.86 $\pm$ 13.13 & 0.89 $\pm$ 0.07 & 0.96 $\pm$ 0.05\\\hline
VGG-16 & 89.26 $\pm$ 8.80 & 92.83 $\pm$ 6.24 & 85.24 $\pm$ 14.45 & 0.90 $\pm$ 0.08 & 0.96 $\pm$ 0.03\\\hline
VGG-19 & 90.16 $\pm$ 7.72 & 87.34 $\pm$ 10.36 & 93.33 $\pm$ 5.77 & 0.90 $\pm$ 0.08 & 0.97 $\pm$ 0.03\\\hline
DenseNet-201 & 96.20 $\pm$ 4.95 & 95.78 $\pm$ 5.27 & 96.67 $\pm$ 4.59 & 0.96 $\pm$ 0.05 & 0.98 $\pm$ 0.03\\\hline
MobileNet-v2 & 95.97 $\pm$ 7.18 & 96.71 $\pm$ 6.04 & 95.14 $\pm$ 8.55 & 0.96 $\pm$ 0.07 & 0.97 $\pm$ 0.05\\\hline
NasNet-Mobile & 89.26 $\pm$ 8.14 & 91.56 $\pm$ 5.12 & 86.67 $\pm$ 13.27 & 0.90 $\pm$ 0.07 & 0.95 $\pm$ 0.06\\\hline
NasNet-Large & 88.59 $\pm$ 7.59 & 90.51 $\pm$ 0.90 & 86.43 $\pm$ 17.17 & 0.90 $\pm$ 0.06 & 0.96 $\pm$ 0.03\\\hline
\textbf{Ensemble} & 96.51 $\pm$ 6.34 & 96.96 $\pm$ 4.79 & 96.00 $\pm$ 8.17 & 0.97 $\pm$ 0.06 & 0.99 $\pm$ 0.03\\\hline
\end{tabular}
\label{reswithoutau}
\end{table}

\begin{table}[!ht]
\centering
\caption{Comparison with other methods.}
\begin{tabular}{cc}
\hline Method & Accuracy\\
\hline
Zhao et al. \cite{zhao2020relation}  &  78.6 \\\hline
Bernheim et al. \cite{bernheim2020chest} & 88.0 \\ \hline
Gozes et al. \cite{gozes2020rapid} &  95.0 \\ \hline
Zheng et al \cite{zheng2020deep} & 90.1 \\ \hline
Ai et al. \cite{ai2020correlation} & 97.0 \\ \hline
Fang et al. \cite{fang2020sensitivity} & 98.0 \\ \hline
Shah et al. (Ctnet-10) \cite{shah2020diagnosis}  & 82.1 \\ \hline
Shah et al. (VGG-19)  \cite{shah2020diagnosis} & 94.5 \\ \hline
\textbf{Ensemble (data augmentation)} & 96.3\\ \hline
\textbf{Ensemble (no data augmentation)} & 96.5 \\ \hline
\end{tabular}
\label{comparison}
\end{table}

\section{Conclusions and Future Works}
\label{conc}

The challenge in COVID-19 detection is especially interesting and, not only, when the data comes from visual information. The complexity of the task is linked to different factors such as the constant increase and variation of data, given that the challenge is in full swing. In support, the convolutional neural networks give a big hand for understand the meaning of information inside the images with the consequent goal of their classification. In this context, we have proposed a framework that combines convolutional neural networks, adapted to COVID-19 detection task, through a transfer learning approach, using an ensemble criteria. The results produced certainly support the theoretical thesis. A multiple model, based on different deep neural networks, compared to a single one is a high discrimination factor. The extensive experimental phase has shown how the proposed approach is competitive, and in some cases surpassing, with respect to state of the art methods. Certainly, the main weak point concerns the computational complexity relating to learning phase, as it is known, takes a long time especially when the data to be processed grows. Future work will certainly concern the study and analysis of convolutional neural networks still unexplored for this type of problem and the application of the proposed framework to additional datasets with purpose to definitively win the challenge to COVID-19.

\section*{Acknowledgements}
Our thinking is for Alfredo Petrosino. He surely would have made a fundamental contribution to winning the battle versus COVID-19 as a man who loved challenges. Finally, we cannot be close to all the families in the world who unfortunately have lost their beloved members.

\bibliographystyle{elsarticle-num}


\section{Appendix}\label{appendix}

\subsection{Network architecture}

In this section we give further details of our neural network architectures. For simplicity we ignore the batch dimension.

Our encoder is an approximation of the AlexNet architecture. Input is a 3-channel colored image with size $128\times128$. At each layer we reduce the image dimensions using strided convolutions, avoiding the need for max pooling operations. We use a Rectified Linear Unit (ReLU) for all non-linear activations. Whilst a deeper network, or one with more modern operations (e.g. Residual connections), would likely have helped overall performance, we aimed to keep are architecture as simple as possible to explore the relative benefits of our proposed learning approach.

\begin{table}[h!]
    \caption{AlexNet encoder architecture details.} \label{tbl:encoder}
    \vspace{.1cm}
    \setlength{\tabcolsep}{0.15cm}
    \centering
    \begin{tabular}{ c c c c c c}
         \hline
         & Layer & Out Size & Batch Norm & Activation \\
         \hline
         C1 & Conv 2D & [64, 30, 30] & True & ReLU \\
         C2 & Conv 2D & [192, 14, 14] & True & ReLU \\
         C3 & Conv 2D & [384, 7, 7] & True & ReLU \\
         C4 & Conv 2D & [256, 4, 4] & True & ReLU \\
         C5 & Conv 2D & [64, 1, 1] & False & ReLU \\
         \hline
         \multicolumn{5}{r}{Parameters: 2,029,056}

    \end{tabular}
    \includegraphics[width=\linewidth]{figures/EncoderArch.png}
\end{table}

Our critic architecture is fully convolutional, and therefore contains no fully connected layers. Unlike the image encoder we find Leaky ReLU activation functions to be more effective. The critic also contains no max pooling layers, using convolutions with either stride=2 or stride=4 to reduce the image size. Input is a 3-channel colored image with size $128\times128$. In particular with the critic, we found normalization to be important. Whilst both Instance Normalization and Batch Normalization proved to be effective, all results in the paper were obtained using Batch Normalization.

\begin{table}[h!]
    \caption{Critic architecture details} \label{tbl:critic}
    \vspace{.1cm}
    \setlength{\tabcolsep}{0.15cm}
    \centering
    \begin{tabular}{ c c c c c c}
     \hline
     & Layer & Out Size & Batch Norm & Activation \\
     \hline
     C1 & Conv 2D & [64, 30, 30] & True & Leaky ReLU \\
     C2 & Conv 2D & [192, 7, 7] & True & Leaky ReLU \\
     C3 & Conv 2D & [384, 4, 4] & True & Leaky ReLU \\
     C4 & Conv 2D & [256, 2, 2] & True & Leaky ReLU \\
     C5 & Conv 2D & [64, 1, 1] & False & Sigmoid \\
     \hline
     \multicolumn{5}{r}{Parameters: 1,883,713}
    \end{tabular}
    \includegraphics[width=\linewidth]{figures/CriticArch.png}
\end{table}

The parameter prediction MLP contains a single shared trunk layer followed by respective parameter branches. No weights are shared on any branch layers. After the shared layers (trunk) we resize the vector to size $n\times d$ where $n$ is the number of object predictions and $d$ number of feature dimensions. We find this effective as it enables each prediction to have its own independent feature transformation. $\mathtt{Center}$ predicts the per-object translation offset in camera space. $\mathtt{RGB}$ predicts an per-object RGB color. $\mathtt{Confidence}$ predicts a per-object confidence which is subsequently used as an opacity value for differentiable rendering with varying objects. $\mathtt{Light}$ predicts a per-scene light direction along a hemisphere above the scene.

For experiments where a canonical mesh orientation, and ground plane are assumed, $\mathtt{Rotation}$ predicts the per-object $i,j$ rotation vector direction. The final rotation $r$ around the $y$-axis is calculated as $r=\mathrm{arctan2}(j, i)$. For experiments where this is not the case (\textsc{Bunny}, \textsc{LineMOD}) we instead using a 6D rotation representation as proposed by \cite{zhou2019continuity}. We choose this representation as it is continuous unlike other representations (Euler angles, rotation matrix, Quaternions etc.), which is more suitable for gradient decent optimizations.

We do not use any normalization in between MLP layers.

\begin{table}[h!]
    \caption{Parameter MLP network architecture.} \label{tbl:mlp}

    \setlength{\tabcolsep}{0.15cm}
    \centering
    \begin{tabular}{ l c c c}
         \hline
          & Layer & Out Size & Activation \\
         \hline
         Trunk & Fully Connected & [64, 1] & ReLU \\
        \hline
         \multirow{2}{5em}{$\mathtt{Center}$} & Fully Connected & [n, 64] & ReLU \\
          & Fully Connected & [n, 3] & Linear \\
        \hline
         \multirow{2}{5em}{$\mathtt{Rotation}$} & Fully Connected & [n, 64] & ReLU \\
         & Fully Connected & [n, 2] & TanH \\
        \hline
         \multirow{2}{5em}{$\mathtt{RGB}$} & Fully Connected & [n, 64] & ReLU \\
         & Fully Connected & [n, 3] & Sigmoid \\
        \hline
         \multirow{2}{5em}{$\mathtt{Confidence}$} & Fully Connected & [n, 64] & ReLU \\
         & Fully Connected & [n, 1] & Sigmoid \\
        \hline
         \multirow{2}{5em}{$\mathtt{Light}$} & Fully Connected & [1, 64] & ReLU \\
         & Fully Connected & [1] & Linear \\
         \hline
         \multicolumn{4}{r}{Parameter range: 4,387 - 48,266}
    \end{tabular}
    \includegraphics[width=\linewidth]{figures/MLPArch.png}
\end{table}

\subsection{Training details}\label{sec:training_details}
All networks were implemented in the PyTorch v1.7 framework. A list of hyperparameters are shown in Tab. \ref{tab:HyperParams}. We trained all models until a convergence threshold based on the image loss is reached. Due to the nature of generator-critic min-max games, we found the critic loss to be an unreliable indicator of convergence. Network training times took between 12-36 hours on a single Nvidia 2080ti depending on task. To allow for a high batch size we utilize virtual batch sizes to avoid running into memory constraints. We achieve this by running multiple forward passes with a smaller batch size and accumulating the gradients before running a backwards pass through our network.

\begin{table}[]
    \centering
    \begin{tabular}{r|c}
        \hline
         Hyperparameter & Value  \\
         \hline
         Batch size & 128 \\
         Optimizer & Adam ($\beta_1$: 0.9, $\beta_2$: 0.999) \\
         Generator Learning rate & 1e-4 \\
         Critic learning rate & 1e-6 \\
         Image Loss $\lambda$ & 0.01 \\
         Critic loss $\lambda$ & 10 \\
         Gradient L2 clipping & 0.5 \\
    \end{tabular}
    \caption{List of hyperparameters used to train networks.}
    \label{tab:HyperParams}
\end{table}

\subsection{Real data experimental setup}\label{sec:exp_setup}
\mycfigure{RealSetup}{Real dataset acquisition pipeline. Images are first captured from either mobile (\textsc{Bunny}) or fixed (\textsc{Sphere Real, Shapes}) cameras (\textbf{a}). Camera alignment is calculated using SfM (\textbf{b}). Using reconstructed geometry and camera exterior orientations we can quickly model the scene (\textbf{c}). Finally, we pass all the parameters into a Differentiable Renderer (DR) for scene parameter learning (\textbf{d}). Note, camera alignment is only used for validation and not during the training or online phase.}

To enable real world experimentation, we recreate physical versions of our synthetic world setups. To enable object learning requires 3 steps, a) data acquisition, b) cameras alignment, c) object and scene modeling (\refFig{RealSetup}).

\textbf{Data acquisition} is performed using two standard Canon DSLR cameras fixed on tripods. We set the validation (novel-view) camera at an $\sim$90\textdegree
angle from the training camera view. To create the validation dataset objects are then placed at random within the scene and the scene is captured from each camera. A training dataset is then created by randomly placing objects into the scene and only capturing data from the training camera. We finally take an image from each camera with no objects in the scene for scene background modeling.

\textbf{Camera alignment}
To get relative exterior orientations between the two cameras we utilize a Structure-from-Motion (SfM) pipeline. Each camera is lifted from the tripod and multiple images take from various viewpoints. We ensure we densely cover the space between the two cameras. The images are input into the SfM software which produces the exterior orientations of the cameras in a common (arbitrary) coordinate system. We further generate a dense point cloud using multi-view stereo to aid object modeling. All SfM processing was performed in MetaShapes v1.6.5 using default settings.

\textbf{Object and scene modeling}
Each shape is modeled inside Blender software, based on the dense point cloud. This produces mesh objects which we pass to the DR. By doing so we ensure our mesh has the same scale as the scene.

We assume a flat ground plane and that objects are in a canonical position, with the only rotation axis being about the objects up-axis. We further transform the camera exterior orientations such that the reconstructed ground plane lies on the $XY$ plane. This is done manually as we have no ground control points to align the bundle adjustment. If a real co-ordinate system is required, ground control points could be placed into the scene avoiding the need for manual transformation (rotation and scaling) of the exterior orientations.

Once the camera alignment is complete, we can then create a ground plane mesh which we texture with the image containing no objects for training and validation cameras respectively. To ensure the same view of the ground plane we model the ground plane to the view frustum of the camera. This can be done by importing the intrinsic matrix estimated from the SfM to a virtual camera. Camera exterior and intrinsic properties are exported along with all mesh objects.

In scenarios where 3D CAD objects are available prior to data capture (e.g. industrial applications where parts are designed in CAD software before being manufactured), the object modeling process is not required.

\textbf{Bunny}
The \textsc{Bunny} dataset demonstrates a more casual data capture method. A single 3D-printed Stanford bunny is filmed from a moving mobile camera. Training data is generated by extracting frames every $n$-frames where $n=10$ in our case. Camera alignment is then computed for a sequence of frames. The first frame in the sequence is used as the input frame and the last frame is used as the novel view. As the motion is continuous, we can get relative camera positions using SfM. Object is then the same as for the \textsc{Real Sphere} and \textsc{Shapes} datasets. In order to model the scene, we require a background image, which contains no object. We found Adobe Photoshop content-aware delete tool to quickly and effectively remove the foreground object, leaving a clean background image which can be used to texture the ground-plane mesh.

\subsection{Validation Experiment}\label{sec:ValidationExperiment}
In this section we propose an experiment to verify the effect the critic has on the network optimization. Our key aim to construct an environment where the input data distribution is from a known distribution $\mathcal{D}$, and therefore the behavior of the critic (to ensure all network outputs belong to this distribution) can be evaluated. We therefore seek to find a differentiable analytical solution which calculates the distance between distribution of the current network output and $\mathcal{D}$. Such a solution would act as an Oracle critic $\mathcal{O}$ (i.e., a perfectly behaved critic).

\textbf{Setup} To aid the experiment, we first reduce our problem space to 2D, where our network must predict a position $t \in [0, 1]$ along the $x$-plane and a luminance $l \in [0, 1]$ (see \refFig{ParamSpace}).

\textbf{Data Distribution} To generate our training dataset, we sample parameters $p\in \mathbb{R}^{t, l}$ from a uniform distribution $\mathcal{D}: \{t, l\} \sim \mathcal{U}(0,1)$. Our optimization then becomes $$\argmin{} D(\hat{p} - \mathcal{D}), $$ where $D$ is a distance metric defined below.

\textbf{Oracle Critic} We define $D$ such that it allows for gradient-based parameter optimization. First, we note that a uniform distribution $\mathcal{U}$ would have a constant density across the parameter space. Therefore, we aim to penalize $\hat{p}$ varying from a constant density. This variation we refer to a discrepancy $D$.

We can measure $D$ by testing for a low discrepancy pattern between all possible subsets of $\hat{p}$ and a random sampling $S \sim \mathcal{D}$. The size $N$ of sampling $S$ here is not strictly important, so long as it represents $\mathcal{D}$ sufficiently. In our experiments $N=300$.

First we calculate the density with a Kernel Density Estimate (KDE) $E$. Whilst it would be tempting to compute $E$ using common methods such as a box or Epechnikov approach, these would not be differentiable. Instead, we use a steep Gaussian kernel $$K = \exp(-(||\hat{p}-S||)^2).$$ Next, we compute KDE estimates $E$ for each point $e \in K$ as $$\frac{1}{2}(e_t^i + e_l^i).$$ Finally, we calculate the variance in the KDE response (discrepancy) $E$ as $$D = \frac{1}{N-1} \sum_{i=0}^{N}{(E_i - \bar{E})^2}.$$ This is the scalar value we minimize in the optimization.

When a simple L2 loss is used and no distribution is enforced (Fig. \ref{fig:ValExpIterations} top) the optimization almost always falls into the local-minima of turning the circle black ($l=0$). This is the same phenomenon we experience in the main paper for $\mathtt{NonCur}$ experiments. We observe this because the majority of gradients for $l$ point towards the black background. As a result, the overall gradient is dominated by this incorrect force. This is visible in the main paper Fig. 6 (left) as parameters falling towards the bottom of the chart. By enforcing $\hat{p} \sim \mathcal{D}$, we see that when the parameters start to fall, they are push back up, where they find more dominant L2 gradients main paper Fig. 6 (right).

A range of optimization visualizations can be found in the supplemental folder under $\mathtt{./val\_experiment}$.

\myfigure{ParamSpace}{Visualization of the parameter space when rendered.}

\subsection{Visualizing training}\label{sec:VisTraining}

In this section we visualize various training scenarios to show the effect of a real critic during training on simple problem. We present two scenarios where the parameter space is 5D ($x,y,r,g,b$) and 6D ($x,y,z,r,g,b$) respectively. For each scenario we demonstrate learning using image loss ($5D_{\mathtt{IMG}}$ / $6D_{\mathtt{IMG}}$), b) curiosity ($5D_{\mathtt{CUR}}$ / $6D_{\mathtt{CUR}}$) and image loss and curiosity ($5D_{\mathtt{OURS}}$ / $6D_{\mathtt{OURS}}$). We generate a dataset by randomly sampling position from a uniform distribution within a fixed bound and each color from a pre-computed color map. This allows us to observe if the critic is able to learn to enforce both the position and color distribution. By visualizing the step-by-step optimization, we note the following observations.

In the 5D scenario, $5D_{\mathtt{IMG}}$ follows dominant and incorrect gradients to match the background color. Furthermore, the gradients push the object towards the edge of the scene, and ultimately out of view of the camera which results in perfectly matching the background. $5D_{\mathtt{CUR}}$ learns to make predictions within the position and color distributions. As no image loss is present, there is no incentive for the network to make a prediction on a given single sample, but instead generate proposals that match the entire distribution. $5D_{\mathtt{OUR}}$ combines both the image loss (i.e., loss on a single given sample) and the critic loss (i.e., loss on the entire distribution). This results in avoiding the local and incorrect minima and results in convergence with respect to the position and color parameters.

Similar observations are made in the 6D scenario. With the additional dimension $6D_{\mathtt{IMG}}$ also finds a local minima by reducing the size of the sphere. This is achieved in practice by pushing the sphere further away from the camera (i.e., $> z$). $6D_{\mathtt{CUR}}$ learns the new position distribution and generates plausible samples that comply with the dataset generation constraints. $6D_{\mathtt{OURS}}$ again finds the correct parameter solution for both position and color. Full training sequences can be found in the supplementary files under $\mathtt{./training\_vis/\{5d, 6d\}}$.

\myfigure{VisTraining}{Example outputs of training visualization experiments for 5D scenario. Columns 1-7 show progress during training. Columns 8 and 9 show the final models network output (col. 9) given an input (col. 8). Full videos can be found at: $\mathtt{./training\_vis/\{5d, 6d\}}$.}

\mycfigure{ValExpIterations}{Renderings of parameters at regular intervals during optimization for \textbf{a)} without curiosity and \textbf{b)} with curiosity. When image loss is used without a curiosity term the parameters degrade into an out-of-distribution solution where luminance $l$ tends towards 0. Optimization visualizations can be found in the supplementary files under $\mathtt{./val\_experiment}$.}

\mycfigure{Results2}{Additional qualitative results.}

\mycfigure{Results3}{Additional qualitative results.}

\mycfigure{Linemod1}{Additional qualitative results on LineMod dataset.}
\mycfigure{Linemod2}{Additional qualitative results on LineMod dataset.}

\section{#1}\label{sec:#2}}
\newcommand{\mysubsection}[2]{\subsection{#1}\label{sec:#2}}
\newcommand{\mysubsubsection}[2]{\subsubsection{#1}\label{sec:#2}}

\newcommand{\scalebox{0.8}{\ding{108}}}{\scalebox{0.8}{\ding{108}}}
\newcommand{\scalebox{0.8}{\ding{117}}}{\scalebox{0.8}{\ding{117}}}
\newcommand{\scalebox{0.8}{\ding{116}}}{\scalebox{0.8}{\ding{116}}}
\newcommand{$\blacktriangleleft$}{$\blacktriangleleft$}
\newcommand{\scalebox{0.8}{\ding{108}}}{\scalebox{0.8}{\ding{108}}}
\newcommand{\ding{72}}{\ding{72}}
\newcommand{\ding{89}}{\ding{89}}
\newcommand{\ding{88}}{\ding{88}}
\newcommand{\ding{70}}{\ding{70}}
\definecolor{colorA}{HTML}{4285f4}
\definecolor{colorB}{HTML}{ea4335}
\definecolor{colorC}{HTML}{fbbc04}
\definecolor{colorD}{HTML}{34a853}
\definecolor{colorE}{HTML}{ff6d01}
\definecolor{colorF}{HTML}{46bdc6}
\definecolor{colorG}{HTML}{000000}
\definecolor{colorH}{HTML}{777777}

\newcommand{\colorIcon}[2]{\textcolor{color#1}{\csname icon#2\endcsname}}
\newcommand{\colorIconAndName}[2]{\colorIcon{#1}{#2}~\text{#2}}
\newcommand{\nameInColor}[1]{\textcolor{color#1}{\texttt{#1}}}
\newcommand{\shortNameInColor}[2]{\textcolor{color#1}{\texttt{#2}}}
\newcommand{\nameAndIcon}[1]{``#1'' (denoted~{\csname icon#1\endcsname})}


\section{INTRODUCTION}
It has building evidence that the Galaxy experienced several major
merger events and left imprints on the stellar halo, in the forms
of streams and overdensities in spatial position
and/or clumps and groups in integrals of motion (e.g., energy, angular
momenta, actions). Two famous mergers are the recent
Sagittarius (Sgr) and the ancient Gaia-Sausage-Enceladus (GSE) dwarf
galaxies. The former is witnessed by the existence of leading and
trailing arms in a large sky coverage \citep{Ibata94}, while
the latter was discovered as the so-called ``Gaia-Sausage" structure
in the radial versus azimuthal velocity space \citep{Belokurov18,Myeong18}
without any dense stream detected so far.
It is proposed that the GSE galaxy collided nearly head-on with
the Galaxy and was disrupted about $6-11$ Gyr ago, leaving a giant cloud of
intermediate-metallicity ($\feh \sim -1.3$) stars on highly radial orbits
\citep{Helmi18,Haywood18,Fattahi20}.
Tracing these events and searching for their imprints are
our long-term goals in order to understand the building blocks of
the Galactic halo, and perhaps also the Galactic disk,
within the context of $\Lambda$CDM cosmology.

As a recent merger with the Galaxy at a relatively low inclination, the debris of Sgr
are found all over the sky
with Galactocentric distances ranging from  15 kpc to
130 kpc \citep{Majewski03,Belokurov14,Sesar17}.
Being on a high eccentricity orbit, the ancient GSE merger left its debris
in a limited region on the sky within 20-30 kpc \citep{Iorio19},
and delivered large amounts of stars
into the solar neighborhood \citep{Belokurov18,Myeong19}.
Interestingly, \cite{Gallart19} suggested that low-$\alpha$ stars found
by \cite{Nissen10} within 0.3 kpc of the Sun
are also debris of the GSE merger. The
Sgr merger event provides a dominant component
in the outer halo ($R>30$ kpc), while the accreted component in the inner halo (and the
Solar neighborhood) mainly comes from the GSE merger
\citep{Naidu20}.

Meanwhile, \cite{Donlon19} suggested that the GSE merger may be
the same as the ``Virgo Radial Merger (VRM)," which includes the Virgo Overdensity,
the Virgo Stellar Stream, and many halo streams. They claimed that the
Virgo Overdensity, the Hercules-Aquila Cloud and
the Eridanus-Phoenix overdensity are probably
unmixed portions of the GSE tidal debris. A systematic search for
this connection has not been carried out and will become possible
with the launch of the future 2-m Chinese Space Survey Telescope (CSST) project \citep{Zhan11} taking into account
of its advantage
of both deep field and high spatial resolution observations.
Since stars belonging to the same structure share similar
abundance patterns \citep{Freeman02},
the chemical mapping of halo stars clustering in both spatial position and
integrals of motion by the CSST project will open the door to reconstruct
the Galactic past merging history.

Last but not least, it is proposed that the GSE merger event perturbed
the proto-disk and produced
a special metal-rich halo-like population in the biggest splash
\citep{Belokurov19}, which heats disk stars up to the Galactic halo.
In support for this scenario, a low angular momentum ($\Vphi<100 \kmprs$)
stellar population with $\feh>-1$ was
found by \cite{Bonaca17}, and was classified as
a ``heated thick disk", rather than an ``in-situ halo."
Later, more works \citep{Haywood18,DiMatteo19,Amarantea20} reported
the existence of Sausage-kinematic stars with thick-disk chemistry.
They were named ``splash" stars and become important
imprints of the GSE merger event on the Galactic disk.

In this paper, we aim to search for additional imprints of the GSE merger event
using the LAMOST low-resolution survey, in combination with Gaia DR2,
based on the chemical and kinematical properties of disk-metallicity stars.
Since the LAMOST survey provides the largest dataset of stars in
the Galactic disk,
it is favorable to investigate the influence of the GSE accretion
event on the evolution of the Galactic disk. It is expected that more
``splash" stars with disk chemistry can be found so that new
merging imprints may be detected from a good statistic analysis on big data.
Meanwhile, it is interesting to probe if there is
a metal-rich tail of the GSE galaxy extending from the main component
at $\feh \sim -1.3$, which will provide an important implication for its
connection with the VRM because the Virgo Overdensity
is metal rich and extends right into the Galactic disk.

\section{Selection of star samples}
The first stage of the LAMOST spectroscopic survey \citep{Zhao06,Zhao12,Cui12,Deng12}
was carried out in a low-resolution (R$\sim$2000) mode with a wide
spectral wavelength coverage of $3800-9000\,\AA$.
Based on abundances from APOGEE DR14 by \cite{Ting19},
\cite{Xiang19} developed a data-driven neural network and
derived stellar parameters and $\afe$ abundances for 8,162,566 stars
from LAMOST DR5. This is an internally consistent
dataset (hereafter DDPayne-LAMOST DR5) and thus becomes a good
sample for a statistic study.

From this catalog, we select stars with $\feh >-1.5$ and require
the errors in both $\feh$ and $\afe$ to be less than 0.1 dex.
We divide the dataset into three samples according to the gravity range,
$2.6<logg<3.3$ (LGB, low red-giant branch), $2.2<logg<2.6$ (RC, red clump)
and $0.0<logg<2.2$ (UGB, upper red-giant branch). With this division,
the analysis is based on more internally-consistent abundances within each
sample than in the whole dataset by reducing the
dependence of abundance on stellar parameters. Meanwhile, results from different samples can be compared and checked independently to draw a reliable conclusion. We do not include dwarf stars
from DDPayne-LAMOST DR5
because their abundances are less reliable
due to the uncertainty in the trained dataset of APOGEE DR14
\citep{Zasowski19}.

APOGEE (Apache Point Observatory Galactic Evolution Experiment)
\citep{Majewski17} is a medium-high resolution (R$\sim$22,500)
spectroscopic survey in the near-infrared spectral range
(H band, 15700-17500\AA). The most recent version is
APOGEE DR16, which provides more reliable parameters and
elemental abundances than APOGEE DR14 thanks to the updated pipeline {\rm ASPCAP} \citep{Garca16} using a grid of only MARCS stellar
atmospheres \citep{Gustafsson08} and a new H-band line list \citep{Smith20}.

We cross-match DDPayne-LAMOST DR5 with APOGEE DR16
and carry out calibrations of $\feh$ and $\afe$ for the selected three
samples. We limit stars with signal-to-noise larger than 80
and a good fitting quality of $\chi_{ASPCAP} < 25$ from APOGEE DR16 to get reliable abundances
for calibration. Fig.~1 shows the comparison of $\feh$, $\afe$
and calibrations by linear fitting between the two datasets for the three samples.
Corrections for $\feh$ are applied to
the three samples of DDPayne-LAMOST DR5 based on the individual
calibration of each sample so that the results from
DDPayne-LAMOST DR5 and APOGEE DR16 can be compared when necessary.
We do not correct $\afe$ ratios because there are
significant scatters in the calibrations, and we prefer to have internally
consistent $\afe$ within each sample from DDPayne-LAMOST DR5 for separating the thin and thick
disks in the following analysis.

We make use of distances from the StarHorse catalog for LAMOST
DR5 by \cite{Queiroz20} with the relative error in distance
less than 10\%. Radial velocities with errors less than $10 \kmprs$
are taken from LAMOST DR5 and proper motions with errors less than $0.25\,mas/yr$ are from {\it Gaia} DR2 \citep{Babusiaux18}.
We check the quality of the {\it Gaia} astrometry
with the renormalized unit weight error of $RUWE<1.44$ \citep{Lindegren18}.
Spatial velocities and orbital parameters are calculated based
on the publicly available
code {\it Galpot}. We employ the default potential of {\it MilkyWayPotential}
provided by \cite{McMillan17}, the solar position of
$R=8.2$ kpc, and the circular
speed $V_{c} = 233.1 \kmprs$ \citep{McMillan11}.
The peculiar velocity of the Sun is ($U_{\odot}, V_{\odot}, W_{\odot}) = $ (11.1, 12.24, 7.25) km/s \citep{Schonrich10}.

With abundances and velocities available, we select stars
with $-100<\Vphi<50 \kmprs$ and $\feh>-0.8$ as our targets
based on the characteristics of the Gaia-Sausage structure,
and we name them ``metal-rich Sausage-kinematic (MRSK)" stars in this work.
The choice of the upper limit of $50 \kmprs$ aims to
avoid stars at $\Vphi \sim 100 \kmprs$, which is
the transition between the
lower boundary of the thick disc and the upper
boundary of the splash box \citep{Belokurov19}.
We adopt a lower limit of $-100 \kmprs$ in order to include
more halo-kinematic disk-metallicity stars. This selection will not
include stars from the Sequoia dwarf galaxy, which has
retrograde orbits with $\Vphi<-100 \kmprs$ and is relatively metal poor
with $\feh < -1.5$ as shown in \cite{Belokurov19}.

Fig.~2 shows the $\VR$ versus $\Vphi$ diagrams and the contour maps
of $\feh$ versus $\afe$ for all stars
with $\feh>-0.8$ and $R>5$ kpc in the three samples. The latter limit is adopted
to avoid stars from the bar/bulge region ($R<5$ kpc). Within the two red
lines, the existence of MRSK stars is obvious.
For normal disk stars with $\feh > -0.8$ and $\Vphi>150 \kmprs$,
the separation between the
high-$\alpha$ thick disk and the low-$\alpha$ thin disks
is shown by the red line of $\afe=-0.08\feh+0.12$ in the $\feh$
versus $\afe$ diagram.
This line will be adopted to pick out
low-$\alpha$ MRSK stars from the high-$\alpha$ group in
the following analysis.
Here we try to adopt a strict criterion to select
the low-$\alpha$ group by lowering the red line close to the boundary
of the low-$\alpha$ contours for normal thin-disk stars.

\begin{figure*}
\centering
\includegraphics[scale=1.0]{f1.eps}
\caption{Comparison of $\feh$ and $\afe$ between DDPayne-LAMOST DR5 and APOGEE DR16 for the three samples (LGB, RC and UGB). Red dash lines are the one-to-one relations and solid lines are calibrations derived from linear fits to the data.}
\label{f1}
\end{figure*}

\begin{figure*}
\centering
\includegraphics[scale=1.0]{f2.eps}
\caption{Left: $\VR$ versus $\Vphi$ diagram for stars with $\feh>-0.8$ and $R>5$ kpc in the three samples. The two red lines show the definition of metal-rich Sausage-kinematic (MRSK) stars within $-100 <\Vphi< 50 \kmprs$. Right: Contour maps of $\feh$ versus $\afe$ for the three samples. Red lines are defined to separate between the thin and thick disks, which will be used to pick out the low-$\alpha$ group from the high-$\alpha$ group among MRSK stars.}
\label{f2}
\end{figure*}

\section{Two groups of MRSK stars with high-$\alpha$ and low-$\alpha$ ratios}
Fig.~3 shows the $\feh$ versus $\afe$ diagrams for MRSK
stars ($-100 <\Vphi< 50 \kmprs$ and $\feh > -0.8$)
in the three samples.
The high-$\alpha$ group above the red line is expected to be
the heated thick disk population, consisting of
so-called ``splash" stars, due to the massive merger event of the GSE galaxy
as suggested by \cite{Belokurov19}.
To our interest, a significant fraction of
MRSK stars have low-$\alpha$ ratios, and they exist in
three different samples. The low-$\alpha$ group extends to a
higher metallicity of $\feh\sim0.0$ (or above) than the
high-$\alpha$ group in each sample.
The existence of low-$\alpha$ MRSK stars is not previously
reported in the literature, and there are totally 1534 such stars
in the three samples. We check that their Sausage-kinematic
characteristics persists when we adopt the distance
by inverting the parallax of {\it Gaia} DR2.
It seems unlikely that errors in $\afe$ and/or kinematics
can produce similar distributions in the $\feh$ versus $\afe$ diagram
for three different samples. Therefore, we suggest that this population
of low-$\alpha$ MRSK stars is real.

As an independent check, we investigate if the low-$\alpha$ MRSK group
exists in the APOGEE survey by applying the same procedure
to APOGEE DR16 for stars in the gravity range of $0<logg<3.8$.
We adopt an extra criterion of $DEC > -10 \deg$, as similar as
the sky coverage of the LAMOST survey, and with the purpose to avoid stars
from Magellanic Clouds and Sagittarius Streams in the southern sky.
Due to the small number of this population, we do not divide
this dataset into three samples (LGB, RC and UGB) accordingly.
As seen in Fig.~4, MRSK stars in APOGEE DR16 show both
high-$\alpha$ and low-$\alpha$ ratios with an obvious gap,
according to the same separation line as that of DDPayne-LAMOST DR5.
The majority of MRSK stars belong to the high-$\alpha$ group, but
the low-$\alpha$ group (consisting of 410 stars) is significant.
Also, the low-$\alpha$ group in this sample
has an similar extension toward super solar metallicity as in DDPayne-LAMOST DR5.
If the systematic shift of 0.03 dex in $\afe$ to a lower value
in APOGEE DR16 (see Fig.~1) is taken into account, the star
number of the low-$\alpha$ group is slightly reduced (398 stars).
In addition, the separation between the high-$\alpha$ and
the low-$\alpha$ groups is more clear in APOGEE DR16
than in DDPayne-LAMOST DR5.
Thus, the existence of a low-$\alpha$ MRSK group is
significant and real in this independent dataset.

Finally, with close inspection on the data in \cite{Belokurov19}, we find that
the $\feh$ versus $\Vphi$ trend in their Figure 2 also shows a hint on
the existence of a metal-rich tail of the GSE galaxy. The distribution of stars
with $\Vphi \sim 0 \kmprs$ seems to cross the splash box and
extends to the metal-rich end until $\feh \sim 0$. Moreover, in
the $\VR$ versus $\Vphi$ diagram as coded with $\afe$,
Sausage-kinematic stars have both high-$\alpha$ and low-$\alpha$ ratios.
This supports the reality of the low-$\alpha$ MRSK group.

\begin{figure*}
\centering
\includegraphics[scale=1.0]{f3.eps}
\caption{$\feh$ versus $\afe$ diagrams for the three samples in DDPayne-LAMOST DR5 data. The red lines are the same as in Fig. 2.}
\label{f3}
\end{figure*}

\begin{figure*}
\centering
\includegraphics[scale=1.0]{f4.eps}
\caption{Upper: $\VR$ versus $\Vphi$ (upper) diagram for stars with $0.0 < \logg < 3.8$, $\feh>-0.8$ and $Dec>-10$ in APOGEE DR16. Lower: Contour map of $\feh$ versus $\afe$ for MRSK stars (with $\feh>-0.8$ and $-100 < \Vphi < 50 \kmprs$) in APOGEE DR16. Lines are the same as in Fig.~2.}
\label{f4}
\end{figure*}

\section{Galactic locations and orbital properties between
high-$\alpha$ and low-$\alpha$ MRSK stars}
\subsection{The $R$, $|Z|$, $\Rapo$ and $\Zmax$ distributions}
The comparison of Galactic locations, in terms
of $R$ and $|Z|$, between the high-$\alpha$ and
low-$\alpha$ groups for the three samples is shown in Fig.~5.
The distributions of the high-$\alpha$
group are reduced by a factor of three in order to have similar peaks
as the low-$\alpha$ group for comparison. Generally,
there is no obvious difference in both $R$ and $|Z|$ distributions
between the high-$\alpha$ and low-$\alpha$ groups for the LGB and
the RC samples. As the luminosity increases (and includes more
distant stars) from the LGB/RC to UGB samples,
we would expect to observe more
old high-$\alpha$ stars than young low-$\alpha$
stars in the UGB sample.
However, this is not the case
and we observe an opposite trend for the UGB sample.
If the low-$\alpha$ group belongs to the locally-born thin-disk population,
it would be located below the splashed high-$\alpha$ thick disk.
Moreover, the GSE merger event happened about 9 Gyr ago,
before the formation of the thin disk. Thus, it is unlikely that these
low-$\alpha$ MRSK stars belong to the young thin-disk population,
and were splashed up to the halo during the GSE merger event. We suggest that
its low-$\alpha$ ratio indicates an origin from the
accreted halo (i.e. the GSE galaxy itself),
rather than from the thin disk.

It is interesting to investigate how far the accreted low-$\alpha$
group would go in the Galaxy and how different in the maximum distance
as compared with the splashed high-$\alpha$ group. In Fig.~6, the histograms of $\Rapo$ and $\Zmax$
are shown for both groups in the three samples. Obviously, the low-$\alpha$
group shows a systematic shift toward a
larger $\Rapo$ by $1-2$ kpc. Specifically, the LGB sample shows
a peak at 8.5 kpc for the high-$\alpha$ group and at 10.5 kpc
for the low-$\alpha$ group despite of their similar $R$ and $|Z|$
distributions. The $\Rapo$ distributions of the low-$\alpha$
group have wide peaks around 10 kpc for RC and UGB samples, extending
up to 15 kpc. The $\Zmax$ distributions generally have peaks around
4 kpc for both groups but the low-$\alpha$ group has extensions toward
higher values in the three samples. These distributions show that
the low-$\alpha$ group is not from the thin disk, which can not
reach a vertical distance as high as the typical thick disk and
would lie far below the splashed thick disk.
Finally, in view of the head-on collision
of the GSE galaxy with the Galaxy, similar peaks at 4 kpc in the
$\Zmax$ distributions
between the low-$\alpha$ and high-$\alpha$ groups may suggest
that they are related, and probably are connected by the same event.
We will discuss this possibility later.

\begin{figure*}
\centering
\includegraphics[scale=1.0]{f5.eps}
\caption{Histograms of $R$ (left) and $|Z|$ (right) for MRSK stars with high-$\alpha$ (black) and low-$\alpha$ (red) ratios in the three samples of DDPayne-LAMOST DR5.}
\end{figure*}

\begin{figure*}
\centering
\includegraphics[scale=1.0]{f6.eps}
\caption{Histograms of $\Rapo$ (left) and $\Zmax$ (right) for MRSK stars with high-$\alpha$ (black) and low-$\alpha$ (red) ratios in the three samples of DDPayne-LAMOST DR5.}
\label{f6}
\end{figure*}

\subsection{The $\Rm-\Zmax$ and $e-\Zmax$ planes}
As compared with its present location $R$, the mean Galactic distance
$\Rm$, as defined to be $(\Rapo+\Rperi)/2.0$, provides
more information on a star's origin.
Fig.~7 shows the $\Rm$ versus $\Zmax$ and the $e$ versus $\Zmax$ planes
for MRSK stars.  Interestingly, there are three sections
for both the high-$\alpha$ (small dots) and low-$\alpha$
groups. In particular, there is a clear gap along the black solid line
of $\Zmax=0.26\Rm+0.4$.
We note that, normal disk stars with
$\Vphi > 150 \kmprs$ are all located in the lowest panel of $\Zmax<2$ kpc.
That is, stars from the thin-disk population has
$\Zmax<2$ kpc and never show up at $|Z| > 2$ kpc.
This is a further support for the suggestion that the low-$\alpha$ MRSK group
does not belong to the thin disk population. Most MRSK stars from
the low-$\alpha$ group occupy the upper two sections with only a few exceptions.
The upper two sections are further separated by the dash line
of $\Zmax=0.9\Rm+1.4$, with blue open circles for
the middle sections and red circles for the upper sections.
The division lines are arbitrarily drawn by eye to separate the
three sections, and we compare the results within each section.
These features share similar characteristics to those in the $\Rapo-\Zmax$ diagram
of \cite{Amaranteb20}, which are caused by resonant trapping effects of the Galactic bar \citep{Moreno15}.
Interestingly, \cite{Schuster20} had shown that 3D star orbits of the high-$\alpha$ star G18-39 and
the low-$\alpha$ star G21-22 in \cite{Nissen10} are trapped by 2D resonant families V and IX of \cite{Moreno15}. Different orbital families in \cite{Moreno15} would reach different $\Zmax$, which explains the appearance of the three sections in Fig.~7.

The distribution of MRSK stars in three sections was previously reported by
\cite{Haywood18} in the $RD_{max}$ versus $\Zmax$ diagram (their Fig.~8)
for the high-$\alpha$ group. In this work, we found similar three-section
distributions for the low-$\alpha$ group in the $\Rm$ versus $\Zmax$ diagram.
This feature is unique to the GSE merger as a result from its high
eccentricity of $e>0.7$. As we limit stars
with $e>0.9$, five sections appear, while stars with $e<0.2$ show only one
section at the lowest part of $\Zmax<2$ kpc.
In view of this dependence, we show the $e$ versus $\Zmax$ diagrams for
stars in the three sections on the right panels of Fig.~7.
Interestingly, stars in the middle sections show a
clump at $\Zmax=3-5$ kpc along the blue line of $\Zmax=8.5e-3.5$, while the
upper sections show sparse distributions without any clump.
There are 82\%, 64\% and 66\% stars in the middle sections located
within 1 kpc (two blue dash lines) along the solid blue lines
for the three samples.

We investigate the distribution of the vertical velocity $V_z$
for stars within two blue dash lines. There is almost zero $V_z$
with a small dispersion of $45 \kmprs$, significantly lower than that
of $\sim 90 \kmprs$ for stars above the upper blue dashed line and
of $\sim 110 \kmprs$ for the Galactic halo.
In particular, the $V_z$ dispersion of $45 \kmprs$ is typical
for the thick disk. Therefore, these clump stars from
the accreted halo could be regarded as the thick disk in kinematics.
We thus propose a scenario that a pile-up of stars at $\Zmax
=3-5$ kpc by the GSE merger event contributes significantly
to (or may be responsible for) the division
between the thick disk and the halo at $|Z|\sim4$ kpc, widely adopted in
the literature \citep{Kinman11,Fernandez19}.
In this scenario, due to the negligible velocity in both $\Vphi$ and $V_z$,
these clump stars, with either low-$\alpha$ or high-$\alpha$ ratio, will stay
at the maximum vertical distance ($\Zmax$) for a longer time than
in other positions, which leads to a pile-up of stars at this distance.
Specifically, the pile-up at $\Zmax=3-5$ kpc (due to the GSE merger event)
produces an observational feature that there is a high density of
stars at $|Z|\sim4$ kpc, beyond which the star number decreases. Thus
this clump is regarded as the division between the thick disk and the halo.
This suggestion is inspired by \cite{Deason13} who proposed that
a rapid transition in structural properties of the Galactic stellar halo
at the break radius of 20-30 kpc \citep{Watkins09,Deason11}
is due to the pile-up in $\Rapo$ of the tidal debris from a small number
of significant mergers, e.g. the Sgr galaxy. Note that the intermediate-eccentricity orbit
of the Sgr merger event will lead to the pile-up of stars in $\Rapo$.
With high-eccentricity
stars, the GSE head-on merger event tends to cause a pile-up in $\Zmax$, rather than $\Rapo$.
It is important to point out
that both the low-$\alpha$ and
high-$\alpha$ groups contribute to the formation of
the high density of stars at $|Z|\sim4$ kpc, and
thus they could be related to the same merger event.
This suggestion is also consistent with the result by \cite{Helmi18} that the GSE
merger led to the formation of the Galactic inner halo and the thick disk.

\begin{figure*}
\centering
\includegraphics[scale=1.0]{f7.eps}
\caption{$\Rm$ versus $\Zmax$ (left) and $e$ versus $\Zmax$ (right) diagrams for MRSK stars with high-$\alpha$ (black dots) and low-$\alpha$ (others) ratios in the three samples of DDPayne-LAMOST DR5. Low-$\alpha$ MRSK stars in the three sections are shown by black triangles (lowest), blue circles (middle) and red circles (upper).}
\label{f7}
\end{figure*}

\subsection{The $\Lz$ versus $\Etot$ diagram}
Fig.~8 shows the $\Lz$ versus $\Etot$ diagrams for
the low-$\alpha$ group in the three sections, as well as the
the high-$\alpha$ group (small dots). With significant overlapping,
the lowest section lies both sides, and the upper section
distributes more close to the central part around $\Lz\sim 0$ and
has a higher energy than the middle section. Specifically,
the lowest parts of Fig.~8 are occupied by stars in the middle
sections (blue circles) and the highest parts are mainly stars
from the upper sections (red circles). However, the overlapping
of both groups is dominated at $\Etotn \sim -1.7\, km^2s^{-2}$.
Therefore, we favor for the suggestion that no obvious
difference between the low-$\alpha$ and high-$\alpha$ groups is found
in this diagram.

The energy of MRSK stars is above the low boundary of the GSE
galaxy at $\Etotn \sim -1.9 \, km^2s^{-2}$, although its metal poor
component lies at $\Etotn \sim -1.5 \, km^2s^{-2}$ according to
\cite{Horta20}. In addition, the energy of the GSE galaxy is significantly
lower than the Sgr galaxy at $\Etotn \sim -1.0 \, km^2s^{-2}$
according to \cite{Naidu20} (see their Fig.~23).
Due to its higher energy, the Sgr merger event
has significant influence on the Galactic halo, since the majority of
member stars can reach up to the outer Galaxy (beyond $20-30$ kpc).
The lower energy of the GSE merger event makes it contribute
significantly to the formation of the Galactic disk, for example,
producing a clump of stars at $\Zmax=3-5$ kpc.

\begin{figure*}
\centering
\includegraphics[scale=1.0]{f8.eps}
\caption{$\Lz$ versus $\Etot$ diagrams for MRSK stars with high-$\alpha$ (black) and low-$\alpha$ (red) ratios in the three samples of DDPayne-LAMOST DR5. Symbols are the same as in Fig.~7.}
\label{f8}
\end{figure*}

\section{The origin of the low-$\alpha$ MRSK group}
\subsection{Connecting the low-$\alpha$ group to the GSE galaxy}

It is well known that Magellanic Clouds (MCs) and Sgr streams have metal rich
components. In APOGEE DR16, it is possible to compare
the low-$\alpha$ MRSK stars with member stars from MCs and Sgr
in the $\feh$ versus $\afe$ diagram. For this purpose,
stars with ``programnames" of ``magclouds" and ``Sgr" in
APOGEE DR16 are picked out, and we limit stars
within the range of $-100<\Vphi<50 \kmprs$, i.e. Sausage-kinematics.
In Fig.~9, these low-$\alpha$ MRSK stars (red circles)
are compared with Sausage kinematic stars from MCs (green crosses) and Sgr
(blue squares), as well as the metal poor GSE members (black dots).
It shows that the low-$\alpha$ MRSK group is
distinct from MCs at $\feh < -0.8$, but has
some overlapping at $-0.8<\feh < -0.4$. Meanwhile, GSE has higher $\afe$
than Sgr for the whole metallicity range of $-1.4<\feh < -0.2$, and
thus the low-$\alpha$ MRSK group does not belong to the Sgr galaxy.
In short, the low-$\alpha$ MRSK stars in APOGEE DR16 are distinct from the member stars of MCs and Sgr. They belong to the GSE galaxy as addressed in Sect. 4.1.

In Fig.~10, stars with $\feh > -1.5$ and $-100<\Vphi<50 \kmprs$
in DDPayne-LAMOST DR5 data are shown in the $\feh$ versus $\afe$
diagrams with the same symbols and colors for the three sections as in Fig.~7.
Three features in Fig.~10 are interesting and useful for understanding the origin of
the low-$\alpha$ group.
Firstly, there is no difference in $\afe$ among the three sections for
the low-$\alpha$ group. All of them belong to the same
population. Secondly, the contour maps show the main component of
the GSE galaxy at $\feh \sim -1.3$ and the second component
at $\feh \sim -0.5$ with a low-density gap at $\feh \sim -0.8$.
Thirdly, the second component becomes weaker from the LGB, through
the RC, to the UGB samples, as
the weak low-$\alpha$ tail becomes more
significant.
In the UGB sample, there is a clear trend of $\afe=-0.20\feh-0.05$
that connects the low-$\alpha$ MRSK group with the main GSE component
at $\feh \sim -1.3$. This connection provides a strong support for
the suggestion that the low-$\alpha$ MRSK group in DDPayne-LAMOST DR5 is
the metal-rich tail of the GSE galaxy.
By shifting the red line upward by 0.05 dex, the trend of
$\afe=-0.20\feh$ indicates the existence of
the metal-rich tail of the GSE galaxy in the LGB and RC samples.

\begin{figure*}
\centering
\includegraphics[scale=1.0]{f9.eps}
\caption{$\feh$ versus $\afe$ diagrams for Sausage-kinematic stars with $\feh > -1.5$, $R>5$ kpc and $Dec > -10$ in APOGEE DR16. Sausage-kinematic stars from Magellanic Clouds (green crosses) and Sagittarius streams (blue squares) are overplotted.
}
\label{f9}
\end{figure*}

\begin{figure*}
\centering
\includegraphics[scale=1.0]{f10.eps}
\caption{$\feh$ versus $\afe$ diagrams and their contour maps for Sausage-kinematic stars with $\feh > -1.5$ and $R>5$ kpc in the three samples of DDPayne-LAMOST DR5. Symbols are the same as in Fig.~7.}
\label{f10}
\end{figure*}

\subsection{A scenario for the origin of the low-$\alpha$ MRSK group}
As described in Sect. 4, similar properties of the low-$\alpha$ and high-$\alpha$ MRSK groups in the $\Rm$-$\Zmax$ and $e-\Zmax$ diagrams indicate that
they probably result from the same process. However, the proposed
splash process by the GSE merger event for the origin of
the high-$\alpha$ group by \cite{Belokurov19} can not explain
the existence of the low-$\alpha$ group. Since $\alpha$ is an age-sensitive abundance \citep{Delgado19}, the
low-$\alpha$ ratio of this MRSK group indicates a young population, which is
not compatible with the splash process of the GSE merger event happened at an
ancient time.

Without a splash process, \cite{Amarantea20} developed a hydrodynamical
simulation, in which a clumpy Milky-Way-like analogue succeeds to produce
a bimodal disk chemistry with Sausage-kinematics, similar as our data.
In their Fig.~1, stars with $-100 <\Vphi< 50 \kmprs$ have metallicity
range from $\feh=-1.5$ to solar metallicity. Moreover,
both high-$\alpha$ and low-$\alpha$ stars exist in the disk metallicity range
of $\feh > -0.8$ (their Fig.~2). In particular, this scenario
does not need to exclude the GSE merger event. According to \cite{Mandelker14},
clumps in a galaxy can have an ex-situ origin associated to mergers.
It is possible that the GSE merger event had brought
the gas-rich cloud
into the Galaxy at an early time, and produced clumps that developed the
bimodal chemistry such as the high-$\alpha$ and low-$\alpha$ MRSK stars
in later evolution. This is a promising solution to the existence of
low-$\alpha$ MRSK stars, but it requires further study
to fit into the case of the GSE dwarf galaxy in more details.

\section{Summary and Conclusions}
A large sample of stars with disk metallicity in the LAMOST survey
provides us a good opportunity to find additional imprints of
the GSE merger event and to extend the search for its member stars
toward solar metallicity. With this purpose, metal-rich stars with Sausage-kinematics,
i.e. $\feh>-0.8$ and $-100 <\Vphi< 50 \kmprs$, in LAMOST DR5
are analyzed in spatial, chemical and kinematical planes.
We divide the stars into three samples (LGB, RC and UGB) so that our analysis is based on higher internally-consistent abundances within each sample of the low-resolution LAMOST survey.

A group of low-$\alpha$ MRSK stars is found, and it is thought to
be the metal-rich component of the GSE merger. In the UGB sample,
there is a continuous trend of $\afe=-0.20\feh-0.05$ that connects it
to the main body of the GSE galaxy at $\-1.3<\feh<-1.1$.
Moreover, it has larger $\Rapo$ (and $\Zmax$) distributions than
those of high-$\alpha$ MRSK stars, which excludes the possibility
to be the low-$\alpha$ thin disk population.
This is the first report on the existence of a metal-rich tail of
the GSE galaxy taking advantage of the large number of disk-metallicity stars in the LAMOST
survey. It provides important implication on the connection between the GSE galaxy and
the Virgo Radial Merger.

In the $\Rm$ versus $\Zmax$ plane, both the high$-\alpha$ and
low-$\alpha$ groups show three-section distributions
separated by $\Zmax = 0.26\Rm+0.4$ and $\Zmax = 0.90\Rm+1.0$ lines.
This feature is unique to the GSE galaxy due
to its high eccentricity of $e>0.7$. Most MRSK stars are located in the middle
and upper sections, while disk-metallicity stars with $\Vphi > 150 \kmprs$
limit within the lowest section of $\Zmax < 2$ kpc. According to \cite{Amaranteb20}, the three-section distributions in the $\Rm$ versus $\Zmax$ plane are caused by different orbital families presented in \cite{Moreno15}.

We find an interesting clump of MRSK stars at $\Zmax=3-5$ kpc,
along the line of $\Zmax=8.5e-3.5$ in the $e$ versus $\Zmax$ diagrams. Since this clump
has negligible velocity in $\Vphi$ and $V_z$, its stars would spend a longer time
at $\Zmax$ than in other positions, leading to a pile-up of MRSK stars
at $|Z|=3-5$ kpc. This clump corresponds to the widely-adopted
disk-halo transition at $|Z|\sim 4$ kpc, and thus is an interesting imprint left by
the GSE merger event. This is new evidence for the scenario by \cite{Helmi18} that
the GSE merger event led to the formation of the Galactic inner halo and the thick disk.

The splash process of the GSE merger event, as proposed to explain the high-$\alpha$ MRSK group, fails to account for
the low-$\alpha$ group because the ancient GSE merger event happened
before its formation.
Instead, a clumpy Milky-Way-like analogue in the hydrodynamical
simulation by \cite{Amarantea20} can produce a bimodal disk chemistry,
similar as our data, without the need for a splash process.
Based on this scenario, both high-$\alpha$ and low-$\alpha$ MRSK stars
in the LAMOST survey belong to the GSE galaxy.
The ancient GSE merger event had brought the gas-rich cloud
into the Galaxy at an early time, and produced clumps
that developed this bimodal chemistry for MRSK stars in
later evolution. This is a promising solution to the origin
of the low-$\alpha$ MRSK stars, but further study using simulations is
necessary to confirm this scenario.

\section*{Acknowledgements}
This study is supported by the National Natural Science Foundation of China (Grant Nos. 11988101, 11625313, 11890694), National Key R\&D Program of China (Grant No. 2019YFA0405502) and the 2-m Chinese Space Survey Telescope project. We would like to extend our sincere thanks to Dr. Sarah A. Bird for her kind help.

Guoshoujing Telescope (the Large Sky Area Multi-Object Fiber Spectroscopic Telescope, LAMOST) is a National Major Scientific Project has been provided by the National Development and Reform Commission. LAMOST is operated and managed by the National Astronomical Observatories, Chinese Academy of Sciences.\\
Funding for the Sloan Digital Sky Survey IV has been provided by the Alfred P. Sloan Foundation, the U.S. Department of Energy Oce of Science, and the Participating Institutions. SDSS-IV acknowledges support and resources from the Center for High-Performance Computing at the University of Utah. The SDSS web site is www.sdss.org.

SDSS-IV is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS Collaboration including the Brazilian Participation Group, the Carnegie Institution for Science, Carnegie Mellon University, the Chilean Participation Group, the French Participation Group, Harvard-Smithsonian Center for Astrophysics, Instituto de Astrofsica de Canarias, The Johns Hopkins University, Kavli Institute for the Physics and Mathematics of the Universe (IPMU) / University of Tokyo, the Korean Participation Group, Lawrence Berkeley National Laboratory, Leibniz Institut f\"ur Astrophysik Potsdam (AIP), Max-Planck-Institut f\"ur Astronomie (MPIA Heidelberg), Max-Planck-Institut f\"ur Astrophysik (MPA Garching), Max-Planck-Institut f\"ur Extraterrestrische Physik (MPE), National Astronomical Observatories of China, New Mexico State University, New York University, University of Notre Dame, Observatario Nacional/ MCTI, The Ohio State University, Pennsylvania State University, Shanghai Astronomical Observatory, United Kingdom Participation Group,Universidad Nacional Autonoma de Mexico, University of Arizona, University of Colorado Boulder, University of Oxford, University of Portsmouth, University of Utah, University of Virginia, University of Washington, University of Wisconsin, Vanderbilt University, and Yale University.