Deep Learning Approach for Spoken Digit Recognition in Gujarati Language

— Speech Recognition is an emerging field in the area of Natural Language Processing which provides ease for human machine interaction with speech. Speech recognition for digits is useful for numbers oriented communication such as mobile number, scores, account number, registration code, social security code etc. This research paper seeks to achieve recognition of ten Gujarati digits from zero to nine ( ૦ to ૯ ) by using a deep learning approach. Dataset is generated with total 8 native speakers 4 male 4 female with the age group of 20 to 40. The dataset includes 2400 labeled audio clips of both genders. To implement a deep learning approach, Convolutional Neural Network (CNN) with MFCC is used to analyze audio clips to generate spectrograms. For the proposed approach three different experiments were performed with different dataset sizes as 1200, 1800 and 2400. With this approach maximum 98.7% accuracy is achieved for spoken digits in Gujarati language with 98% Precision and 98% Recall. It is analyzed from various experiments that increase in dataset size improves the accuracy rate for spoken digit recognition. No of epochs in CNN also improves accuracy to some extent.

I. INTRODUCTION
Speech is the most common medium used by humans for communication. Humans express their thoughts in the form of speech. Speech recognition is the emerging field in the area of Natural Language Processing [1]. Performance and Accuracy of Speech recognition process depends upon the speaker's mental state, emotions, accent, style and ease of speaking and environment. Human speech includes various parameters such as words, sentences, digits and fillers. Speaker feels more comfortable when he is communicating in his native language. In general, people have good command and control over their native language as they have learnt it from their childhood. Speech recognition with native language provides better performance and accuracy with the use of algorithms [2]. Generally the speech recognition process focuses upon identification of words and sentences. Very little work has been done by researchers regarding digits recognition in Gujarati language. Digits recognition with deep learning includes zero to nine numbers identification using optimized algorithms.
Gujarati is Indo-Aryan language and used as native language for people living in Gujarat state in India. It is spoken by 55.5 million people in India. Gujarati language includes 34 consonants, 12 vowels, 10 digits and 106 special symbols.
Speech recognition related to the Gujarati language is a challenging task due to morphological structure, language barrier, dialects and limited resources [3]. This research paper mainly focuses upon digit recognition in Gujarati language using deep learning approach.
Deep learning is a subfield of machine learning that includes various methods and functions that imitates human behavior and knowledge. Deep learning algorithms apply a hierarchy of nonlinear transformation to the input values and generate statistical models as output that can predict results in future on their own. Deep learning approaches require large amounts of data for training and require more computation resources [4]. There are several forms of neural networks such as Recurrent Neural Network (RNN), Convolutional Neural Network (CNN), Artificial Neural Network (ANN) and Feedforward Neural Network (FNN). Most popular deep neural network is Convolutional Neural Network (CNN) to understand and analyze data in visual form.
In the proposed approach for digit recognition for Gujarati language CNN model was implemented. For the analysis process, dataset was prepared with 8 speakers with 2400 labeled audio clips. Different dataset sizes are used as input to CNN to analyze the impact on accuracy.
II. LITERATURE REVIEW Gunawan, A. (2010) [5] developed a digit recognition system from zero through nine using HMM and MATLAB for GUI. Spoken digits were in an isolated structure with 34 male and female speakers. They have created two different modules for isolated and continuous speech recognition. To check effect of variation of environment on accuracy they have selected noisy and clean environment. In a clean environment, in multispeaker mode for isolated and continuous speech average accuracy found is 96%. For speaker independent mode, average accuracy is 78%. In noisy environments, in multispeaker mode for isolated and continuous speech average accuracy found is 80%. For speaker independent mode average accuracy is 61.6%. They have concluded that recognition accuracy depends upon some factors such as noise rate and microphone quality.
Chapaneri and Jayaswal (2013) [6] have proposed a speaker-independent system for isolated digits using WMFCC. They have used SOLA based techniques and IFDTW to reduce time complexity for recognition. The overall accuracy of the proposed system for recognition 0 to 9 was 99.16%. They have *Corresponding Author. www.ijacsa.thesai.org concluded that the system performs 22 times faster with higher accuracy and lower computational overhead achieved using WMFCC and FIFDTW. Datta et al., (2021) [7] have done comparative analysis for methods to recognize zero to nine digits. For comparison they have analyzed cepstral and MFCC. They have concluded that the feature extraction method MFCC performed with better results compared to cepstral in command based systems. Kurian [10] proposed a framework for digit recognition using neural network. For feature extraction they have used Mel Frequency Cepstral Coefficient (MFCC) and Filter Banks (FB) coefficients. They have applied a proposed model on two different datasets one for spoken digits and another for spoken TV commands for different age groups for a total of 88 Arabic native speakers from which 44 were male and 44 were female. Total first digit database entries were 8800 tokens and the second TV command database includes 10000 instances. It is analyzed that the delta-delta feature used to characterize speech recognition signals gives better accuracy over 96% compared to MFCC and FB.
Dhanashri D. & Dhonde S. (2017) [11] have proposed isolated word speech recognition using HMM and DNN for acoustic modeling. Dataset for digit sequences samples were collected in a quiet environment at 20 kHz with 55 men and 57 women. System accuracy is measured for three different hidden layers such as 300,400 and 500 units. Accuracy achieved for 300 hidden layer units is 83.67%, for 400 its 85.44% and for 500 units its 86.06%. They have analyzed that proposed model accuracy can be improved by increasing the number of hidden layer units. [4] have proposed Arabic digits speech recognition model using Recurrent Neural Network (RNN). Feature extraction is done with MFCC and extracted features are fed into the Long Short Term Memory (LSTM) network. They have collected total 1040 data points from different countries and dialects. For the training phase 840 data points were used from 1040 data points and achieved 94% with minimum loss. In the testing phase total 200 data points were included and 65% accuracy was achieved. They have found that the most recognized digit is '0' with 80% accuracy and the worst recognized digit is '6' with 60% accuracy.
Hamidi M. et al. (2020) [12] have presented analysis of Amazigh digit recognition with IVR system for speakerindependent systems in noisy environments. They have used open source software Asterisk 1.6, Ekiga softphone and CMU Sphinx Tools. They have selected 22 Moroccan native speakers of both genders male and female. They have concluded that minor degradation in accuracy found for SNR as low 3 and 9 lb. Major degradation was found with decoded speech signal and best recognition rate at 3 lb.  [14] have performed comparative analysis for four deep learning architectures such as DBN-DNN, LSTM, TDNN and CNN. For feature extraction MFCC and LPC techniques are used. They have used TIDigits dataset which includes 226 speakers from which 111 men and 114 women each uttered 77 digit sequences. Kaldi speech recognition toolkit is used for training and decoding processes. They have concluded that CNN gives the best recognition performance in a noisy and reverberant environment compared to other three architectures. There are also research instances where the authors worked for printed Gujarati characters, for instance [15]. A comparative study of Literature review is presented in Table I. www.ijacsa.thesai.org

III. METHODOLOGY
The methodology to recognize digits in Gujarati language using CNN includes 6 major steps depicted in Fig. 1.

A. Data Collection
For data collection, 8 native speakers, comprising 4 males and 4 females, of the Gujarati language were selected. All 8 speakers were between the age group of 20 to 40. For recording purposes Audacity software was used and recording is performed in a lab environment. Total 30 occurrences of 0 to 9 digits were spoken by each speaker. Dataset consists of 0 to 9 digits in Gujarati with English language representation listed in Table II. Total 2400 occurrences (8 speakers * 10 digits * 30 occurrences) were used for the proposed speech recognition model. Collected data is split in two parts: Training and Testing. For data division 80:20 ratios are followed from total collected data. From 2400 total occurrences 1920 (80% of 2400) were used for training the model and 480 (20% of 2400) occurrences used for testing purposes. All the audio files are collected and stored in wav format.

B. Signal Processing
Speech signals are non-stationary signal combination of voiced and unvoiced signals. From actual speech signals, it is required to extract voiced segments without silence and noise. Data segments with silence are not useful for feature analysis. To remove noise from signals, audacity has a provision with noise reduction option through which we can control by specifying how much dB value we want to consider as noise from 0 to 24 (maximum) dB values.
For processing, first we have removed noise from recorded audio files. Through Audacity noise reduction parameters can be set with residue data that is removed. It provides ease to select a voice sample window which contains useful information about speech signals. The dataset containing a total of 2400 occurrences of the 10 digits, from 0 to 9, were recorded by the speakers in Gujarati. Noise parameters are set to preferences and selected tracks are stored with reduced noise.

C. Feature Extraction with MFCC
Raw data is not used directly in the model training process because it contains noise and other unnecessary parameters that affect accuracy of the model. Feature extraction is an important step in the speech recognition process. Extracted features are used as input to the model to generate text sequences. In the model, extracted features as input, gives better performance compared to raw audio signals. To extract spectral frames from speech signals, we have used the MFCC technique which is a widely used technique for feature extraction. First audio signal is converted from analog to digital with a sampling frequency 16 kHz. To improve model performance, pre-emphasis increases the energy in higher frequency for better phone detection.
To differentiate phones from audio signals, segments are divided into window width 25ms. From each segment, we have extracted 39 features. Hamming window is used to divide signal into segments without noise in a high frequency region. By applying Discrete Fourier Transform (DFT) to convert the time domain of a signal into frequency domain for easy analysis.
By modeling human hearing property into a model at the feature extraction step can improve the performance of the model. For mapping actual frequency with human hearing frequency, mel scale parameter is used. At lower energy humans can easily identify increased levels compared to higher energy. To adapt the same human behavior in the system, log is applied to Mel-filter.
After DFT transformation, the time domain and frequency domain are inverted. The lowest frequency is the highest frequency in the time domain. So inverse transform IDFT (inverse Discrete Fourier Transform) of log of magnitude of signal is generated which is called a cepstrum. After applying IDFT, MFCC includes 12 coefficients to extract signal sample energy as features.
In the dynamic feature generation step, with 13 features, MFCC technique generates first order derivatives and second order derivatives for each feature so total 26 features generated. Derivatives are calculated with differences between coefficients of audio signals to understand transition. Calculated 39 features generated by MFCC for each audio signal that was included as input to the speech recognition system.

D. Deep Learning Approach
Deep learning includes various sub domains such as computer vision, digital signal processing, music classification and automatic speech recognition. With the popularity of virtual assistants like Alexa, Google Home and Siri, the field machine learning systems are promoted at a higher level for speech to text or action. All these systems work upon audio signal processing using deep learning concepts.

E. CNN
Neural networks include a major part as Convolutional Neural Network. CNN is made up of neurons that have learnable weights and biases. It includes various application implementations such as image processing, object identification and classification. CNN can also be used for speech recognition as it can be applied on acoustic models and define convolution layers to each overlap acoustic frames. It extracts phones, pitch, frequency and gender features from acoustic frames according to weights of the network. Convolution Layer is the first layer of CNN that is the core building block and performs major computational tasks. It is a mathematical operation that merges two sets of information. CNN includes three main aspects: locality, weight sharing and pooling. Before start processing input data it is necessary to reduce non-white noise through locality. Weight sharing is used to improve robustness of model and overfitting and reduce weights. Pooling or subsampling process reduces the size. The input image is resized to optimal level to input into the convolutional layer. CNN Architecture is shown in Fig. 2.
In CNN, a filter or kernel is available that is placed in a sliding window over some of the pixels of the image depending upon kernel size to produce feature maps in different convolutional layers. In this step, through element wise multiplication each original pixel value of the image is multiplied with values in the filter. The sum of multiplication values is generated. If the value found is too low then there is nothing to compare in the activation map by curve detector filter. By adding more filters input data can be verified in depth of the activation map. The second layer called the Activation layer, www.ijacsa.thesai.org ReLu (Rectified Linear Unit) function is applied to increase non-linearity in the CNN as images consist of objects that are not linear. In the third pooling layer, that involves downsampling of features. Common pooling layer uses 2 X 2 max filters with stride of 2 that is non-overlapping filter. Max filter returns max value in the feature from selected region. Last layer is known as a fully connected layer that pooled a full feature map matrix into a single column that was used for processing in a neural network. All the features are combined in a network through fully connected layers.

F. Preparing Model
The dataset consists of a total 2400 occurrences both from male and female speakers. It consists of 10 spoken digits recording from 0 to 9. Recording of the dataset was performed using audacity software with 22050 Hz monophonic 16-bit audio files with .wav format. Google colab is used to implement and execute neural networks to get GPU and TPU as a runtime environment. It is free and easy to use. Librosa python library is used to extract spectrograms from audio files. Total dataset is splitted into training and testing data with 80:20 ratios. From a total 2400 audio files, 1920 files were used in the training session. The algorithm is written in python with TensorFlow library using Keras open source library as python interface. Each spectrogram file is stored with PNG format. Spectrogram images of different digits in Gujarati are depicted in Fig. 3. All spectrogram images belong to classes ranging from 0 to 9. Then data is split into training and testing sets. We have used sequential types to build a layer wise model in keras.
In the training phase, input images of the spectrogram are passed to a 2D convolutional layer (Conv2D) with filter size of 32. There are 2 layers as Conv2D layers. First layer consists of 64 nodes and the second layer consists of 32 nodes. With kernel size of 3, 3 X 3 filter matrix is used. For two layers activation function ReLu (Rectified Linear Activation) is used. A 'Flatten' layer is added between conv2D and dense layer to transform the resultant vector into 1 dimensional vector. Dense is the standard output layer used in neural networks. Mainly 10 different nodes were added in layers ranging from 0 to 9. Resultant vector is passed to a fully connected dense layer with a Softmax function that concludes the output in sum up to 1 to represent it in the form of probabilities. Based on the highest priority, prediction series can be generated. After setting all the parameters for the model, compilation is done with three parameters. Compiling process includes optimizer, loss and metrics. SGD (Stochastic Gradient Descent) optimizer is added from keras as it is an iterative method for optimizing objective function. The 'categorical_crossentropy' is used for loss function which indicates the model performance. Third parameter Accuracy metric is used to show accuracy of model on validation data.

G. Training Model
To train a model, fit () function is used with training data and validation data as test data with number of epochs. Epochs specify the total number of model cycles for data. We can get better model performance by increasing the number of epochs up to some extent. For proposed model, 500 epochs were set to analyze neural network performance.

H. Testing Model
To predict the output based upon test data, predict_generator function is used. It generates an array with 10 digit numbers. Each digit range represents the digit classification according to probability.

IV. RESULT AND DISCUSSION
The experiment was executed with training and testing of data to create a neural network. The proposed approach includes a total 2400 occurrences of spoken digits 0 to 9 in Gujarati language. The model is developed with Python programming language with Keras and Tensorflow libraries. The model is executed three times with varied dataset sizes as 1200, 1800 and 2400 digits. The result analysis with various data sets is depicted in Table III. In the first iteration with 1200 dataset values, received accuracy was 58.6% with 58% precision and 58% recall. In the second iteration, accuracy achieved is 78% with 78% precision and 78% recall. In the third iteration with 2400 occurrences, 98.7% accuracy is achieved with 98% precision and 98% recalls.
It is analyzed that accuracy of recognition improves with dataset size. The result can be viewed in graphical form depicted in Fig. 4.  It is observed that CNN gives better performance with speech recognition compared to other related approaches described in literature review.

V. CONCLUSION
Speech recognition is a very important and useful application to interact with machines using speech. This emerging field includes many challenges like dataset size, noisy environment, age, gender, dialect and mental state of the speaker and many more. Speech wav files for Gujarati digits are stored and features extracted with MFCC technique. In this research paper, digit recognition for Gujarati language is performed using CNN.
Speech is passed and processed as input into multiple convolutional layers. For model training Tensorflow and Keras are used with python. Maximum Accuracy rate of recognition out of three experiments with different datasets size for all digits 0 to 9 is found 98.7%. It is observed that accuracy for speech recognition models can be improved by increasing the size of the dataset. Deep learning approach implemented with CNN trained with a large dataset gives better performance. In the future the proposed model can be enhanced to observe results for major parameters such as age, gender and dialect of the speaker that affect the speech recognition result.


extfieldsofstudy:
Computer Science
Engineering

text:
Multi-task UNet: Jointly Boosting Saliency Prediction and Disease Classification on Chest X-ray Images

Human visual attention has recently shown its distinct capability in boosting machine learning models. However, studies that aim to facilitate medical tasks with human visual attention are still scarce. To support the use of visual attention, this paper describes a novel deep learning model for visual saliency prediction on chest X-ray (CXR) images. To cope with data deficiency, we exploit the multi-task learning method and tackles disease classification on CXR simultaneously. For a more robust training process, we propose a further optimized multi-task learning scheme to better handle model overfitting. Experiments show our proposed deep learning model with our new learning scheme can outperform existing methods dedicated either for saliency prediction or image classification. The code used in this paper is available at https://github.com/hz-zhu/MT-UNet.

Introduction
Recent work in machine learning and computer vision have demonstrated advantages of integrating human attention with artificial neural network models, as studies show that many machine vision tasks, i.e., image segmentation, image captioning, object recognition, etc., can benefit from adding human visual attention (Liu and Milanova, 2018).
Visual attention is the ability inherited in biological visual systems to selectively recognize regions or features on scenes relevant to a specific task (Borji et al., 2012), where "bottom-up" attention (also called exogenous attention) focuses on physical properties in the visual input that are salient and distinguishable, and "top-down" attention (also called endogenous attention) generally refers to mental strategies adopted by the visual systems to accomplish the intended visual tasks (Paneri and Gregoriou, 2017). Early research on saliency prediction aims to understand attentions triggered by visual features and patterns, and thus "bottom-up" attention is the research focus (Borji et al., 2012). More recent attempts, empowered by interdisciplinary efforts, start to study both "bottom-up" and "top-down" attentions, and therefore the terms, saliency prediction and visual attention prediction, are used interchangeably (Sun et al., 2021). In this paper, we use the term saliency prediction as the prediction of human visual attentions allocations when viewing 2D images, containing both "bottom-up" and "top-down" attentions. 2D heatmap is usually used to represent human visual attention distribution. Note that saliency prediction studied in this paper is different from neural network's saliency/attention which can be visualized through class activation mapping (CAM) by Zhou et al. (2016) and other methods (Simonyan et al., 2013;Fu et al., 2019;Selvaraju et al., 2016). With the establishment of several benchmark datasets, data driven approaches demonstrated major advancements in saliency prediction (review in Borji (2019) and ). However, saliency prediction for natural scenes is the primary focus, and more needs to be done in the medical domain. Hence, we intend to study the saliency prediction for examining chest X-ray (CXR) images, one of the most common radiology tasks worldwide.
CXR imaging is commonly used for the diagnosis of cardio and/or respiratory abnormalities; it is capable of identifying multiple conditions through a single shot, i.e., COVID-19, pneumonia, heart enlargement, etc. (Ç allı et al., 2021). There exists multiple public CXR datasets (Irvin et al., 2019;Wang et al., 2017). However, the creation of large comprehensive medical datasets is labour intensive, and requires significant medical resources which are usually scarce (Castro et al., 2020). Consequently, medical datasets are rarely as abundant as that for non-medical fields. Thus, machine learning approaches applied on medical datasets need to address the problem of data scarcity. In this paper, we exploit the multitask learning for solution.
Multi-task learning is known for its inductive transfer characteristics that can drive strong representation learning and generalization of each component task (Caruana, 1997). Therefore, multi-task learning methods partially alleviates some of the major shortcomings in deep learning, i.e., high demands for data sufficiency and heavy computation loads (Crawshaw, 2020). However, to apply multi-task learning methods successfully, challenges still exist, which can be the proper selection of component tasks, the architecture of the network, the optimization of the training schemes and many others (Zhang and Yang, 2021;Crawshaw, 2020). This paper investigates the proper configuration of a multi-task learning model that can tackle visual saliency prediction and image classification simultaneously.
The main contributions of this paper are: 1) development of a new deep convolutional neural network (DCNN) architecture for CXR image saliency prediction and classification based on UNet (Ronneberger et al., 2015), and 2) proposal of an optimized multi-task learning scheme that handles overfitting. Our method aims to outperform the state-of-theart networks dedicated either for saliency prediction or image classification.

Saliency prediction with deep learning
DCNN is the leading machine learning method applied to saliency prediction (Pan et al., 2016;Kümmerer et al., 2016;Jia and Bruce, 2020;Kroner et al., 2020). Besides, transfer learning with pre-trained networks was observed to boost the performance of saliency prediction (Oyama and Yamanaka, 2017;Kümmerer et al., 2016;Oyama and Yamanaka, 2018). A majority of DCNN approaches are for natural scene saliency prediction, and so far, only a few studied the saliency prediction for medical images. By Cai et al. (2018), the generative adversarial network is used to predict expert sonographer's saliency when performing standard fetal head plane detection on ultrasound (US) images. However, the  saliency prediction is used as a secondary task to assist the primary detection task, and thus, the saliency prediction performance failed to outperform benchmark prediction methods in several key metrics. Similarly, by Karargyris et al. (2021), as a proof-of-concept study, the gaze data is used as an auxiliary task for CXR image classification, and the performance of saliency prediction is not reported in the study.

CXR image classification with deep learning
Public datasets for CXR images enabled data driven approaches for automatic image analysis and diagnosis (Serte et al., 2020;. Advancements in standardized image classification networks, i.e., ResNet (He et al., 2016), DenseNet (Huang et al., 2017), and EfficientNet (Tan and Le, 2019), facilitate CXR image classification. Yet, CXR image classification remains challenging, as CXR images are noisy, and may contain subtle features that are difficult to recognize even by experts (Ç allı et al., 2021;Khan et al., 2021).

Multi-task Learning Method
As stated in Section 1, component task selection, network architecture design, and training scheme are key factors for multi-task learning. We select the classification task together with the saliency prediction based on the fact that attention patterns are task specific (Karessli et al., 2017). Radiologists are likely to exhibit distinguishable visual behaviors when different patient conditions are shown on CXR images (McLaughlin et al., 2017). This section introduces our multi-task UNet (MT-UNet) architecture, and derives a better multi-task training scheme for saliency prediction and image classification. Figure 1 shows the architecture of the proposed MT-UNet. The network takes CXR images, x ∈ R 1×H×W , where H and W are image dimensions, as input, and produces two outputs, predicted saliency y s ∈ R 1×H×W , and predicted classification y c ∈ R C , where C is the number of classes. As the ground truth for y s is human visual attention distribution, represented as a 2D matrix whose elements are non-negative and sum to 1, y s is normalized by Softmax before output from MT-UNet. Softmax is also applied to y c before output so that the classification outcome can be interpreted as class probability. For the simplicity of notation, batch dimensions are neglected. The proposed MT-UNet is derived from standard UNet architecture (Ronneberger et al., 2015). As a well-known image-to-image deep learning model, the UNet structure has been adopted for various tasks. For example, the UNet is appended with additional structures for visual scene understanding (Jha et al., 2020), the features from the bottleneck (middle of the UNet) are extracted for image classification tasks (Karargyris et al., 2021), and by combining UNet with Pyramid Net (Lin et al., 2017), features at different depth are aggregated for enhanced segmentation (Moradi et al., 2019). What's more, the encoder-decoder structure of UNet is utilized for multi-task learning, where the encoder structure is used to learn representative features, along with designated decoder structures or classification heads for image reconstruction, segmentation, and/or classification (Zhou et al., 2021;Amyar et al., 2020). In our design, we apply classification heads (shaded in light green in Figure 1), which are added not only to the bottleneck but also the ending part of the UNet architecture. This additional classification specific structures aggregates middle and higher-level features for classification, exploiting features learnt at different depths. The attention heads perform global average pooling operations to the 4D tensors, followed by concatenation, and two linear transforms (dense layers) with dropout (rate=25%) in the middle to produce classification outcomes. The MT-UNet belongs to the hard parameter sharing structure in multi-task learning, where different tasks share the same trainable parameters before branched out to each tasks' specific parameters (Vandenhende et al., 2021). Having more trainable parameters in task specific structures may improve the performance for that task at a cost of introducing additional parameters and increasing computational load (Crawshaw, 2020;Vandenhende et al., 2021). In our design, we wish to avoid heavy structures with lots of task specific parameters, and therefore, task specific structures are minimized. In Figure 1, we use yellow and green shades to denote network structures dedicated for saliency prediction and classification, respectively.

Multi-task Training Scheme
Balancing the losses between tasks in a multi-task training process has a direct impact on the training outcome (Vandenhende et al., 2021). There exist multi-task training schemes (Kendall et al., 2018;Chen et al., 2018;Guo et al., 2018;Sener and Koltun, 2018), and among which, we adopt the uncertainty based balancing scheme (Kendall et al., 2018) with the modification proposed in (Liebel and Körner, 2018). Hence, the loss function is: where L s and L c are loss values for y s and y c , respectively; σ s > 0 and σ c > 0 are trainable scalars estimating the uncertainty of L s and L c , respectively; σ s and σ c are initialized to 1; ln(σ s +1) and ln(σ c +1) are regularizing terms to avoid arbitrary decrease of σ s and σ c . With Equation (1), we know that σ values can dynamically weigh losses of different amplitudes during training, and loss with low uncertainty (small σ value) is prioritized in the training process. L > 0. Given y s and y c with their ground truthȳ s andȳ c , respectively, the loss functions are: where H(Q, R) = −Σ n i Q i ln(R i ) stands for cross entropy of two discrete distributions Q and R, both with n elements; H(Q) = H(Q, Q) stands for the entropy, or self cross entropy, of discrete distribution Q. L s is the Kullback-Leibler divergence (KLD) loss, and L c is the cross-entropy loss. By observing Equation (2) and Equation (3), we know that only the cross entropy terms, H(·, ·), generate gradient when updating network parameters, as the term −H(ȳ s ) in L s is a constant and has zero gradient. Therefore, we extend the method in (Kendall et al., 2018), and use 1 σ 2 to scale a KLD loss (L s ) as that for a cross-entropy loss (L c ).
Although the training scheme in Equation (1) yields many successful applications, overfitting for multi-task networks still can jeopardize the training process, especially for small datasets (Wang et al., 2020). Multiple factors can cause overfitting, among witch, learning rate, r > 0, shows the most significant impact (Li et al., 2019). Also, r generally has significant influences on the training outcome (Smith, 2018), making it one of the most important hyper-parameters for a training process. When training MT-UNet, r is moderated by several factors. The first factor is the use of an optimizer. Many optimizers, i.e., Adam (Kingma and Ba, 2014) and RMSProp (Tieleman et al., 2012), deploy the momentum mechanism or its variants, which can adaptively adjust the effective learning rate, r e , during training. As a learning rate scheduler is often used for more efficient training, it is the second factor to influence r. The influence of r from a learning rate scheduler can be adaptive, i.e., reduce learning rate on plateau (RLRP), or more arbitrary, i.e., cosine annealing with warm restarts (Loshchilov and Hutter, 2016). By observing Equation (1), we know that an uncertainly estimator σ for a loss L also serves as a learning rate adaptor for L, which is the third factor. More specifically, given a loss value L with learning rate r, the effective learning rate for parameters with a scaled loss value L σ 2 is r σ 2 . Decreasing r upon overfitting can alleviate its effects (Smith, 2018;Duffner and Garcia, 2007), but Equation (1) leads to increased learning rate upon overfitting, further worsening the training process. This happens because training loss decreases when overfitting occurs, reducing its variance at the same time. Thus, σ decreases accordingly, which increases the effective learning rate, thus creating a vicious circle of overfitting. More detailed mathematical derivation is presented in Appendix A. This phenomenon can be observed in Figure 2, where changes of losses and σ values during a training process following Equation (1) are presented. We can see from Figure 2(a), at epoch 40, after an initial decrease in both the training and validation losses, the training loss start to acceleratedly decrease while the validation loss start to amplify, which is a vicious circle of overfitting. A RLRP scheduler can halt the vicious circle by resetting the model parameters to a former epoch and reducing r. Yet, even with reduced r, a vicious circle of overfitting can remerge in later epochs.  (1) To alleviate overfitting, we propose the use of the following equations to replace Equation (1): The essence of Equations (4) and (5) is to fix the uncertainty term for one loss in Equation (1) to 1, so that the flexibility in changing effective learning rate is reduced. With the uncertainty term fixed for one component loss, Equations (4) and (5) demonstrate the ability to alleviate overfitting and stabilize the training processing. It is worth noting that Equations (4) and (5) cannot be used interchangeably. We need to test both equations to check which can achieve better performances, as depending on the dataset and training process, overfitting can occur of different severity in all component tasks. In this study, training process with Equation (5) achieves the best performance. Ablation study of this method is presented in Section 5.

Dataset and Evaluation Methods
We use the "chest X-ray dataset with eye-tracking and report dictation" (Karargyris et al., 2021) shared via PhysioNet (Moody et al., 2000) in this study. The dataset was derived from the MIMIC-CXR dataset (Johnson et al., 2019a,b) with additional gaze tracking and dictation from an expert radiologist. 1083 CXR images are included in the dataset, and accompanying each image, there are tracked gaze data; a diagnostic label (either normal, pneumonia, or enlarged heart); segmentation of lungs, mediastinum, and aortic knob; and radiologist's audio with dictation. The CXR images in the dataset are in resolutions of various sizes, i.e., 3056 × 2044, and we down sample and/or pad each image to 640 × 416. A GP3 gaze tracker by Gazepoint (Vancouver, Canada) was used for the collection of gaze data. The tracker has an accuracy of around 1°of visual angle, and has a 60 Hz sampling rate (Zhu et al., 2019). Several metrics have been used for the evaluation of saliency prediction performances, and they can be classified into location-based metrics and distribution-based metrics (Bylinskii et al., 2018). Due to the tracking inaccuracy of the GP3 gaze tracker, location-based metrics is not suited for this study. Therefore, in this paper, we follow suggestions in (Bylinskii et al., 2018) and use KLD for performance evaluation. We also include histogram similarity (HS), and Pearson's correlation coefficient (PCC) for reference purposes. For the evaluation of classification performances, we use the area under curve (AUC) metrics for multi-class classifications (Hand and Till, 2001;Fawcett, 2006), and the classification accuracy (ACC) metrics. We also include the AUC metrics for each class: normal, enlarged heart, and pneumonia, denoted as AUC-Y1, AUC-Y2, and AUC-Y3, respectively. In this paper, all metrics values are presented as median statistics followed by standard deviations behind the ± sign. Metrics with up-pointing arrow ↑ indicates greater values reflect better performances, and vise versa. Best metrics are emboldened.

Benchmark comparison
In this subsection, we compare the performance of MT-UNet, with benchmark networks for CXR image classification and saliency prediction. Detailed training settings are presented in Appendix B.
For CXR image classification, the benchmark networks are chosen from the top performing networks for CXR image classification examined in (El Asnaoui et al., 2021), which are ResNet50 (He et al., 2016) and Inception-ResNet v2 (abbreviated as IRNetV2 in this paper) (Szegedy et al., 2017). Following Karargyris et al. (2021), we also include a state-of-the-art general purpose classification network: EfficientNetV2-S (abbreviated as EffNetV2-S) (Tan and Le, 2021) for comparison. For completeness, classification using standard UNet with additional classification head (denoted as UNetC) is included. Results are presented in Table 1, and We can see that MT-UNet outperforms the other classification networks.
For CXR image saliency prediction, comparison was conducted with 3 state-of-the-art saliency prediction models, which are SimpleNet (Reddy et al., 2020), MSINet (Kroner et al., 2020) and VGGSSM (Cao et al., 2020). Saliency prediction using standard UNet (denoted as UNetS) is also included for reference. Table 2 shows the result, where MT-UNet outperforms the rest. Visual comparisons for saliency prediction results are presented through Table 4 in Appendix C.

Ablation study
To validate the modified multi-task learning scheme, ablation study is performed. The multi-task learning schemes following Equations (1), (4) and (5) are compared, and they are denoted as MTLS1, MTLS2, and MTLS3, respectively. Please note that the bestperforming MTLS3 is used for benchmark comparison in Section 5.1. Figure 3 in Appendix C shows the training process for MTLS2 and MTLS3. With Figures 2 and 3, we can see that overfitting occurs both for MTLS1 and MTLS2, but the overfitting is reduced in MTLS3.
The training processes shown in Figures 2 and 3 are with optimized hyper-parameters. The resulting performances are compared in Table 3 in Appendix C. We can see that MTLS3 outperforms the rest learning schemes both in classification and in saliency prediction.
To validate the effects of using classification head that aggregates features from different depths, we create ablated versions of MT-UNet that use features from either the bottleneck or the top layer of the MT-UNet for classification, denoted as MT-UNetB and MT-UNetT, respectively. Results are presented in Table 3 in Appendix C. We can see that MT-UNet generally performs better than MT-UNetT and MT-UNetB.

Discussion
In this paper, we build the MT-UNet model and propose a further optimized multi-tasking learning scheme for saliency prediction and disease classification with CXR images. While a multi-task learning model has the potential of enhancing the performances for all component tasks, a proper training scheme is one of the key factors to fully unveil its potentiality. As shown in Table 3, MT-UNet with the standard multi-task learning scheme may barely outperform existing models for saliency prediction or image classification.
Several future work could be done to improve this study. The first would be the expansion of the gaze tracking dataset for medical images. So far, only 1083 CXR images are publicly available with radiologist's gaze behavior, limiting extensive studies of gazetracking assisted machine learning methods in the medical field. Also, more dedicated studies on multi-task learning methods, especially for small datasets, can be helpful for medical machine learning tasks. Overfitting and data deficiency are the lingering challenges encountered by many studies. A better multi-task learning method may handle these challenges more easily.
published the multi-modal chest X-ray dataset for this research. This research is supported by Compute Canada and the Natural Sciences and Engineering Research Council of Canada (NSERC).

Appendix A. Mathmatical deriviation of vicious circle for overfitting
Let L ≥ 0 be the loss for a task, T , and σ > 0 be the variance estimator for L used in Equation (1). Therefore, the loss for T following Equation (1) can be expressed as: The partial derivative of L with respect to σ is: During a gradient based optimization process, to minimize L, σ converges to the equilibrium value (σ remains unchanged after gradient descend) which is achieved when ∂L ∂σ = 0. Therefore, the following equation holds when σ is at its equilibrium value, denoted asσ: which is calculated by letting ∂L ∂σ = 0. Let f (σ) = L,σ > 0, we can calculate that: Therefore, we know that f (σ) is strictly monotonically increasing with respect toσ, and hence the inverse function of f (σ), f −1 (·), exists. More specifically, we have: As a pair of inverse functions share the same monotonicity, we know thatσ = f −1 (L) is also strictly monotonically increasing. Thus, when L decreases due to overfitting, we know thatσ will decrease accordingly, forcing σ to decrease. The decreased σ leads to an increase in the effective learning rate for T , forming a vicious circle of overfitting.

Appendix B. Training settings
We use the Adam optimizer with default parameters (Kingma and Ba, 2014) and the RLRP scheduler for all the training processes. The RLRP scheduler reduces 90% of the learning rate when validation loss stops improving for P consecutive epochs, and reset model parameters to an earlier epoch when the network achieves the best validation loss. All training and testing are performed with the PyTorch framework (Paszke et al., 2019). Hyper-parameters for optimizations are learning rate r, and P in RLRP scheduler. The dataset is randomly partitioned into 70%, 10% and 20% subsections for training, validation and testing, respectively. The random data partitioning process preserves the balanced dataset characteristic, and all classes have equal share in all sub-datasets. All the results presented in this paper are based on at least 5 independent trainings with same hyper-parameters. NVIDIA V100 and A100 GPUs (Santa Clara, USA) were used.    Table 4: Visualization of predicted saliency distributions. The ground truth and predicted saliency distributions are overlaid over CXR images. Jet colormap is used for saliency distributions where warmer (red and yellow) colors indicate higher concentration of saliency and colder (green and blue) colors indicate lower concentration of saliency.


A case study of the profit-maximizing multi-vehicle pickup and delivery selection problem for the road networks with the integratable nodes

This paper is a study of an application-based model in profit-maximizing multi-vehicle pickup and delivery selection problem (PPDSP). The graph-theoretic model proposed by existing studies of PPDSP is based on transport requests to define the corresponding nodes (i.e., each request corresponds to a pickup node and a delivery node). In practice, however, there are probably multiple requests coming from or going to an identical location. Considering the road networks with the integratable nodes as above, we define a new model based on the integrated nodes for the corresponding PPDSP and propose a novel mixed-integer formulation. In comparative experiments with the existing formulation, as the number of integratable nodes increases, our method has a clear advantage in terms of the number of variables as well as the number of constraints required in the generated instances, and the accuracy of the optimized solution obtained within a given time.

Introduction
Recently, the profit-maximizing multi-vehicle pickup and delivery selection problem (PPDSP) has gained a lot of attention in the field of practical transportation and logistics Liu et al. (2018); Riedler & Raidl (2018); Asghari & Mirzapour Al-e-hashem (2020); Huang et al. (2020). This problem was first proposed in Qiu et al. (2017), which involves three classical problem models: routing optimization, pickup and delivery, and selective pickup (a.k.a. knapsack ). Solving this problem quickly and optimally can both help improve the operational efficiency of the carriers and contributes to more eco-friendly transportation.
To the best of our knowledge, in the graph-theoretic models constructed in the existing studies Ting et al. (2017); Gansterer et al. (2017); Al-Chami et al. (2018); Ahmadi-Javid et al.   Figure 1c can be regarded as a feasible solution of the request-based model, but cannot be corresponded to in the location-based model.
(2018), the definitions of nodes are based on the pickup and delivery location from the requests. This means that one request needs to correspond to two nodes. However, in application scenarios, there are often plural requests coming from or to arrive at the same location. The number of variables to be required in the existing mixed-integer formulation heavily depends on the number of nodes in the model, and since PPDSP is an N P-hard problem, its computational complexity is exponential with respect to the number of variables.
The motivation of this study is to provide a more reasonable and effective mathematical model for practical logistics and transportation problems. Considering the road networks with the integratable nodes as above, we proposes a new concise model in which the definition of nodes is based on the locations rather than the requests, and give a novel mixed-integer programming formulation that reduces the number of required variables. It is necessary to claim that the solution set of our proposed location-based method is a subset of the solution set of the request-based method. This is because the integrated node can only be passed at most once under the Hamiltonian cycle constraint, where "at most" is due to the objective function of profit-maximizing (i.e., it is possible that any node will not be traversed). In Figure 1, we take the delivery of single truck as an example, and assume that the truck is transported according to the route given in Figure 1a as a feasible solution in this scenario, where the pickup coordinates of request 1 (abbr. Req. 1 pickup coord.) and Req. 2 pickup coord. are the identical locations. We can integrate them as location 1 as shown in Figure 1b, and after integrating all the same coordinates into one location respectively, our model becomes more concise. This integration changes the nodes of the model from repeatable coordinates based on the request to the unique locations, which results in each location can only be visited once. Therefore, the location-based model can also correspond to the feasible solutions in Figure 1a, but not to some of the feasible solutions corresponded to the request-based model, such as the solution given in Figure 1c. Whereas in real logistics application scenarios, making multiple retraces to the same location is rare (e.g. insufficient vehicle capacity, goods cannot be mixed, etc.), and also non-efficient. For such locations with a large number of requests or high loading volume, it is a more common strategy to assign multiple vehicles to them. Therefore, we argue that the location-based model can improve the optimization efficiency although it reduces the range of feasible solutions.

Preliminaries
Given a directed graph G = (V, E), where V = {0, 1, 2, . . . , |V |} is the set of nodes representing the location points (0 is depot) and E is the set of arc, denoted as > od, and given a set of trucks T = {1, 2, . . . , |T |}, we define a type of Boolean variables x t od that is equal to 1 if the truck t passes through the arc > od and 0 otherwise. We denote the load capacity of truck t as c t , and the cost of the truck t traversing the arc > od as l t od , where t ∈ T and o, d ∈ V .
Let R = {1, 2, . . . , |R|} be a set of requests. Each request r (r ∈ R) is considered as a tuple r = w r , q r , f (r), g(r) , where • w r is the the payment that can be received for completing the shipping of request r; • q r is the volume of request r; • f (r) is the loading point of request r; • g(r) is the unloading point of request r, and f : R → V \ {0} (resp. g : R → V \ {0}) is a function for mapping the loading (resp. unloading) point of requests r. We also define another type of Boolean variables y t r that is equal to 1 if request r is allocated to truck t and 0 otherwise, where r ∈ R and t ∈ T . Besides, we denote the number of location points visited by truck t before it reaches Definition 1 (Delivery of Truck t). Delivery of truck t is denoted by D t = r∈R {r | y t r = 1}, where t ∈ T .
Definition 2 (Route of Truck t). Route of truck t is denoted by } as the set of location points visited and departed exactly once by truck t, where the predecessor and successor are the same node is not counted (e.g., even if > vv | x t vv = 1, neither this time can be included in the number of visits or departures of truck t to/from location v.); • Ensure that no subtour exists in S t . 1 Loading Before Unloading ∀i ∈ D t , u t f (i) < u t g(i) .
Capacity Limitation At any time, the total volume of cargo carried by truck t cannot exceed its capacity c t .
Definition 3 (Delivery Routing Solution). Delivery routing solution DS = t∈T {(D t , S t )} that is a partition of R into disjoint and contained D t with the corresponding S t : Denote the set of all possible delivery routing solutions as Π(R, T ).
Definition 4 (Profit-Cost Function). A profit-cost function assigns a real-valued profit to every D t : w : D t → R and a cost to every S t : l : S t → R, where w(D t ) = r∈Dt w r and l(S t ) = > od∈St l t od . For any delivery routing solution DS ∈ Π(R, T ), the value of DS is calculated by In general, a delivery routing solution DS that considers only maximizing profits or minimizing costs is not necessarily an optimal DS. Therefore, we have to find the optimal delivery routing solution that maximizes the sum of the values of profit-cost functions. We define a delivery routing problem in profit-cost function as follows.
Definition 5 (PPDSP). For a set of requests and trucks (R, T ), a profit-maximizing multivehicle pickup and delivery selection problem (PPDSP) is to find the optimal delivery routing solution DS * such that Here we show an example of PPDSP.
Example 1. Assume that two trucks T = {t 1 , t 2 } are responsible for three requests R = {r 1 , r 2 , r 3 } with the following conditions.
• The information about the requests: • The information about the trucks: 4 -The capacities of the trucks c t 1 = 6 and c t 2 = 3.
-The cost matrices of each truck through each arc (δ is depot).
In this example, the optimal delivery routing solution DS * for the given PPDSP is cδ})}, and its value is 11.

Problem Formulation
In this section, we present the following mixed-integer programming (MIP) formulation of the location-based model for PPDSP and prove its correctness as well as the space complexity 5 of the generated problem.

Correctness
The objective function in Eq. (1) maximizes the profit-cost function for all possible delivery routing solutions. Eq. (2) introduces two types of Boolean variables x t od and y t r . Eq. (3) guarantees that, each request can be assigned to at most one truck. Proof. According to Eq. (4) (resp. Eq. (5)), if request r is assigned to truck t, then truck t reaches the pickup (resp. dropoff) point of r at least once via a point that is not the pickup (resp. dropoff) point of r. Eq. (6) ensures that the number of visits and the number of departures of a truck at any location must be equal. Eq. (7) restricts any truck departing from a location to at most one other location. It is clear from the above that, all nodes of r∈Dt {f (r), g(r)} are ensured to be visited and departed by truck t exactly once.
Eq. (8), a canonical MTZ subtour elimination constraint Miller et al. (1960), resrticts the integer variables u t v , where v ∈ V \ {0}, such that they form an ascending series to represent the order that truck t arrives at each location v. Eq. (12) gives the domain of such variables u t v . Furthermore, Eq. (8) associated with Eqs. (4)-(7) also enforces that depot 0 must be visited and departed by truck t exactly once if D t = ∅.
Eq. (9) ensures that, for any request r, the arrival of its pickup point must precede the arrival of its dropoff point by truck t if y t r = 1. Theorem 2. For PPDSP, the capacity constraint can be guaranteed by the Eqs. (10) and (11). (10), whose domains are between 0 and c t (i.e., the capacity of truck t) as given in Eq. (11). Such a variable h t v can be considered as the loading volume of truck t when it departs at location v. We define the amount of change in the loading of truck t at location v as Γ, which equals the loaded amount at v (i.e., r∈f −1 (v) (q r · y t r )) minus the unloaded amount at v (i.e., r∈g −1 (v) (q r · y t r )). The main part of Eq. (10) guarantees that h t d − h t o is exactly equal to Γ when x t od = 1, which is written as a big-M linear inequality, where we specify M as c t + r∈R q r .

Space Complexity
Consider that the number of trucks is m (i.e., |T | = m) and the number of requests is n (i.e., |R| = n). Assume that the average repetition rate of the same locations is k, where k ≥ 1, then we can estimate the number of unique locations as 2nk −1 , which is also the number of integrated nodes (viz., |V \ {0}| = 2nk −1 ). Proof. The number of Boolean variables x t od (resp. y t r ) required in proposed formulation is m(1 + 2nk −1 ) 2 (resp. mn); while both the number of required integer variables u t d and h t d are m(1 + 2nk −1 ). In Table 1, we list the bounded numbers of linear (in)equations that correspond to Eqs. (3)-(10) involved in the proposed PPDSP formulation. Therefore, the number of required variable as well as the number of corresponded linear (in)equations are of Θ(mn 2 k −2 ).

Experiment
We compare the performance of a MIP optimizer in solving randomly generated PPDSP instances based on the proposed formulation (i.e., location-based model), and based on the existing formulation (i.e., request-based model), respectively. Please refer to Appendix A for the specific description of the existing formulation used for the comparative experiments. 7

Instances Generation
The directed graph informations for generating instances are set based on the samples of TSPLIB benchmark with displaying data as the coordinates of nodes. 2 We denote the number of nodes contained in the selected sample as |V |, and set the first of these nodes to be the depot node. The pickup and dropoff points of each request are chosen from the non-depot nodes. Therefore, given an average repetition rate k for the same locations, we need to repeatably select n pairs of non-depot nodes (i.e., a total of 2n repeatable non-depot nodes) to correspond to n requests, where n = round( k(|V |−1) 2 ). In Algorithm 1, we construct the repeaTimeList of length |V | − 1 to mark the number of times of each non-depot node being selected. In order to ensure that each non-depot node is selected at least once, each value on such list is initialized to 1. We continuously generate a random index of repeaTimeList and add one to the value corresponding to the generated index (i.e., the number of times of the non-depot node being selected increases by one) until the sum of these numbers reaches 2n. The pickup and dropoff nodes for each request are randomly paired up in Algorithm 2. We first rewrite repeaTimeList as a list of length 2n, shuffList, which consists of all selected non-depot nodes, where the number of repetitions indicates the number of times they are selected. For example, we have shuffList = [0, 0, 0, 1, 2, 2] for repeaTimeList = [3, 1, 2]. Next, we shuffle shuffList and clear pairList, which is used to store pairs of nodes indicating the pickup and dropoff points of all requests. We then divide shuffList into n pairs in order.
2 http://comopt.ifi.uni-heidelberg.de/software/TSPLIB95/tsp/ 8 If two nodes of any pair are identical (i.e., the pickup and dropoff are the same point), or the pair with considering order is already contained in pairList, then we are back to Line 4. Otherwise, we append the eligible pair of nodes into pairList until the last pair is added.
Algorithm 2: Randomly pair up the pickup and dropoff nodes for each request if i = n − 1 then reshuffle ← false 12 return pairList Since we intend to generate instances corresponding to different k for each selected sample of TSPLIB benchmark, and the number of requests n gets smaller as k decreases, we need to produce a decrementable list of node pairs such that we can remove some of the node pairs to correspond to smaller k, but the remaining part of the node pairs still contains all locations in the sample other than that as depot (i.e., they still need to be selected at least once). Therefore, in Algorithm 3, we assume that a sorted list of node pairs pops the end elements out of it one by one, yet it is always guaranteed that every non-depot node is selected at least once. We insert the node pair in which both nodes are currently selected only once into the frontmost of head in Line 7, and then append the node pair in which one of them is selected only once into head in Line 11. Then we keep inserting the node pair in which the sum of the repetitions of the two nodes is currently maximum into the frontmost of tail in Line 22. Note that such insertions and appends require popping the corresponding node pair out of pairList and updating L. We merge head and tail to obtain a sorted list of node pairs sortedPairs at the end of Algorithm 3. Finally, we can generate a list containing n requests for each selected sample of TSPLIB benchmark, as shown in Algorithm 4.
In addition, we generate data for the three types of trucks recursively in the order of maximum load capacity of 25, 20, and 15 until the number of generated trucks reaches m. And these three types of trucks correspond to their respective cost coefficients for each traversing arc of 1.2, 1, and 0.8. For example, for a truck t of load capacity 15, l t od , i.e., the Algorithm 3: Sort the node pairs by the sum of the repeat times of each node in the node pair Input: repeaTimeList, pairList We generate a total of 3 × 5 × 5 = 75 PPDSP instances for our comparative experiments based on the proposed formulation and on the existing formulation, respectively, according to the following parameter settings.
• The selected TSPLIB samples are {burma14, ulysses16, ulysses22 }; 3 Distance( > od) is the Euclidean distance between the coordinates of node o and the coordinates of node d.

Experimental Settings
All experiments are performed on an Apple M1 Pro chip, using the Ubuntu 18.04.6 LTS operating system via the Podman virtual machine with 19 GB of allocated memory. A single CPU core is used for each experiment. We implemented the instance generators for both the proposed formulation and the existing formulation by using Python 3. Each generated problem instance is solved by the MIP optimizer-Cplex of version 20.1.0.0 CPLEX (2020) within 3,600 CPU seconds time limit. Source code for our experiments is available at https://github.com/ReprodSuplem/PPDSP.

Results
Tables 2-4 show the performance of the existing formulation-based method and our proposed formulation-based method in terms of the number of generated variables (#Var.), the number of generated constraints (#Con.), and the optimal values (Opt.) for the various average repetition rates of the same locations (k) and the different numbers of trucks (m), corresponding to TSPLIB samples burma14, ulysses16 and ulysses22, respectively. Each cell recording the left and right values corresponds to a comparison item, where the left value refers to the performance of the item based on the existing formulation-based method, while the right value corresponds to the performance of the item based on the proposed formulation-based method. We compare the performance of these two methods in such cells and put the values of the dominant side in bold.
We can see that either the number of variables or the number of constraints generated by our proposed method is proportional to m, while neither the number of variables nor the number of constraints generated by our proposed method increases significantly as k becomes larger. Such a result is consistent with Theorem 3. In contrast, although the number of variables and the number of constraints produced by the existing method is Table 2: Comparison of the existing formulation-based method with the proposed formulation-based method for TSPLIB sample burma14 in terms of the number of generated variables, the number of generated constraints, and the optimal values. m item k = 1 k = 1.5 k = 2 k = 2.5 k = 3 (n = 7) (n = 10) (n = 13) (n = 16) (n = 20) also proportional to m, with k becoming larger, exponential increases in both of them are observed. Furthermore, it is interesting to note that even when k = 1, the proposed method generates fewer variables than the existing method. As k increases, the optimal value obtained by the proposed method increases within the time limit; while there is no such trend in the optimal values obtained by the existing method. In addition, as m grows up, in theory, the upper bound of the optimal value cannot be smaller. However, the optimal values obtained by the proposed method and the existing method do not maintain the trend of increasing, instead, they both have inflection points. This is due to the fact that the search space of the problem becomes huge, which makes it inefficient for the solver to update the optimized solution. There are even cases on the existing method where the initial solution is not obtained until the end of time.
Last but not least, we can clearly see that for those problems with larger k (i.e., more nodes that can be integrated), the method based on our proposed formulation generates fewer variables and constraints, as well as achieves larger optimal values, than the method based on the existing formulation.

Conclusion
In this paper, we revisit PPDSP on road networks with the integratable nodes. For such application scenarios, we define a location-based graph-theoretic model and give the corresponding MIP formulation. We prove the correctness of this formulation as well as analyze its space complexity. We compare the proposed method with the existing formulation-based method. The experimental results show that, for the instances with more integratable nodes, our method has a significant advantage over the existing method in terms of the generated problem size, and the optimized values.

A. Existing Formulation
Here we provide a detailed description of the existing formulation proposed in Qiu et al. (2017) that were used to perform the comparative experiments in Section 4.
As we stated before, in order to compare the most essential differences between the request-based model and the location-based model, the constraints on time windows are not considered in this study. However, because in the existing formulation, the time windows constraints also act as the subtour elimination constraint. For this reason, we replace those constraints on time windows in the existing formulation with the same MTZ subtour elimination constraint as in our proposed formulation to ensure the fairness of the comparative experiments.
In the existing formulation, n requests are represented as a directed graph G = (V, E), whose node set V = {0} ∪ P ∪ D ∪ {2n + 1}, where P = {1, 2, . . . , n} is the set of pickup nodes, D = {n + 1, n + 2, . . . , 2n} is the set of dropoff nodes, and node 0 and 2n + 1 indicate two depots. 4 A request r is formalized by using the combination (r, r + n) to refer to the pickup and dropoff nodes, and each node v (v ∈ P ∪ D) corresponds to a change amount q v of the truck's load, where q r and q r+n respectively correspond to the loading and unloading of request r, q r > 0, q r+n < 0, q r = |q r+n | and r ∈ {1, 2, . . . , n}. In addition, we particularly specify q 0 = q 2n+1 = 0 for the two depot nodes. The definitions and notations of the other variables are consistent with that previously described in this paper. The existing formulation used for comparative experiments are as follows.
It should be noted that we found three incorrectnesses in the existing formulations in Qiu et al. (2017), which correspond to Eqs. (A.3), (A.4) and (A.12) in the above formulation, and we have rectified each of them.