Ismael Bada

Juan Karsten

Dominique Fohr

Irina Illina

Data Selection in the Framework of Automatic Speech Recognition

Keywords: speech recognition, neural networks, acoustic model, data selection

Training a speech recognition system needs audio data and their corresponding exact transcriptions. However, manual transcribing is expensive, labor intensive and error-prone. Some sources, such as TV broadcast, have subtitles. Subtitles are closed to the exact transcription, but not exactly the same. Some sentences might be paraphrased, deleted, changed in word order, etc. Building automatic speech recognition from inexact subtitles may result in a poor models and low performance system. Therefore, selecting data is crucial to obtain a highly performance models. In this work, we explore the lightly supervised approach, which is a process to select a good acoustic data to train Deep Neural Network acoustic models. We study data selection methods based on phone matched error rate and average word duration. Furthermore, we propose a new data selection method combining three recognizers. Recognizing the development set produces word error rate that is the metric to measure how good the model is. Data selection methods are evaluated on the real TV broadcast dataset.

Introduction

Automatic speech recognition (ASR) is one of the sub field of natural language processing with many practical applications: automatic closed captioning for hearing-disabled persons, taking notes of conversations between doctors and patients, voice control and many more. Despite of the rapid development of speech recognition, there are still many challenges in the field. One of the challenges is the training of a speech recognizer, which requires a huge amount of transcribed training data. The transcribed training data consists of audio data and the corresponding text transcriptions. However, transcribing audio manually is labor intensive and also time consuming. There exist many unlimited supply of audio data from internet, TV broadcasts, radio, as well as video streaming websites, but there is no available exact transcription. However, some TV broadcasts, such as CNN headline news, ABC world news tonight, BBC, have subtitles that can be used for training a speech recognition system.

To train a speech recognition system, one possibility is to use TV broadcasts data that have subtitles. These subtitles are close, but not exactly the same as what people uttered. Some sentences might be paraphrased, deleted, changed in word order, etc. There are some examples of approximate subtitles:

• Real transcription: Russia started badly with the dropping at the hands of Spain. But, they got better and better. Spain looked unstoppable to start with but since then they have looked a little.

•

Corresponding subtitle: Russia started badly with at beating at the hands of Spain. Spain looked then they have looked a little. Furthermore, subtitles are often badly aligned with the audio. Some segments in training audio can contain unconstrained conversational speech, use of foreign words, high out-ofvocabulary rates, channel noise and simultaneous speech from more than one speaker. Even, thus audio data is sometimes difficult to be recognized by humans. These facts make hard to use subtitles for ASR.

The idea of using untranscribed audio data (or unsupervised training, no subtitles) has been proposed firstly in [START_REF] Zavaliagkos | Utilizing Untranscribed Training Data to Improve Performance[END_REF] and [START_REF] Kemp | Unsupervised Training of a Speech Recognizer: Recent Experiments[END_REF]. Authors of [START_REF] Zavaliagkos | Utilizing Untranscribed Training Data to Improve Performance[END_REF] proposed an iterative training procedure: decode untranscribed data and keep only the segments with high confidence score for the next training iteration. Even an 80% error rate system can improve itself automatically, but the system performance is limited. [START_REF] Lamel | Lightly Super-vised and Unsupervised Acoustic Model Training[END_REF] were the first to propose lightly supervised training with a large amount of training data. Instead of using untranscribed training data, they trained speech recognition system using audio data with subtitles. Lightly supervised approach allows selecting "good" training data. First, an acoustic model from another task (or another corpus) is used to recognize audio data. The decoding results are compared with the subtitles and removed if they disagree. These selected data are used to train a new acoustic model. [START_REF] Chan | Improving broadcast news transcription by lightly supervised discriminative training[END_REF] proposed the confidence measure metric to remove the bad audio segments. When decoding acoustic inputs, an ASR produces word hypothesis and their corresponding confidence measure. The confidence measure value is used to remove potentially bad segments where the confident value is lower than a threshold. [START_REF] Mathias | Discriminative training of acoustic models applied to domains with unreliable transcripts[END_REF] applied lightly supervised approach on medical conversation data.

Very recently, a new point of view on the data selection has been proposed. [START_REF] Levit | Don't Count on ASR to Transcribe for You: Breaking Bias with Two Crowds[END_REF] suggest an original two-stage crowdsourcing alternative. First, iteratively collects transcription hypotheses from the web and, then, asks different crowds to pick the best of them. [START_REF] Li | Large-Scale Domain Adaptation via Teacher-Student Learning[END_REF] proposed an approach to domain adaptation that does not require transcriptions but instead uses a corpus of unlabeled parallel data, consisting of pairs of samples from the source domain of the well-trained model and the desired target domain.

In the present paper, the same problem of data selection for acoustic modeling training using a huge data corpus is considered. We want to select a good acoustic data to train Deep Neural Network acoustic models. The scientific contributions of this paper are:

-We study the impact of data selection on the word error rate.

-

We explore different variations of slightly supervised training of acoustic models.

-We present a comparison of different data selection approaches in the context of TV broadcast news speech transcription.

Methodology

Lightly supervised data selection

To generate an accurate speech recognition, a very large training audio corpus with its exact corresponding transcription is required. This is particularly true for Deep Neural Network (DNN) based systems, having millions of parameters to train. However, transcribing audio is labor intensive and time consuming. There are unlimited supply of audio data in the internet, television, radio and other sources. But very few have available transcription. However, some TV broadcasts have subtitles. By utilizing these audio data with the corresponding subtitles, we hope to produce a high performance speech recognizer with less supervision. Nevertheless, some problems exist when using the data with subtitles as training dataset.

Training using the subtitles faces several disadvantages compared to the manual transcriptions: indication of non-speech events (coughing, speaker turn) and acoustic conditions (background noise, music, etc.) are missing.

The main idea of lightly supervised approach is to use the automatic speech recognizer to transcribe training audio data. After this, only well transcribed segments (segments where automatic transcription corresponds to subtitles) will be used as training data [START_REF] Lamel | Lightly Super-vised and Unsupervised Acoustic Model Training[END_REF].

We assume that we have a massive amount of training audio data and corresponding subtitles. In general, the lightly supervised approach operates as follow:

1. Randomly select a subset of the training set.

2.

Train an acoustic model on a small amount of manually annotated data or use an acoustic model from another task. 3. Using ASR, recognize all training audio data. 4. Align the automatic transcriptions with the subtitles of the training data. Some transcriptions and subtitles might disagree. We can remove or correct these segments.

5.

Retrain a new acoustic model using the data we selected in the previous step. 6. Optionally reiterate from step 3. These steps can be iterated several times as long as the error rate is decreasing. This method uses the idea of training acoustic models in less supervised manner because the training dataset (subtitles) is not the actual transcription. Using subtitles as training data greatly reduces the manual transcription effort (20-40 time less).

Revisited lightly supervised data selection

In the lightly supervised approach presented previously, a very important step is the step 3. In the case of a disagreement between automatic transcriptions and subtitles, which part of subtitles to keep and which part to remove or correct? Can we use additional criteria to better choose the training data? How many training data to keep? In this section, we propose to study some of these questions.

Using AWD and PMER

According to [START_REF] Lanchantin | Selection of multi-genre broadcast data for the training of automatic speech recognition systems[END_REF], using of Average Word Duration (AWD) and Phone Matched Error Rate (PMER) during the data selection step (step 3) allows increasing greatly the quality of the selected training data. AWD is used as metric to detect if errors occur in aligning the start and end time of a segment or if something went wrong in the recognition process.

= ℎ

Usually, duration of a word cannot exceed an upper limit threshold and the duration cannot be lower than a bottom limit threshold. If it is the case, this means that the corresponding transcription or subtitle is wrong.

Phone Error Rate (PER) and Word Error Rate (WER) are usually used to measure the performance of a speech recognition system:

= 100 * + + _ _
Error rate is obtained by comparing exact transcriptions and decoding transcriptions produced by the speech recognition system. Word error rate is obtained by the comparison at the word level, phone error rate at the phone level. Our training set has only subtitles. So we can only compare subtitles and recognized transcriptions. To avoid the confusion, we will use Phone Matched Error Rate (PMER). High PMER shows that at phone level the corresponding subtitle is very different compared to recognized transcription. This means possible problems in audio signal (noise, music) or in subtitle. In this case, it is better to discard this segment from the training set.

We chose to use PMER and not WMER because we use phone acoustic models. During acoustic training we interested by the phone sequence and not by the word sequence. For example, the words "too" and "two" have the same sequence of phones: / t uw /. If we misrecognized "two" instead of "too", it will be sad to reject corresponding subtitle since these two words have a same phone sequence.

In our work, we propose to use these measures to increase the quality of the data selection. The proposed iterative methodology is as follow:

1. Randomly select a subset of the training set. This set is used to train an initial acoustic model. Choose N hours of the top PMER segments to make a new training set to train the next acoustic model. 6. Continue the step 2-5 until the data selection does not improve anymore. At each iteration, the number N of selected hours can be augmented. To measure the improvement of the approach at each iteration, a development set recognition could be performed.

System combination

Usually, different ASR systems (with different acoustic models and/or language models) will make different errors. Thus, if several systems provide the same transcription as the original subtitle for one segment, it is very likely that the subtitle corresponds exactly to what has been uttered. We can use it reliably for training acoustic models. The general idea of system combination approach is to combine different ASR systems by varying the language models or acoustic models or both. We have chosen to vary the language model because an acoustic model variation is a very time consuming task when we use a huge training set. The language models can be built with different constraints. A constrained language model is trained only with the sentences of the training corpus used for selection. A less constrained LM is trained also on data from different sources. The idea is as follows: if the recognition result of one ASR and the recognition result of another ASR are the same, we can trust this recognition result. The proposed combination approach works almost in the same way as the method of section 2.3: training acoustic model with the subset of training data (audio data and their subtitles), recognizing and selecting from the full training data and repeating these steps to do better data selection. However, the difference lies in the recognition and data selection steps.

We built three speech recognition systems and each recognition system is used to perform recognition of the same training set. Consequently, we had three transcriptions which have the same amount of segments. We average the value of PMER and AWD from three corresponding decoding transcriptions. After this, we select the training data (step 4 of the approach presented in section 2.2.1) with the proposed combination algorithm:

Select the subtitles from the training set based on AWD as follow: threshold1 < AWD < threshold2 If a segment has zero PMER with one ASR, select the segment and corresponding subtitles Else If a segment have the same phone sequence using two ASRs and PMER< thresholdPMER, select the segment and the corresponding subtitle.

Else sort by average PMER. Choose top N hours segments with the lowest PMER. These subtitles will be chosen to train the next acoustic model.

We hope that using the recognition results when two ASRs agree will help. If a development corpus is available, the thresholds can be chosen to minimize word error rate.

Experiments

Audio corpus

We used the data from the Multi Genre Broadcast (MGB) challenge [START_REF] Bell | The MGB Challenge: Evaluating Multi-Genre Broadcast Media Recognition[END_REF], [START_REF] Woodland | Cambridge University Transcription Systems for the Multigenre Broadcast Challenge[END_REF]. MGB is a challenge to automatically transcribe TV broadcasts. TV broadcast data are recorded in highly diverse environments, speech with background music, non-speech events and sounds, etc. The challenge organizers only provided TV broadcast audio data and their corresponding subtitles. As presented previously, subtitles may be different compared to the actual transcription due to deletion, insertion, substitution and paraphrasing. Thus, MGB data recognition is a very difficult task.

MGB challenge data consists of: Table 1 shows the statistics of the training and the development sets. We can see that, in average, each show contains about 2/3 of speech and 1/3 of non-speech events. These non-speech events are difficult to recognize.

Transcription system

KATS (Kaldi-based Automatic Transcription System) speech recognition system is based on Context Dependent HMM-TDNN phone models [START_REF] Peddinti | A time delay neural network architecture for efficient modeling of long temporal contexts Proceedings of Interspeech[END_REF] [2] [START_REF] Deng | Recent Advances in Deep Learning for Speech Research at Microsoft[END_REF]. We used Kaldi toolkit for training and for recognition [START_REF] Povey | The kaldi speech recognition toolkit[END_REF]. The TDNN architecture has 6 hidden layers, each hidden unit utilizes Rectified Linear Unit (RELU) activation function. The TDNN has around 9100 output nodes (senones) with softmax activation function. The feature vectors are MFCC with 40 feature values. The baseline phonetic lexicon contains 118k pronunciations for 112k words. Using the pocolm [START_REF] Povey | [END_REF], 3gram language model is estimated on text corpora of about 640 million words. 1 present some statistics of the training set of MGB. From the figure we observe that the training set has a high number of segments (subtitles) of average duration of 4.96 seconds and about 14.4 words per segment. Figure 1 shows that there is a large number of subtitles with only few words, so they correspond to a very short speech duration. These short segments can be not easy to recognize by an ASR. AWD values were given by the MGB organizers and obtained after ASR recognition. We can see that the majority of speech segments have an AWD between 0.25 seconds and 0.6 seconds. If one segment has a very small AWD or a very high one, it means that something went wrong and this segment corresponds rather to non-speech events. For safety reasons, for data selection we have extended this interval and we chose AWD between 0.16 and 0.6 for the following experiments.

Experimental results

Study of MGB training data

Table 3 presents the number of speech hours in function of PMER on the training set. PMER values for each segment were given by the MGB organizers. We can observe that one third of the training data have PMER greater than 30%. This means that if we want to keep only a very good training data with very low PMER (so, with a very good quality subtitles), we will have a small training set. In contrast, if we keep all data, a large number of subtitles do not correspond to what was uttered and a negative impact on training is observed.

Impact of the data selection

In order to assess the influence of the data selection, we trained different ASR systems with different amount of training data. The amount of training hours is selected according to the PMER provided by the MGB organizers. For these experiments, we kept only training data with the PMER below some threshold. Figure 3 presents the WER on the development set according to the number of hours selected for training the ASR system (see Table 3). For example, for PMER below 50 we kept only 305 hours of corresponding training data. The trained system is used to recognize the development data and the obtained WER is presented in Figure 3. From this figure we can observe that, at first, WER decreases when the amount of training data increases. But selecting too much data (i.e. using subtitles with high PMER) the WER begins to increase. Therefore, it is important to find a compromise between the quantity of training data and the quality of training data. In conclusion, data selection is important to train an efficient ASR.

Results of data selection methods

We studied and evaluated the presented data selection approaches on the development corpus.

The execution time of one iteration of data selection takes about 43 hours. This is very time consuming. To speed up the experiments, at each iteration of data selection, we decided to select a different number of hours of data (parameter N in the selection algorithm). We hope that a strong selection at the first iteration and less constrained selection at the next iterations will improve the recognition results.

To build an initial acoustic model (called ASR-AM0), we selected randomly 100 hours because we do not have any information about the quality of available subtitles (in real life, only subtitles are available with no addition information, neither PMER nor AWD).

According to the algorithm of section 2.2.1, during the first iteration of data selection, we decode the whole training corpus with ASR-AM0 with KATS system. We kept only segments whose AWD is inside [0.16, 0.6]. We excluded all other segments. We sorted the remaining segments according to PMER and we select N=150 hours. With these 150 hours, we trained ASR-AM1. For second iteration N=200h (ASR-AM2) and for the last iteration N=300h (ASR-AM3).

Table 4 presents the recognition results on the development set for each iteration of the data selection algorithm. Results are presented in terms of percent of correct recognition and in term of WER. The best results are highlighted in bold.

System combination

For system combination, we designed three different recognition systems. They share the same acoustic models (ASR-AM3), but the language models are different. For the first one only the subtitles of the training corpus are used to train the LM. This model is the most constrained. For the second one, the LM is trained using 640 million words of TV subtitles provided by the organizers of the MGB challenge (least constrained model). The last one is a combination of the two previous language models. The thresholdPMER was chosen experimentally and its value is 30.

Using these three ASRs for system combination, a relative improvement of 2% on WER was observed compared to ASR-AM3 (cf. table 5). This improvement is significant. It could be interesting to combine systems using different acoustic models, for instance different acoustic features or different neural networks architecture (Long Short Term Memory, Highway networks).

ASR

Conclusion

In this article, we explored different methods of data selection for building an automatic speech recognition system. The methods are inspired by lightly supervised technique. We studied data selection methods based on phone matched error rate and average word duration. Furthermore, we proposed a new data selection method combining three recognizers. The experiments are conducted on a TV broadcast corpus with subtitles. We have shown that selecting data is crucial for obtaining accurate acoustic models. We have studied the influence of PMER on data selection. The proposed system combination is beneficial to select better data and to obtain an efficient acoustic model.

2 .

 2 Train an acoustic model using the audio and the subtitles of this training set. 3. Decode the full training set with the obtained acoustic model. This will produce new decoding results and new values of PMER and AWD. The new values of PMER are obtained from comparing the subtitles and the decoding results. 4. Select the subtitles from the training set based on AWD as follow: threshold1 < AWD < threshold2 5. Sort the obtained segments according to PMER.

Figure 1 :

 1 Figure 1: Histogram of number of words in function of number of subtitles. MGB training set.

Figure 2

 2 Figure 2 displays the histogram of number of hours of speech in function of AWD for the training set.AWD values were given by the MGB organizers and obtained after ASR recognition. We can see that the majority of speech segments have an AWD between 0.25 seconds and 0.6 seconds. If one segment has a very small AWD or a very high one, it means that something went wrong and this segment corresponds rather to non-speech events. For safety reasons, for data selection we have extended this interval and we chose AWD between 0.16 and 0.6 for the following experiments.Table3presents the number of speech hours in function of PMER on the training set. PMER values for each segment were given by the MGB organizers. We can observe that one third of the training data have PMER greater than 30%. This means that if we want to keep only a very good training data with very low PMER (so, with a very good quality subtitles), we will have a small training set. In contrast, if we keep all data, a large number

Figure 2 : 3 :

 23 Figure 2: Distribution of number of hours of speech in function of AWD (in seconds) for the training set PMER <3 <15 <30 <50 <80 All # hours 112 210 260 304 311 349 Duration % 32 60 74 87 89 100Table 3: Number of hours of training speech according to PMER. Duration (%) as percentage of the total train set.

Figure 3 :

 3 Figure 3: WER on the development set according to the number of hours selected for training the ASR system.

Figure 4 :

 4 Figure 4: Number of hours of speech in function of PMER for the train set. For example, to PMER between 5 and 10 corresponds about 36 hours of speech.

Figure 4

 4 Figure 4 gives more details about the training data distribution in function of PMER values and for different acoustic models. Firstly, this figure shows that a large amount of training data have a PMER below 15 and, so, have good quality subtitles. Secondly, PMER of 0 corresponds about 38 hours of speech at initial iteration, 57 hours for ASR-AM1 and about 60 hours for ASR-AM2. This shows that from one iteration to another the acoustic model performance increases and the acoustic model choose better training data.

Table 1 :

 1 1. A training set contains audio data with their corresponding subtitles. This training set is used for training speech recognition systems. 2. A development set contains around 8 hours of audio data and their corresponding manual transcriptions (exact transcriptions). This dataset can be used to evaluate studied approaches. 3. A text corpus: 640 million words of TV subtitles are provided. These data can be used to train ASR language model. MGB challenge datasets.

	Datasets	# of shows	# hours of shows	# hours of speech
	Training	751	470	349
	Development	16	8.8	6.8

Table 2 :

 2 MGB challenge train set statistics.

	Number of segments (subtitles)	253 K
	Average segment duration	4.96 sec
	Average number of words per segment	14.4
	Vocabulary size	52 K
	Total number of words	3 650 K

Table 2 and figure

 2figure

Table 4 :

 4 Table shows that each data selection iteration improves the ASR system. The best result of 35% WER is obtained at the last iteration. Performing one more iteration gives the same result and is not presented in the table. Results presented in table 4 are not comparable with those in figure 3 because in figure 3, we used PMER given by the organizers. WER on development set for different data selection iterations. Language model is estimated on text corpora of about 640 million words.

	ASR	#hours selected	PMER	Corr (%)	WER (%)
	ASR-AM0	100	10	65.2	40.2
	ASR-AM1	150	15	69.0	36.1
	ASR-AM2	200	21	69.3	35.7
	ASR-AM3	300	49	69.7	35.0

Table 5 :

 5 WER on development set.

		#hours	Corr	(Sub, Del, Ins)	WER
		selected	(%)	(%)	(%)
	ASR-AM3	300	69.7	(14.3, 16.0, 4.7)	35.0
	System combination	300	70.2	(13.8, 16.0, 4.5)	34.3

Acknowledgements

This work is funded by the CPER LCNH project supported by the Great East region and by AMIS (Access Multilingual Information opinionS) Chist-Era project. Experiments presented in this paper were carried out using the Grid'5000 testbed, supported by a scientific interest group hosted by Inria, CNRS, RENATER and other Universities and organizations (https://www.grid5000.fr).


Xavier Nicolay
email: xavier.nicolay@univ-reunion.fr

Rehan Noordally
email: rehan.noordally@univ-reunion.fr

Pascal Anelli
email: pascal.anelli@univ-reunion.fr

Nour Mohammad Murad

Tahiry Razafindralambo
email: tahiry.razafindralambo@univ-reunion.fr

Where is My Next Hop ? The Case of Indian Ocean Islands

Keywords: Latency, Metrology, Active measurement, End-to-end delay, Mapping, Peering, Route, Network

Internet has become a foundation of our modern society. However, all regions or countries do not have the same Internet access regarding quality especially in the Indian Ocean Area (IOA). To improve this quality it is important to have a deep knowledge of the Internet physical and logical topology and associated performance. However, these knowledges are not shared by Internet service providers. In this paper, we describe a large scale measurement study in which we deploy probes in different IOA countries, we generate network traces, develop a tool to extract useful information and analyze these information. We show that most of the IOA traffic exits through one point even if there exists multiple exit points.

Introduction

Internet is fundamental to our society as it provides important services ranging from safety and security services to entertainment. Internet was designed to carry applications data with no time constraint and limited user interactions. Nowadays, several applications focus on users' interactions and timely content delivery is critical. As the Internet connectivity is expected to improve, the user experience is also expected to improve. However, all countries or regions are not equal from the Internet access and performance point of view.

The Internet end-to-end connectivity is not improving, as demonstrated by Lee et al. in [START_REF] Lee | Has internet delay gotten better or worse?[END_REF] and Cardozo et al. in [START_REF] Cardozo | On the endto-end connectivity evolution of the Internet[END_REF]. Within the last decade, changes in the Internet path and bufferbloat issues have worsen the TCP's congestion control [START_REF] Gettys | Bufferbloat: Dark buffers in the Internet[END_REF]. However, beyond bandwidth, low latency is required for new Internet applications. Some geographic area such as The Indian Ocean Area (IOA) including Madagascar, Mauritius, Mayotte, Reunion Island and Seychelles have poorly meshed Internet topology and low performance. In [START_REF] Noordally | Analysis of internet latency: The reunion island case[END_REF][START_REF] Noordally | How long delays impact tcp performance for a connectivity from reunion island ?[END_REF], Noordally et al. focus on Reunion Island Internet connectivity and performance. In this paper, we focus on the whole Indian Ocean Area (IOA) region.

Improving the Internet access for the IOA is important since it can break the Internet silo in this region and help the development of these countries. Briscoe et al. wrote a survey in [START_REF] Briscoe | Reducing internet latency: A survey of techniques and their merits[END_REF] that describes the factors of latency and proposes some solutions to reduce it. It includes exploiting path diversity to select the shortest path and load-balancing to prevent congestion. To evaluate the possible implementation of these solutions, we need a deep understanding of the physical and logical Internet topology in this area.

To the best of our knowledge, no such large scale study has been carried out by the scientific community in the IOA region. This paper tries to fill this gap by:

• Deploying 16 probes in different countries located in IOA: Madagascar (MG), Mauritius (MU), Reunion Island (RE), Seychelles (SC) and Mayotte (YT).

• Generating 4, 480, 000 traceroute traces using randomly selected IPv4 addresses during a one month measurement campaign.

• Developing a tool to extract the logical topology of the Internet in the IOA region based on IP localization and ICMP variant of Paris-Traceroute.

• Analyzing the traces and provide some insight on issues of Internet Access in the IOA regarding exit point of each country, path length and geographical distances.

Our main contribution is to analyze and understand the traffic information from the IOA islands to identify the bottleneck of the Internet traffic in this region.

The remainder of this paper is organized as follows. Section 2 describes the topology of the submarine cables connecting Indian Ocean's Island to the Internet as well at the Internet eXchange Point (IXP). The results are analyzed in section 5. Section 6 reviews the related work. Finally, we conclude in section 7.

Background

The map in figure 1 1 shows that each Island is connected to the Internet with one or more submarine cables. We can notice that LION/LION2 cable provides a link between Mayotte, Madagascar, Reunion and Mauritius. These 4 islands have also an IXP. This equipment can be used by each Internet Service Provider (ISP) connected to it to exchange their traffic.

We know that:

• Each Country / Island is connected to one or more Submarines cables (see figure 1), so we know the real topology.

• There is 4 IXP (Mayotte, Reunion, Mauritius, Madagascar). See AXIS Project [START_REF]African internet exchange system project[END_REF] for more information.

• There are a few of ISP, none are present everywhere: ComoresTelecom, Emtel, CEB FiberNET Co Ltd, Blueline, Telma, Canal + Telecom, Orange, SRR (SFR Réunion), STOI, Telco OI (Only), Zeop, Airtel, Cable & Wireless, Intelvision and Kokonet.

We do not know:

• The exact interconnection between the IXPs. This information is very useful to evaluate logical topology of the network.

• The logical path of a TCP/IP session. This information is related to the core objective of our paper. We aim at analyzing the Internet access of IOA islands.

• The regional traffic in percent of international traffic. This information could help us in our analysis. Indeed, Internet access performance is strongly correlated with traffic shapes.

• The capacity of each Internet Service Provider. This information could provide us some intuition on the traffic and peering policy of each ISP.

In this paper, based on our knowledge of the IOA islands Internet architecture, we aim at analyzing and understanding the traffic information from the IOA islands to identify the bottleneck of the Internet traffic in this region.

Measurement operations

We study the Internet connectivity of Islands localized in the Indian Ocean Area according to the delay and the network paths. To do so we collected traceroute traces between some islands of Indian Ocean and 10, 000 destinations distributed worldwide.

Our active measurements made from Indian Ocean involve 16 raspberry pi [START_REF]Raspberry Pi, official website[END_REF] probes distributed over the 5 countries: 2 hosted at Madagascar, 1 at Mauritius, 1 at Mayotte, 11 at Reunion Island, and 1 at Seychelles. Our trace includes measurements performed from March 22 nd 2017 to April 22 nd 2017. We created a random set of 1, 000, 000 public IPv4 addresses among which only 83, 850 responded to ICMP Echo request.

This new set was geo-referenced by country. The second column of table 1 shows the geographical distribution of these IPv4 addresses and the third column shows the actual distribution of the IPv4 addresses provided by the website https://www.countryipblocks.net 2 . The two distributions are distant from one another. To respect the actual representation, we have decided to use the second one. Among these 83, 850 IPv4 addresses, we selected a subset of 10, 000 addresses that fits the actual geographical distribution. Each of our local probes was configured to perform a traceroute toward all of the IP of our data set within one day. A probe started a new measurement every 8.64s which lasted for an average of 28s. The number of traceroutes running simultaneously has been limited to 4, resulting in a maximum bit-rate of 5, 06Kb/s, which is negligible compared to the available bandwidth which is at least of 128.33Kb/s in Reunion Island [START_REF] Vergoz | Classement des opérateurs Internet réunionnais[END_REF]. To further prevent the congestion induced by our measurements on the destination, the sequence of destinations to visit was randomized on each probe. Our final data set contains a total of 4, 480, 000 traceroute traces.

Tools involved

Traceroute tool

The original traceroute [START_REF] Malkin | Traceroute Using an IP Option[END_REF] developed by Malkin is known to produce inconsistent results in the context of load-balancing.

To circumvent this issue, Paris-Traceroute was created by the authors of [START_REF] Augustin | Avoiding traceroute anomalies with paris traceroute[END_REF]. Thus, TCP packets are sent instead of ICMP. In [START_REF] Wenwei | On evaluating the differences of TCP and ICMP in network measurement[END_REF], the authors compared the ICMP and TCP techniques. While they found that in most of the cases the results are similar, when the ratio between the mean Round-Trip Time (RTT) and the minimum RTT tends to be large (beyond 20), the results of the TCP variant tend to be less stable. For this reason, we use the ICMP version of Paris-Traceroute protocol in our experiment.

Geolocation tool

The coordinate of the IPv4 address were obtained with the database of RIPE NCC [START_REF]RIPE NCC[END_REF]. We used their API [START_REF]Ripestat data api[END_REF] to retrieve information such country, latitude, longitude and AS about each the 83, 850 IPv4 addresses and each of the routers found during the traceroute measurements.

In order to update the localization and enhanced the performance, our own MySQL database is used. We have two main tables: one with {IP, Latitude, Longitude, Countries} and one other with {IP, mapx, mapy}.

We found out that some of the IPs were not properly geo-localized. We inferred an approximate geo-localization of the node according to the minimum delay from several probes distributed worldwide. Then, an IP was considered to be part of the same continent as the probes with which it had the closest delay. The database used for the geolocalization are used in our tool.

Our original tool: Rtraceroute

As we want to handle more than 1 million traces, we develop our tool [START_REF] Nicolay | rtraceroute: Réunion traceroute visualisation[END_REF] with C and threadpool [START_REF] Heber | Simple c thread pool[END_REF]. Maximum of parallelization are implemented and now the tool can parse about 4.5 millions in about 1 hour on a computer with 8 Cores (which is a part of an IBM 3650 on VMWare). All traces are read and IP are geolocated. When the country change between 2 IP address, a link is created. Bogon IP ( * or private IP) change nothing. The tool generate two maps from an empty map of the world [START_REF]Carte vierge du monde[END_REF]. The first one draws all links to a new map. The second map is also created: one point of the link is the country studied (for example: Mauritius which has its coordinates {x,y} = {2611,1569} for the map and stored in the MySQL database). Filtering the links with an extremity permit to show the next country-hop of a country. In other words, we can simply discover the real connectivity of a country.

Our tool is available at http://t.univ-reunion.fr/414

Data filtering tool

Our measurement campaign lead to a raw data set of 4, 480, 000 traceroute traces. But some traces are useless and need to be sanitized. We removed traces that met one of the following criteria:

• the destination has not been reached

• there is 3 following stars at the end of the trace

• the presence of ' !N' (network unreachable) or ' !H' (host unreachable) marks due to Paris-Traceroute

• some corrupted trace (exception probe, empty trace, ...)

• the loops (more than 200 hops)

• the presence of IP whose countries are not present in the RIPE NCC database This criterium was only applied for the geographic path analysis performed in section 5.2

We obtain a new data set after filtering of 1, 053, 894 clean traceroute traces.

Our Traces are available at http://t.univ-reunion.fr/411.

Results

Path length, geographical distance and Round Trip Time

The first set of results we obtained from our large scale experiments with the nature of traffic from the IOA islands to the rest of the world in figure 2. This figure shows the distribution of the distance between IOA islands and their random destination in kilometers based IP location. The distribution result confirms the experiments settings given in table 1. The two pics in the figure at 10,000 km and 15,000 km represents Europe plus Asia and North America respectively. In the second set of results, we plot the average path length in numbers of hops depending of the geographical distance, see figure 3. The measure shows that the number of hops is stable depending relatively to the distance except for Madagascar. It is important to notice here that on average even if the distance between the source IP address (an IOA island) is geographically close to the destination IP address (between 0 and 5,000 km) the number of hops remains high on average. A first interpretation of figure 3 is that the number of hops between an IOA island and another IOA island is on average the same as number of hops between an IOA island and an European country. This is even worse in the case of Madagascar where reaching an IOA island needs on average more hops than reaching an European country.

Figure 4 plots the RTT depending on the number of hops. As expected, the RTT increases with the number of hops. However as shown in Figure 5 the RTT is stable depending on the distance. This result shows that from an IOA island point of view, geographically close (small distance in kilometers) are harder to reach in terms of RTT and number of hops.

The results in this subsection tend to show that geographical distances are not related to hop distances from the IOA islands point of view. They indicate that on average when traffic exits from an IOA island wherever the destination is (in terms of geographical position), this destination will be reached in a constant number of hops.

Path Analysis

In Figure 6, we plot the number of exit points for each IOA island. By exit point, we mean the first hop for each island outside its own country. This than 20 exit points. This figure seems to indicate some path diversity for the Internet traffic of each island. In table 2 the details of the exit points of each country are given. This table gives the ISO-3166 country code 3 for the exit point and the percentage of traffic using this exit point. This table confirms the diversity of the exit points, except for Mayotte, but also shows that for most islands in the IOA more than 90% of our data exits through only three countries (written in bold in the tabular). Table 2 also shows some asymmetrical behavior in routing and peering policies. For example, Mauritius (MU) appears in Madagascar's (MG) exit point but Madagascar does not appear in Mauritius'. 7b give the repartition of exit points by country and by continent for Madagascar. In these figure we only plot the exit points that is used for at least 1% of our data. We can see from these figures that only three countries are used as exit points for Madagascar and that more than 94% of the traffic exists through Europe. When combined with the previous results, we can say that Europe has a well meshed network with all other countries since from IOA island point of view the entire world is at constant number of hop on average. These results also explain the decreasing behavior between the distance in kilometers and the number of hops. Since all traffic exists through Europe, reaching IOA island from Madagascar needs more hop count than reaching an European country. Figures 8a and8b show the exit points for Mauritius. For Mauritius, the number of countries holding more that 1% of the traffic is higher that Madagascar's. These figures show that more than 80% of the traffic of Mauritius exit through Europe. However, it is is important to notice that unlike Madagascar, Figures 9a and9b show the results of Reunion Island and we can see that the trend are the same as Madagascar and Mauritius trends. More than 95% of the traffic exit through Europe and especially through France.

Figures 10a and10b show the results for Seychelles. We can see from these figures that more than 90% of the traffic exists through Europe and especially Great Britain.

For Mayotte, none Figure are needed. All our data gone through Europe, and more specifically in France, as showned in Table 2.We can see that Mayotte is a special case of IOA islands. Indeed, Mayotte uses only one exit point which is France. This make the Internet access of Mayotte not resilient/robust in case In the previous results, showed in Figures 6 to 10, we can see that for traffic from IOA islands more than 80% of the traffic use only one Continent exit point. Moreover, even if this exit point is in Europe for Madagascar, Reunion Island, Mauritius, Seychelles and Mayotte they are not all located in the same country which increase the RTT, number of hop and therefore reduce the Internet performance.

Related work

The IOA routing rules increase delays. It is not the unique cases present in the World. Recent studies about routing rules show their impact on the delay. In [START_REF] Zheng | Internet routing policies and round-trip-times[END_REF], the authors work on the notion of Triangle Inequality Violations (TIV) and its impact on the delay. This notion said that the sum of delay between two nodes of a triangle is necessarily higher as the delay between one of the two previous nodes to the last one. If this rule is not respected, it is the case of TIV. One particular TIV is called Boomerang routing. A study of this phenomenon [START_REF] Obar | Internet surveillance and boomerang routing: A call for Canadian network sovereignty[END_REF] as shown that many paths between Canadian ISP's take indirect paths through the USA. This sort of connection was frequent in IOA. Only the presence of the IXP and a real inter-connection of the ISP could resolve this problem.

For [START_REF] Gupta | Peering at the internet's frontier: A first look at isp interconnectivity in africa[END_REF] the main reason for long delay in the African region is due to peering agreements. Despite the numerous IXPs in South Africa or West Africa, some ISPs preferred to inter-connect in an European or Asian IXP. To bypass this rule, AFNIC and private companies, like Google, Akamai, etc... have made some investments in the African continent [START_REF]African internet exchange system project[END_REF]. In [START_REF] Fanou | On the diversity of interdomain routing in africa[END_REF], authors show that new infrastructures have not been correctly used by the African ISP.

These needed to join an IXP based outside of the African continent, and that dependence on submarine cable. Chan and all worked on the impact of failures in submarine cables, in particularly on the SEA-ME-WE-4 [START_REF] Chan | Non-cooperative diagnosis of submarine cable faults[END_REF]. Furthermore adding a new submarine cable or increasing their bandwidth will not reduce latency [START_REF]Reducing Internet Transport Latency[END_REF]25].

Conclusion

Studying path and delay is a very important task in regions where the Internet access is not very fairly distributed. The Indian Ocean Area are connected to the Internet by only one or two submarines cables, depending on the country. From our probes, we used the paris-traceroute tool to create an active metrology measurement, and used our tool rtraceroute to analyze the data produced and enhance the knowledge between topology and logical paths.

Our results shows that the distance have no impact on the path length. For each new adding-node in the path, the increasing delay are variable, depended of the source's country. The surprising results is the decrease of the delay when the distance between the source and destination increases. The major result indicates that most of the islands are world-connected but with a poor regional peering and meshing. It seems that the IXP and regional peering are not really optimized / well configured. We encourage the ISP for a better use of the IXP. We discover that Mayotte is really a special case with only one direct connection to France. We are also surprised to discover that the only interconnection between the five islands {MG, MU, SC, RE, YT} is reduced to: {M G → M U } and {SC → M U }. It is a poor peering.

The future step of our research concern the deployment of probes in the other French overseas department and compare the situation. We leave for the future work the path and delay evolution over time in the IOA. An analysis of the TCP performance of the IOA could also be done, with the help of the different local ISP. From this study, we can imagine to place a closest Regional IXP, to improve the regional peering and enhance the TCP performance.

Figure 1 :

 1 Figure 1: Mascarene islands submarine cable. [8]

Figure 2 :

 2 Figure 2: Distribution of the distance between IOA's sources and world's destinations.

Figure 3 :

 3 Figure 3: Dependence of the path length on geographical location.

Figure 4 :Figure 5 :

 45 Figure 4: Dependence of the RTT and the path length.

Figure 6 :

 6 Figure 6: Number of Internet output links from each countries from IOA

Figure 7 :

 7 Figure 7: Number of Madagascar's Internet output links (up to 1%)

Figure 8 :

 8 Figure 8: Number of Mauritius's Internet output links (up to 1%)

Figure 9 :

 9 Figure 9: Number of Reunion's Internet output links (up to 1%)

Figure 10 :

 10 Figure 10: Number of Seychelles's Internet output links (up to 1%)

Table 1 :

 1 Geographical distribution of our 83, 850 randomly obtained IPv4 addresses and the actual geographical distribution provided by the CountryIp-Blocks website.

	Continent	Random CountryIPBlocks
	Africa (AF)	0.95%	2.59%
	Asia (AS)	32.96%	23.34%
	Europe (EU)	28.99%	20.7%
	North America (NA)	8.89%	47.55%
	Oceanie (OC)	0.7%	1.55%
	South America (SA)	6.30%	4.27%
	Other (bogons)	21.19%	0.0%

Source: www.submarinecablemap.com

https://fr.wikipedia.org/wiki/ISO_3166

Table 2: Table of exit point from IOA MG MU RE SC YT AE 0.0005% AZ 0.001% AE 0.0001% EU 0.0001% FR 100% AU 0.0005% HU 0.001% DK 0.0001% KE 0.0001% FI 0.0005% IN 0.001% TH 0.0001% AZ 0.002% JP 0.0005% JP 0.001% MQ 0.0005% DZ 0.002% MU 0.0005% TW 0.001% RO 0.0005% ID 0.002% SA 0.0005% CA 0.003% AR 0.0007% IL 0.002% TR 0.0005% AE 0.004% TW 0.0007% NA 0.002% AT 0.001% CO 0.004% RS 0.0014% SA 0.002% CH 0.001% CZ 0.004% CR 0.0017% CO 0.01% GR 0.001% RS 0.004% BG 0.0019% PK 0.011% NG 0.001% SE 0.004% SE 0.0019% HK 0.013% HU 0.0015% TH 0.004% AT 0.0020% RO 0.018% IS 0.0015% BG 0.006% HK 0.0030% IT 0.023% MD 0.0015% GR 0.006% JP 0.0049% TH 0.034% PH 0.0015% PL 0.007% SG 0.0073% ZA 0.037% ES 0.002% US 0.007% NL 0.0103% SE 0.038% KE 0.0025% ZM 0.007% ZA 0.0117% VE 0.039% IE 0.003% RU 0.012% IE 0.0121% IE 0.040% IL 0.004% UG 0.013% RU 0.0241% SG 0.043% AZ 0.0045% CR 0.016% DE 0.0628% PL 0.045% CZ 0.005% DE 0.016% ES 0.0690% AE 0.049% ID 0.005% IE 0.042% EU 0.478% RS 0.049% NL 0.006% SA 0.083% US-CO 1.263% CA 0.059% PK 0.0139% GB 0.181% US 3.542% RU 0.11% CO 0.016% US 0.276% GB 39.0897% FR 0.11% HK 0.0174% SG 0.340% FR 55.4101% US-WA 0.12% RO 0.019% ZA 0.574% US 0.16% CR 0.02% FR 1.167% DE 0.17% BG 0.021% KE 1.217% MU 4.41% IT 0.021% MY 13.764% GB 94.4% VE 0.022% IT 32.228% NA 0.04% EU 50.004% SE 0.04% SG 0.042% RS 0.046% CA 0.059% TH 0.060% RU 0.096% ZA 0.125% DE 0.193% EU 0.250% US 0.602% US 4.39% FR 10% GB 83.79%

the most used exit point is not Great Britain.


Sonia Tabti

Julien Rabin

Abderrahim Elmoataz

SYMMETRIC UPWIND SCHEME FOR DISCRETE WEIGHTED TOTAL VARIATION

Keywords: Weighted Total Variation, Upwind scheme, image denoising, isotropy preservation, co-area formula

. We focus on two new symmetric formulations based on the 2 and ∞ norms respectively and propose a dedicated optimization algorithm to solve convex problems based on such TV penalties. We demonstrate the theoretical and practical interest of such formulations for image processing tasks.

INTRODUCTION

Context In this paper, we are interested in the discrete formulation of Total Variation (TV) in image processing and computer vision as a prior for piecewise constant images. This topic has been widely studied since it has been proposed by Rudin-Osher-Fatemi (ROF) [START_REF] Leonid I Rudin | Nonlinear total variation based noise removal algorithms[END_REF] and remains an active research field; as an example, non-local formulations have been proposed in [START_REF] Gilboa | Nonlocal operators with applications to image processing[END_REF], [START_REF] Li | Smoothing the finite differences defining the non-local total variation and application in image restoration[END_REF] and a combination of TV with Non-Local means has been proposed in [START_REF] Sutour | Adaptive regularization of the NL-means: Application to image and video denoising[END_REF], all in order to preserve textures and thin structures. The Total-Generalized-Variation (TGV) studied in [START_REF] Bredies | Total generalized variation[END_REF] and the method in [START_REF] Louchet | Total variation denoising using iterated conditional expectation[END_REF] are variants that reduce the staircasing effect observed with TV.

Total Variation is applied in many imaging and computer vision problems such as segmentation [START_REF] Tony F Chan | Algorithms for finding global minimizers of image segmentation and denoising models[END_REF], inpainting, deconvolution and motion estimation [START_REF] Chambolle | A first-order primaldual algorithm for convex problems with applications to imaging[END_REF], optical flow [START_REF] Ranftl | Non-local total generalized variation for optical flow estimation[END_REF], filtering [START_REF] Gilboa | A total variation spectral framework for scale and texture analysis[END_REF], point cloud processing [START_REF] Hidane | Nonlinear multilayered representation of graph-signals[END_REF] and classification [START_REF] Bresson | Multiclass total variation clustering[END_REF]. Most studies are based on continuous formulations that are then approximated to solve the discrete problem. In this work, we consider directly the discrete setting.

Motivation In [START_REF] Chambolle | An upwind finite-difference method for total variation-based image smoothing[END_REF], a new discrete TV formulation is proposed, referred to in the following as TV u , and inspired by upwind finite difference schemes used in the numerical resolution of Partial Differential Equations. Despite the upwind scheme [START_REF] Chambolle | An upwind finite-difference method for total variation-based image smoothing[END_REF] ability to better preserve isotropy in comparison with the usual ones, Condat shows in [START_REF] Condat | Discrete total variation: New definition and minimization[END_REF] that it is not invariant to contrast inversion. Consequently, some artifacts tends to appear with TV u such as light isolated dots, which are not suitable especially for filtering purpose. Condat also proposes a new dual formulation of the Total Variation, denoted here by TV c , which avoids this issue. However, the price to pay is that this formulation is not explicit. Note in addition that both schemes do not check the discrete co-area formula. This formula, verified by the usual anisotropic TV definition, allows to solve some non-convex problems by convex relaxation such as segmentation, see for example [START_REF] Tony F Chan | Algorithms for finding global minimizers of image segmentation and denoising models[END_REF].

Contributions In this paper, we study two non-local TV formulations on graphs. This type of graph-based formulation was initially proposed in [START_REF] Elmoataz | Nonlocal discrete regularization on weighted graphs: a framework for image and manifold processing[END_REF]. The proposed formulations rely on a symmetric definition combining an upwind scheme with a downwind scheme expressed in an ∞ and 2 norms respectively. We propose an efficient way to express these norms allowing us to solve convex optimization problems with these regularization terms and illustrate their properties on some numerical experiments. These properties are: isotropy preservation, bias reduction in comparison with TV u , better filtering results in comparison with variational approaches of the same type and the versatility of the non-local ∞ formulation that can mimic other formulations.

Outline In Section 2, we first recall the local, non-local and upwind formulations of the Total Variation, then, we introduce the proposed symmetric formulations. In Section 3, the optimization algorithm to solve convex variational problems based on these regularization terms is described. Some comparative experiments in Section 4 illustrate the interest of the proposed framework.

DISCRETE UPWIND FORMULATIONS OF TOTAL VARIATION

Notations

A graph G(V, E, ω) is defined by a set of vertices V connected by edges E ⊆ V × V and a weight function

ω : E → R + . A graph is symmetric, or undirected if ω(v, u) = ω(u, v) ∀(u, v) ∈ E.
In practice, this property is usually satisfied for applications in imagery except in the domain borders because of special conditions such as Neumann or Dirichlet conditions. In the sequel, we consider a scalar image f defined on a regular cartesian grid Ω ⊂ Z 2 , u :

x ∈ Ω → f (x) ∈ R.
The border of the domain is denoted by ∂Ω. The associated n-

connexity graph is G(Ω, E, ω), where: E = {(x, x + t), ∀ x ∈ Ω\∂Ω, t ∈ N n }.
In practice, we use: N 2 = {(1, 0), (0, 1)} or:

N 4 = {(1, 0), (-1, 0), (0, 1), (0, -1)} with ω = 1, or: N 8 = {u = (i, j) s.t. i, j ∈ {0, 1, -1}} with either ω(i, j) = 1/ f (i) -f (j) 2 or ω = 1.
For a field p : E → R, the inner product is defined by: p, q = (v,u)∈E p(v, u)q(v, u) and the canonical norm by: ||p|| = p, p . We also define the composed norm ||.|| 1,k :

||p|| 1,k :=        v∈V u s.t (u,v)∈E |p(v, u)| k 1 k , if 1 k < ∞ v∈V max u s.t (u,v)∈E |p(v, u)|, if k = ∞
(1) which will be useful to unify the different formulations of the total variation studied here. The dual norm is:

||p|| * 1,k = max ||q|| 1,k 1 p, q = ||p|| ∞,k * with 1 k + 1 k * = 1.
(2) We note 1 (respectively 0) the vector fields where each component is equal to 1 (resp. 0). The adjoint of a linear operator K is noted K * , such that: Kp, q = p, K * q . The operators (p) + = max{0, p} and (p) -= min{0, p} apply a thresholding to each component of p. The inequality: p 0, implies that each component of p is non-negative. The characteristic function of a set C is denoted by χ C , and the indicator function 1 C . The proximal operator of an application f is:

Prox f (x) := argmin y 1 2 ||x -y|| + f (y).

Upwind and downwind difference operators

We recall here the definition of difference operators on graphs as proposed in [START_REF] Elmoataz | Nonlocal discrete regularization on weighted graphs: a framework for image and manifold processing[END_REF]. Without loss of generality and for sake of simplicity, we consider in the sequel fully connected graphs E = V × V, where the lack of an edge is represented by a zero weight ω. The centered difference operator of a vectorial function on G:

∇ ω f : v ∈ V → ω(v, u)(f (u) -f (v)) u∈V ∈ R V .
(3) The upwind and downwind uncentered difference operators are then respectively defined as follow:

∇ - ω f : v ∈ V → (∇ ω f (v)) - (upwind) ∇ + ω f : v ∈ V → (∇ ω f (v)) + (downwind) (4)
The adjoint of the difference operator ∇ * ω : R E → R V is defined in the general case by:

∇ * ω p(v) = u∈V ω(u, v)p(u, v) - u∈V ω(v, u)p(v, u).
(5)

Usual formulations of TV and properties

Discrete co-area formula As mentioned in the introduction, an important property satisfied by the continuous formulation of the Total Variation is the co-area formula which allows for instance to establish a link between the notion of perimeter and the TV norm. In practice, it is interesting to check that the discrete formulation of TV also statisfies this property which is written as follows in our framework:

TV(f ) = R TV χ {f t} dt. (6)
TV formulations Most of usual discrete formulations can be written this way:

TV ω|k (.) := ||∇ ω .|| 1,k with ∈ { , +, -}. ( 7
)
If we consider the local definition of the scalar Total Variation on an n-connexity graph G(Ω, E, ω) with ω = 1, we get the following definitions:1

• TV |2 with n = 2 : isotropic Total Variation [START_REF] Leonid I Rudin | Nonlinear total variation based noise removal algorithms[END_REF],

• TV |1 with n = 2 : anisotropic Total Variation,

• TV u = TV - |2 with n = 4: upwind isotropic Total Variation [START_REF] Chambolle | An upwind finite-difference method for total variation-based image smoothing[END_REF].

In the more general case of Non-Local formulations (NL-TV), with any weight function ω, we get:

• TV ω|2 : isotropic Non-Local Total Variation [START_REF] Gilboa | Nonlocal operators with applications to image processing[END_REF][START_REF] Hidane | Nonlinear multilayered representation of graph-signals[END_REF],

• TV ω|k , TV - ω|k and TV + ω|k with k ∈ [1, ∞) [START_REF] Elmoataz | Nonlocal discrete regularization on weighted graphs: a framework for image and manifold processing[END_REF] . Note that, among all previous definitions, only TV |1 and TV ω|1 check the discrete co-area formula [START_REF] Sutour | Adaptive regularization of the NL-means: Application to image and video denoising[END_REF]. Bias of uncentered upwind (and downwind) schemes Next proposition shows that, unlike usual centered formulations TV ω|k , uncentered schemes, that is to say upwind and downwind schemes, are biased except in the case of undirected graphs with k = 1.

Proposition 1 (Contrast inversion). We have the contrastinversion property ∀ k 1 and ∀ ω:

TV - ω|k (-f ) = TV + ω|k (f ) = TV - ω|k (f ). (8)
If k = 1, uncentered schemes verify:

TV ω|1 = 1 2 TV - ω|1 + 1 2 TV + ω|1 (9) TV - ω|1 (f ) = 1 2 TV ω|1 (f ) -1 2 f, ∇ * ω 1 (10)
where the second term is equal to zero in the undirected graph case.

The lack of invariance to contrast inversion is illustrated in [START_REF] Condat | Discrete total variation: New definition and minimization[END_REF] where light pixels have a lower value of TV - |2 than dark pixels, which produces noticeable denoising artifacts.

Symmetric uncentered formulation

Inspired by the equation ( 9) which shows that bias can be removed by defining an uncentered symmetric scheme, we focus on a symmetric definition of Total Variation with the ∞ and the 2 norms. Definition We consider the following formulation:

TV ± ω|k := 1 2 TV - ω|k + 1 2 TV + ω|k for k ∈ {2, ∞}. (11)
This definition has been proposed in [START_REF] Tenbrinck | Solving minimal surface problems on surfaces and point clouds[END_REF] in the case k = ∞ but has never been used in practice. Notice also that the case k = 1 is not interesting as we have, according to equation ( 9): TV ± ω|1 = TV ω|1 . Properties of the proposed formulations The symmetric formulation ( 11) is a semi-norm and verifies the next proposition in the case k = ∞.

Proposition 2 (Co-area [START_REF] Tenbrinck | Solving minimal surface problems on surfaces and point clouds[END_REF]). A constant weight ω(v, u) in each vertex v ∈ V is a sufficient condition for TV ± ω|∞ to verify the discrete co-area formula [START_REF] Sutour | Adaptive regularization of the NL-means: Application to image and video denoising[END_REF].

Using the definitions in section 2.1, and in particular (2), we can show the next proposition: Proposition 3 (Dual formulation of TV ± ω|k [START_REF] Ranftl | Non-local total generalized variation for optical flow estimation[END_REF]).

TV ± ω|k (f ) = max r1 0,r2 0 1 2 ∇ ω f, r 1 -r 2 -χ ||r1|| ∞,k * ≤1 -χ ||r2|| ∞,k * ≤1
where k * verifies 1 k + 1 k * = 1.

APPLICATION AND OPTIMIZATION

Inverse problems

In order to demonstrate the practical interest of the uncentered symmetric formulations proposed, we focus on the resolution of inverse problems as in the ROF model [START_REF] Leonid I Rudin | Nonlinear total variation based noise removal algorithms[END_REF] :

min f 1 2 ||f -g|| 2 + λTV ± ω|k (f ). ( 12
)
where g is an image corrupted with Gaussian noise, f is the image to be estimated and λ the regularization parameter.

Dual formulation

Using proposition 3, we obtain the following primal-dual problem:

min f ∈R V max r=(r1,r2) 1 2 ||f -g|| 2 + Kf, r -χ r∈S k * , (13)
where K = λ 2 ∇ w ; -∇ w so that ||K|| 2 8λ with neighborhood N 2 (see [START_REF] Chambolle | An algorithm for total variation minimization and applications[END_REF]) and ||K|| 2 16λ with neighborhood N 4 . The convex set S k * is defined as follows:

S k * = {(r 1 , r 2 ) : r i ∈ R E + , ||r i || ∞,k * ≤ 1 ∀ i ∈ {1, 2}}.

Algorithm

Problem ( 13) can be solved with the Chambolle-Pock algorithm [START_REF] Chambolle | A first-order primaldual algorithm for convex problems with applications to imaging[END_REF] by following these steps at each iteration t:

           r (t+1) = Proj S k * (r (t) + γ (t) K f (t) ) f (t+1) = Prox τ (t) 2 ||.-g|| 2 (f (t) -τ (t) K * r (t+1) ) f (t+1) = f (t+1) + θ (t) (f (t+1) -f (t) )
where

θ (t) = 1/ 1 + 2ρτ (t) , τ (t+1) = θ (t) τ (t) , γ (t+1) = γ (t) /θ (t) (14)
with these conditions and parameters setting:

τ 0 , γ 0 > 0, θ 0 ∈ [0, 1], f 0 = g, ρ = 0.7/λ, τ 0 = 0.99/||K||, τ (t) γ (t) ||K|| 2 < 1.
The proximal operator is (see for example [START_REF] Chambolle | A first-order primaldual algorithm for convex problems with applications to imaging[END_REF]):

Prox τ 2 ||.-g|| 2 (f ) = 1 1+τ (f + τ g).
Next proposition proves that, in order to perform the projection on S k * , one should perform a projection onto the nonnegative portion of the k * unit ball. Due to lack of space, the proof of this proposition does not appear in this paper.

Proposition 4 (Projection on S k * ). For k * ∈ {1, 2} Proj S k * (r 1 , r 2 ) = (r 1 , r2 ) where ri (v) = Proj ||.|| k * 1 ((r i (v)) + ) ,
For i ∈ {1, 2}. Hence, for k * = 2, we obtain the same result as [START_REF] Chambolle | An upwind finite-difference method for total variation-based image smoothing[END_REF] with: Proj ||.|| 1 (x) = x max{1,||x||} . For k * = 1, the projection on the 1 -ball can be computed with linear complexity (see for example [START_REF] Condat | Fast projection onto the simplex and the l1 ball[END_REF]).

EXPERIMENTS AND DISCUSSION

Study of shape preservation

Figure 2 shows the regularization results of a synthetic image for the various definitions of TV considered previously. Each shape represents an k ball and corresponds to the structural element of the adjoint TV |k * formulation studied here (see e.g. [START_REF] Chan | Recent developments in total variation image restoration[END_REF] for TV |2 ). As expected, the proposed isotropic formulation TV ± |2 provides a satisfying result on the disc shape, like TV u and TV c . In contrast, TV |2 produces some blur in a particular direction, as already reported in [START_REF] Chambolle | An upwind finite-difference method for total variation-based image smoothing[END_REF][START_REF] Condat | Discrete total variation: New definition and minimization[END_REF]. This observation stands for the other symmetric formulation considering the diamond for TV ± |∞ (in comparison with TV |∞ ). Even more interestingly, the result obtained with TV ± ω|∞ for a larger neighborhood (N 8 instead of N 4 ) and a well chosen weight distribution, ω = 1 ||.|| 2 , "mimic" the behavior of TV ± |2 . Last, observe that the contrast is better preserved with TV ± |∞ in comparison with all these approaches.

Denoising application: Bias removal

= 1/ f (u) -f (v) 2 .
The proposed symmetric formulations do not suffer from artefacts, unlike the upwind scheme TV u , which results in a better PSNR espcially with the Non-Local (NL) weight computed computed as in [START_REF] Buades | A non-local algorithm for image denoising[END_REF]. Staircasing effect is reduced by using TV |H , the Huber TV, see for instance [START_REF] Nikolova | Local strong homogeneity of a regularized estimator[END_REF]. Note that the results obtained with TV ± |2 and TV ± |∞ can be improved with a higher value of λ but the PSNR obtained with the NL weight will still be the highest. problem [START_REF] Gilboa | A total variation spectral framework for scale and texture analysis[END_REF]. In fact, as mentioned in Section 2.2, we observe that TV u produces a bias represented by light isolated pixels. The proposed regularization terms as well as TV c , do not suffer from such bias and manage to provide a better PSNR in comparison with the other schemes in most cases.

f TV |2 , N 2 TV |∞ , N 2 TV c , N 4 [15] TV u , N 4 [1] TV - |∞ , N 4 TV ± ω|∞ , N 8 TV ± |2 , N 4 TV ± |∞ , N 4
We can also observe that the result obtained with the Nonlocal formulation is really improved visually and textures are well preserved. Moreover, with Huber TV, the staircasing effect is well reduced. The same benefits are also observed for color images (Figure 3). The vectorial TV definition used is described in [START_REF] Bresson | Fast dual minimization of the vectorial total variation norm and applications to color image processing[END_REF].

CONCLUSION

In this paper, we proposed two uncentered and symmetric discrete TV schemes on weighted graphs with an 2 and an ∞ norm respectively. We also developed an efficient optimization algorithm to adress generic problems that may benefit from such new formulations. These terms are semi-norms with an explicit formula unlike the framework proposed in [START_REF] Condat | Discrete total variation: New definition and minimization[END_REF]. We illustrated in a denoising application that both proposed formulations remove the bias observed with TV u [START_REF] Chambolle | An upwind finite-difference method for total variation-based image smoothing[END_REF]. An extension to color image filtering has also been proposed in the case of the 2 norm with the same advantages as the grayscale version. Our formulations are generic in the sense that the neighborhood N n and the weights ω can be defined arbitrarily and adapted to specific problems. This property allows the ∞ formulation to be able to mimic the properties of other formulations. We also recall that TV ± 1|∞ verifies the co-area formula which allows to solve some non-convex optimization problems.

Figure 1 Fig. 1 .

 11 Figure 1 illustrates the interest of the proposed formulations in the context of a simple application: the denoising

Fig. 2 .Fig. 3 .

 23 Fig. 2. Regularizations of the image f with 4000 iterations and λ = 6.

up to a normalization factor related to the choice of the neighborhood Nn.

The authors thank ANR-14-CE27-0001 GRAPHSIP and FEDER/FSE 2014/2020 (GRAPHSIP project) for funding.