Viet-Anh Tran
email: viet-anh.tran@gipsa-lab.inpg.fr

Gérard Bailly
email: gerard.bailly@gipsa-lab.inpg.fr

Hélène Loevenbruck
email: helene.loevenbruck@gipsa-lab.inpg.fr

Tomoki Toda
email: tomoki@is.naist.jp

Predicting F0 and voicing from NAM-captured whispered speech

Keywords: voice conversion, F 0 estimation, neural network, non-audible murmur, whispered speech

The NAM-to-speech conversion proposed by Toda and colleagues which converts Non-Audible Murmur (NAM) to audible speech by statistical mapping trained using aligned corpora is a very promising technique, but its performance is still insufficient, mainly due to the difficulty in estimating F 0 of the transformed voice from unvoiced speech. In this paper, we propose a method to improve F 0 estimation and voicing decision in a NAM-to-speech conversion system based on Gaussian Mixture Models (GMM) applied to whispered speech. Instead of combining voicing decision and F 0 estimation in a single GMM, a simple feed-forward neural network is used to detect voiced segments in the whisper while a GMM estimates a continuous melodic contour based on training voiced segments. The error rate for the voiced/unvoiced decision of the network is 6.8% compared to 9.2% with the original system. Our proposal benefits also to F 0 estimation error.

Introduction

Speech conveys a wide range of information. Among them, the linguistic content of the message being uttered is of prime importance. However, paralinguistic information such as the speaker's mood, identity or position with respect to what he/she says also plays a crucial part in oral communication [START_REF] Stylianou | Continuous probabilistic transform for voice conversion[END_REF]. Unfortunately, when a speaker murmurs or whispers, this information is degraded.

To solve this problem, Nakajima et al. [START_REF] Nakajima | Non-audible murmur recognition[END_REF] found that acoustic vibrations in the vocal tract can be captured through the soft tissues of the head with a special acoustic sensor called a NAM microphone attached to the surface of the skin, below the ear. Using this stethoscopic microphone to capture non-audible murmur, Toda et al. [START_REF] Toda | NAM-to-Speech Conversion with Gaussian Mixture Models[END_REF] proposed a NAM-to-Speech conversion system based on the GMM model in order to convert "non-audible speech" to ordinary speech. It was shown that this system effectively works but its performance is still insufficient, especially in the naturalness of the converted speech. This is due to the difficulties in F 0 estimation from unvoiced speech. These authors claimed that it is inevitable to improve the performance of NAM-to-Speech systems. Nakagiri et al. [START_REF] Nakagiri | Improving Body Transmitted Unvoiced Speech with Statistical Voice Conversion[END_REF] proposed another system which converts NAM to whisper. F 0 values do not need to be estimated for converted whispered speech because whisper is another type of unvoiced speech, just like NAM, but more intelligible.

Another direction of research consists in using a phonetic pivot by combining speech recognition and synthesis techniques as in the Ouisper project [START_REF] Hueber | Continuous-Speech Phone Recognition from Ultrasound and Optical Images of the Tongue and Lips[END_REF]. By introducing higher linguistic levels, such systems can potentially predict a phonological structure that can be used in speech resynthesis. But no results have been reported yet and such an approach seems unsuitable for applications with open domain.

In this paper, we propose to improve signal-based GMM mapping by a better estimation of the voiced source of the converted speech. Whisper-to-speech was used because of difficulties in getting accurate phonetic segmentation in NAM.

In the training stage, whispered speech and ordinary speech utterance pairs were carefully aligned using phonetic transcription information. The main difference between our system and the original system proposed in [START_REF] Toda | NAM-to-Speech Conversion with Gaussian Mixture Models[END_REF] is that only voiced segments were provided as input to train the GMM model which maps spectral vectors of whisper to F 0 values of converted speech. In the conversion stage, we use a feedforward neural network to detect the voiced segments in whispered utterance and then compute F 0 for these segments only instead of computing these values for all segments.

Another innovative aspect of this paper is the language: we have applied voice conversion techniques to French, which has much more complex syllabic structures and a larger phonemic inventory than Japanese. Degraded performance is thus expected with respect to original published results.

The paper is organized as follows. Section 2 describes some characteristics of whispered speech. Section 3 describes the frameworks of our NAM-to-Speech conversion system. Section 4 describes our experimental evaluations and finally, conclusions are drawn in Section 5.

Whispered speech

In recent years, advances in wireless communication technology have led to the widespread use of mobile phones for private communication as well as information access using speech. Speaking loudly to a mobile phone in public places may be a nuisance to others, however, whispered speech can only be heard by a limited set of listeners surrounding the speaker and can therefore effectively be used for quiet and private communication over mobile phones [START_REF] Ito | Analysis and recognition of whispered speech[END_REF].

Acoustic features

In normal speech, voiced sounds involve a modulation of the air flow from the lungs by vibrations of the vocal folds. However, there is no vibration of the vocal folds in the production of whispered speech. Exhalation of air is used as the sound source, and the shape of the pharynx is adjusted such that the vocal folds do not vibrate. Due to this difference in production mechanism, the acoustic characteristics of whisper differ from those of normal speech. A study on the acoustic properties of vowel sounds [START_REF] Ito | Analysis and recognition of whispered speech[END_REF] has shown an upward shift of the formant frequencies for vowels in whispered speech compared to normal speech. The shift is larger for vowels with low formant frequencies. The authors also found that the cepstral distances between normal and whispered speech for vowels and voiced consonants are higher than those of unvoiced consonants. This means that the vocal tract characteristics of vowels and voiced consonants change more significantly in whisper relative to ordinary speech than those of unvoiced consonants.

The perception of vowel pitch in normal speech is related mainly to the fundamental frequency (F 0 ) which corresponds to periodic pulsing. In whispered speech, however, although there is no periodic pulsing, some pitch-like perception may occur. Higashikawa et al. [START_REF] Higashikawa | Perceived Pitch of Whispered Vowels -Relationship with formant frequencies: A preliminary study[END_REF] have shown that listeners can perceive pitch during whispering and formant frequency could be one of the cues used in perception. More precisely, the authors in [START_REF] Higashikawa | Acousticalperceptual correlates of "whispered pitch" in synthetically generated vowels[END_REF] indicated that "whisper pitch" is more influenced by simultaneous changes in F1 and F2 than by changes in only one of the formants.

NAM microphone

Nakajima et al. [START_REF] Nakajima | Non-audible murmur recognition[END_REF] proposed a new communication interface which can capture acoustic vibrations in the vocal tract from a sensor placed on the skin, below the ear. This position is shown in Fig. 1. This position allows a high quality recording of various types of body transmitted speech such as normal speech and whisper. Body tissue and lip radiation act as a lowpass filter and the high frequency components are attenuated. However, the non-audible murmur spectral components still provide sufficient information to distinguish and recognize sound accurately [START_REF] Heracleous | Audible (normal) speech and inaudible murmur recognition using NAM microphone[END_REF]. Currently, the NAM microphone can record sound with frequency components up to 4 kHz while being little sensitive to external noise.

NAM-to-Speech conversion system

Several approaches have been proposed to convert non audible speech to modal voice. The most trivial system consists in chaining NAM recognition with speech synthesis. Direct signal-to-signal mapping using aligned corpora is also very promising: Toda et al. [START_REF] Toda | NAM-to-Speech Conversion with Gaussian Mixture Models[END_REF] applied statistical feature mapping [START_REF] Stylianou | Continuous probabilistic transform for voice conversion[END_REF] [START_REF] Kain | Spectral voice conversion for text-to-speech synthesis[END_REF] to NAM-to-speech conversion.

Although the segmental intelligibility of synthetic signals computed by statistical feature mapping is quite acceptable, listeners have difficulty in chunking the speech continuum into meaningful words. A large part of this problem is due to impoverished synthetic intonation. In this study, we focus on improving the estimation of pitch and voicing of the converted speech. Fig. 2 shows the conversion used in our system. Before training the models for spectral estimation and F 0 estimation, the pairs of whisper and speech uttered by a speaker must be aligned because of different speaking rates. Transcription information was used for this task to get a better alignment compared to blind dynamic time warping (DTW).

Spectral Estimation

We use the same schema for spectral estimation as the one proposed by Toda [START_REF] Toda | NAM-to-Speech Conversion with Gaussian Mixture Models[END_REF]. As described in [START_REF] Toda | NAM-to-Speech Conversion with Gaussian Mixture Models[END_REF], feature vector X t of whisper consisting of spectral feature vectors of several frames around a current frame t was aligned with a target speech feature Y t = [y t , ∆(y t )] consisting of static and dynamic features. These vectors were then used to train a GMM for representing the joint probability density p(X t , Y t | Θ), where Θ denotes a set of GMM parameters. Another model for representing the probability density of global variance (GV) of the target static features p(v(y) | Θ v ) was trained where Θ v denotes a set of parameters of a Gaussian distribution, v(y) denotes global variance over the time sequence y of target static feature. This global variance information is used to alleviate the over-smoothing, which is inevitable in the conventional ML-based parameter estimation [START_REF] Toda | Voice Conversion Based on Maximum Likelihood Estimation of Spectral Parameter Trajectory[END_REF].

In the conversion, the target static feature y was estimated from source feature X = [X 1 , X 2 ,…, X T ] so that a likelihood L = p(Y|X, Θ) w p(v(y) |Θ v ) was maximize where w is a weight and the vector Y is represented as W y, where W denotes a conversion matrix from the static feature sequence to the static and dynamic feature sequence.

Excitation Estimation

The mixed excitation is defined as the frequency-dependent weighted sum of white noise and a pulse train with phase manipulation. The weight is determined based on an aperiodic component in each frequency band [START_REF] Kawahara | Aperiodicity extraction and control using mixed mode excitation and group delay manipulation for a high quality speech analysis, modification and synthesis system STRAIGHT[END_REF].

Aperiodic estimation was done in the same way as the spectral estimation except that global variance (GV) was not used because GV does not cause any large difference to the converted speech in the aperiodic conversion [START_REF] Ohtani | Maximum Likelihood Voice Conversion Based on GMM with STRAIGHT Mixed Excitation[END_REF].

ML-based conversion method was used the F 0 estimation. Static and dynamic features Y t of F 0 are used while keeping the same feature vector of whisper X t as that used for the spectral conversion. However, instead of using all the segments in each pair of utterance, only voiced segments X t , Y t were extracted to train a GMM on the joint probability in a similar way as the spectral estimation in order to avoid loosing some Gaussian components for representing the zero values of F 0 set for unvoiced segments. A feed-forward neural network is used to predict these segments from X. For synthesis, continuous F 0 values are predicted that are paced by the voicing parameter computed by the network.

Evaluation

In order to show our improvement in F 0 estimation and voicing attribution, two evaluations were done, comparing our system with the original system proposed in [START_REF] Toda | NAM-to-Speech Conversion with Gaussian Mixture Models[END_REF].

The training corpus consists in 200 utterance pairs of whisper and speech uttered by a French male speaker and captured by a NAM microphone and head-set microphone. Respective speech durations are 4.9 minutes for whisper (9.7 minutes with silences) and 4.8 minutes for speech (7.2 minutes with silences). The 0 th through 24 th mel-cepstral coefficients were used as a spectral feature at each frame. The spectral segment feature of whisper was constructed by concatenating feature vectors at current ± 8 frames, and then the vector dimension was reduced to 50 using a PCA technique. Log-scaled F 0 extracted by STRAIGHT [START_REF] Kawahara | Restructuring speech representations using a pitchadaptive time frequency smoothing and instantaneousfrequency-based F 0 extraction: Possible role of a repetitive structure in sounds[END_REF] was used as the target feature.

Voicing estimation

In the original system [START_REF] Toda | NAM-to-Speech Conversion with Gaussian Mixture Models[END_REF], first, the authors took all the frames in each utterance and estimated F 0 value at each frame by using the trained GMM model. Then, they used a threshold to assign the voiced/unvoiced label for this frame. The F 0 values in every unvoiced frames were then set to zero. We applied this original technique to our data.

Then to compare with our technique, we created a feedforward neural network with 50 input nodes, 17 hidden nodes and 1 output node. The segmental features at each frame of the whispered utterances were used as input vector for this network. The voiced/unvoiced label for each segment in the training whispered data was obtained from the voiced/ unvoiced label of the corresponding speech utterance by aligning the two utterances. All the whispered utterances used for training the GMM were also used to train this network. 1 shows the evaluation of this network. Compared with the error in the original system, the result shows that we have a slight improvement of the voiced/unvoiced detection.

Parameter evaluation

We also compared the two systems with different number of mixtures for estimating F 0 on both the training and the test data. The number of mixtures of mel-cepstral mapping function was set to 32. Full covariance matrices were used for both GMMs. The test corpus consisted of 70 utterance pairs not included in the training data. The error was calculated as the normalized difference between synthetic F 0 and natural F 0 in the voiced segments that were well detected by the two systems. The error is given by the following formula: Err = (synthetic_F 0 -natural_F 0 )/ natural_F 0 . Fig. 3 shows that the proposed framework outperforms the original system. In addition, when the number of Gaussian mixtures increase, the errors of both systems on the training data decrease, but these errors on test data are little sensitive to the number of mixtures. This is further illustrated in Table 2 which provide correlation coefficients for the two systems. Fig. 4 shows an example of whispered-, converted-(our system) and natural-speech. As can be seen the formant patterns of converted speech are flatter than those of natural speech. Global variance was used to attenuate this difference.

Perceptual evaluation

Sixteen French listeners who had never listened to NAM participated in our perceptual tests on intelligibility and naturalness of the converted speech from the two systems. We used 20 utterances which were not included in the training.

Evaluation on intelligibility

For this evaluation, we used synthetic speech with the estimated mel-cepstrum, estimated aperiodic component and estimated F 0 by both systems. Each listener heard an utterance pronounced in modal speech and the converted utterances obtained from the whispered speech with both systems. For each utterance, they were asked which one was closer to the original one, in terms of intelligibility (ABX test). Fig. 6a provides the mean intelligibility scores for all the listeners for the converted sentences using the original system and our new system. The intelligibility score is higher for the sentences obtained with our new system (F = 23.41, p <.001).

In the evaluation in [START_REF] Toda | NAM-to-Speech Conversion with Gaussian Mixture Models[END_REF], it was shown that F 0 estimation (compared with using a constant F 0 ) improves intelligibility. This might be caused by a better detection of the voiced/ unvoiced property in the F 0 estimation. In our case, a feedforward neural network was used instead of using GMM.

Evaluation on naturalness

We used here the version of the test utterances produced with the modal voice. These utterances are considered as ideal desired targets of the mapping. We thus further processed the synthetic utterances by warping their time scale to the targets using the warping procedure used for training.

We then conducted an ABX test. Because of the warping, all utterances have almost the same temporal organization. For each sentence, subjects choose A or B as the nearest to the natural X in terms of naturalness. Fig. 6b shows the mean naturalness that all the listeners rate as the nearest. Again the proposed system was strongly preferred to the original one (F = 74.89, p <.001).

Conclusions

This paper described the improvement in F 0 estimation and voicing decision we propose for NAM-to-speech conversion system applied to whispered speech. GMM models were used to estimate the spectra, aperiodic component and F 0 of the converted speech from spectral segments obtained from NAM-captured whispered speech, based on a maximum likelihood criterion. To estimate the F 0 features in the whispered utterances, only voiced segments were used. They were detected using a simple feed-forward neural network. Although the performance of the system is improved compared to that of the original system, the estimated pitch is still flat due to the GMMs. In the future, we will investigate how to obtain audible speech from whisper by using a HMM which is appropriate for modelling a time sequence of speech parameters. Also, we plan to use complementary information such as video in the aim of obtaining other useful parameters.

Figure 1 :

 1 Figure 1: Position of NAM microphone.

Figure 2 :

 2 Figure 2: Conversion process of NAM-to-Speech system.Spectral Estimation.

Figure 3 :

 3 Figure 3: Natural and synthetic F 0 curve.

Figure 5

 5 Figure 5 shows an example of a natural (target) F 0 curve and the synthetic F 0 curves generated by the two systems. It shows that our new system is closer to the natural F 0 curve than the original system.

Figure 5 :

 5 Figure 5: Natural and synthetic F 0 curve for the same utterance : "Armstrong tombe et s'envole".

Figure 6 :

 6 Intelligibility score (a) and Naturalness score (b).

Table 1 :

 1 Voicing error using neural network or GMM.

	Type of error	Feed-fwd NN (%)	GMM (%)
	Voiced error	2.4	3.3
	Unvoiced error	4.4	5.9
	Total	~ 6.8	~ 9.2

Table

Table 2 :

 2 Correlation coefficient between natural F 0 and converted F 0 by the two systems. Whispered speech captured by NAM sensor, converted speech and ordinary speech for the same utterance: "Armstrong tombe et s'envole".

	# of	Our system	Original system
	Gaussian mixtures	train	test	train	test
	8	0.565	0.495	0.494	0.451
	16	0.667	0.493	0.513	0.456
	32	0.746	0.498	0.603	0.460
	64	0.834	0.499	0.682	0.444
	/ Whisper captured by NAM microphone /
	/ Converted speech /
	/ Natural speech /
	Figure 4:

Acknowledgment: The authors are grateful to Coriandre Vilain and Alain Arnal for data acquisition, Prof. Hideki Kawahara of Wakayama University in Japan for the permission to use the STRAIGHT system.


Subramaniyam Kannan
email: kannan5@purdue.edu

Rebecca S Salo
email: rsalo@purdue.edu

Paul C Wood
email: pwood@purdue.edu

Saurabh Bagchi
email: sbagchi@purdue.edu

Somali Chaterji
email: schaterj@purdue.edu

OPTIMISM: Optimal Placement of Intrusion Detection Sensors to Identify Multi-Stage Attacks in Software Defined Networks

A major threat to network security is multi-stage attacks, where an attacker compromises an external-facing server and then penetrates into the servers deeper in the network one by one by using the elevated access privilege at each compromised node. She advances through these steps till she gains access to protected information present in the crown jewel of the network. Intrusion detection sensors (IDS) can detect such attacks, but limited available resources constrain the number of IDS that can be deployed. Software defined networking (SDN) provides network flexibility, and combined with network function virtualization (NFV), it enables agile IDS placement optimizations. We present OPTIMISM, a novel system for placing IDS to maximize network protection or benefit and minimize total deployment costs. The overall aim of OPTIMISM is to select an IDS placement configuration such that the probability of attacker reaching the crown jewel is reduced and the corresponding deployment cost is within the network administrator's budget. A Genetic Algorithm (GA) is used for optimizing the IDS placement that yields the maximum benefit subject to a cost constraint. OPTIMISM provides an option to determine the placement configuration that has the maximum return on investment and also dynamically reconfigures the placement of IDS using NFV depending on runtime alerts. Using a simulation, we show that OPTIMISM selects IDS placement configurations that have either maximum benefit given a budget or maximum return on investment. Using our system, network administrators can maximize network protection and reduce the overall costs of protecting their systems from multi-stage attacks.

I. INTRODUCTION

Multi-stage attacks (MSA) [START_REF] Clark | The problem isn't attribution: it's multi-stage attacks[END_REF] are attacks where the attacker compromises an outward facing server and uses that to gain access to subsequent servers, many of which are not directly accessible from outside the network, until reaching the target server or crown jewel deep inside the network. These threats are especially difficult to defend against because they require placing network security appliances, such as intrusion detection sensors (IDS) [START_REF] Liao | Intrusion detection system: A comprehensive review[END_REF], at points throughout the network instead of just at the edges of a system. Only network administrators with unlimited budgets can afford to place IDS at every location in a network, leaving the majority of defenders with a challenging placement problem. In this paper, we provide a budget-conscious approach for placing IDS in a network to optimize their effectiveness at stopping multi-stage attacks while meeting budget constraints.

The naïve approach for placing IDS to protect against MSA is to place an IDS in every server, providing the maximum probability of attacker detection. However, this is impractical in large networks as each IDS would consume substantial amounts of computing and network resources and in turn significantly i ncrease t otal p rotection c osts. Administrators need to optimize the placement of IDS in the face of MSA, and since each IDS placement point has a different cost and impact on the quality of the defense, there is ample opportunity to improve efficiency. F or e xample, a n I DS t hat i s n ot on a path to the crown jewel provides no value, while an IDS on a highly utilized network link may require very expensive processing hardware to avoid adding latency. Through models and simulations, the number of IDS deployed and on which servers they are deployed can be carefully controlled to minimize deployment costs and maximize attacker detection [START_REF] Abu Sharkh | Resource allocation in a network-based cloud computing environment: design challenges[END_REF], [START_REF] Modelo-Howard | Secure configuration of intrusion detection sensors for dynamic enterprise-class distributed systems[END_REF].

Several existing methods have tried to solve the problem of optimizing IDS placement within a traditional network. Noel et. al [START_REF] Noel | Optimal ids sensor placement and alert prioritization using attack graphs[END_REF] and Modelo-Howard [START_REF] Modelo-Howard | Determining placement of intrusion detectors for a distributed application through bayesian network modeling[END_REF][START_REF] Qian | Information assurance: dependability and security in networked systems[END_REF] use attack graphs to predict network vulnerabilities and perform alert inferencing to decide on IDS placements. In practice, however, it is difficult t o g enerate a ttack g raphs a nd t o k eep them updated as the networks and applications are dynamically changing with new vulnerabilities and attack opportunities. Our system OPTIMISM (OPTIMal Intrusion SysteM) operates without relying on attack graphs and instead continuously monitors the network traffic to dynamically decide on the IDS placements in response to both alerts and changing network topologies and traffic. A nother s ystem d esigned t o optimize the location of IDS to stop Denial of Service (DoS) attacks is BroFlow [START_REF] Lopez | An elastic intrusion detection system for software networks[END_REF], which places IDS to maximize network coverage but only for DoS. Hasan and Mouftah [START_REF] Md | Optimal trust system placement in smart grid scada networks[END_REF] proposed an IDS placement system that uses network segmenting within the larger network and places IDS at the segment border nodes. If a multi-stage attack penetrates the border nodes, however, then there is no way to detect the attacker inside these segments. Our solution is focused on protecting specific a ssets i n the 1 system rather than general network segments, allowing the network administrator to focus resources on targets of interest rather than entire network segments.

Software defined networks (SDN) [START_REF] Kim | Improving network management with software defined networking[END_REF] are becoming increasingly common, and OPTIMISM specifically leverages SDN's ability to optimize the IDS placement problem for MSA. The centralized SDN controller oversees the network, and it allows flexibility in routing thus simplifying network reconfiguration [START_REF] Kreutz | Software-defined networking: A comprehensive survey[END_REF]. Similarly, through network function virtualization (NFV) [START_REF] Han | Network function virtualization: Challenges and opportunities for innovations[END_REF] techniques, network functions such as IDS can be virtualized, liberating them from being pinned to specific hardware and thus enabling flexible deployment in response to network changes. Using a combination of NFV and SDN, IDS can be dynamically deployed, either by moving an IDS to a new virtual machine or by re-routing traffic to the desired virtual machine using the SDN controller. The solution presented in this paper is the first step in developing a dynamic IDS deployment system where placements can evolve over time.

In this work, we develop an algorithm that minimizes the cost of deploying IDS in SDN while maximizing their effectiveness at mitigating multi-stage attacks by scoring and balancing the benefits and deployment costs for a network. The benefit and deployment cost are calculated for each configuration by analyzing all the potential attack paths in the network, and the configuration which produces the highest benefit for a cost constraint is sought in an iterative manner either by using the genetic algorithm or through exhaustive search. Additionally, OPTIMISM provides a best value optimization to find the configuration for which the benefit-tocost ratio is maximized. In order to stop new attacks through newly discovered attack paths, OPTIMISM runs a dynamic version of the algorithm in response to alerts or updated attacker information to dynamically move the IDS, using the functionality of NFV, to maximize protection.

OPTIMISM utilizes a directed acyclic graph to represent all of the server instances (nodes) and network connections (edges) in a protected network. A set of edge nodes connects the attacker and the legitimate traffic to the modeled network, and each edge gets two kinds of weights. The first weight represents the probability that an attacker is able to travel over the link (after penetrating the source server), and the second weight represents the amount of legitimate network traffic that traverses on that edge. We define the benefit as the probability that an attacker reaches a protected asset, a known protected node, and the deployment cost as a numerical function of installation and processing costs. From this information, each IDS placement, i.e., a mapping of a sensor being placed (or not) on a specific node, is scored in both benefit and cost. The optimizer in OPTIMISM is responsible for selecting a placement that maximizes benefit while either meeting a cost constraint or maximizing the return on investment through a multi-objective function. Once a placement is known, it can be dynamically improved by updating the information about the traffic or attacker paths in the system. We model alerts as increasing the probability an attacker is on a particular path and dynamically update IDS placements so that the attacker can be found more quickly.

First we compare our runtime-optimized results to an exhaustive search technique. For networks with more than 5 servers, our system consistently maintains over 85% of the maximum benefit, which corresponds to placement of an IDS on every node, while keeping costs as low as 30% of the naïve complete-coverage cost. The error in the performance of the configuration found by the genetic algorithm is less than 5% when compared to exhaustive search techniques; and for networks larger than 10 servers, the genetic algorithm running time is less than a third of the exhaustive search running time. Our dynamic adaptation further improves the static placement by an average of 38.2% by incorporating alert information into the system model and re-distributing the IDS placements.

In this paper, we make the following novel contributions: 1) OPTIMISM optimizes IDS placement for multi-stage attacks without relying on attack graphs by giving higher priority to placements deeper in the network, both in large and small networks. 2) OPTIMISM accounts for network structure and traffic to accommodate distributed applications. Thus, its IDS placement strategy incorporates both IDS quality and traffic processing costs. 3) OPTIMISM solves two rigorous optimization problems, which provide a choice between maximum performance given a budget (cost constraint optimization) or maximum return on investment (best value optimization). 4) OPTIMISM generates nearly optimal solutions for IDS deployment using genetic algorithms with faster run time than exhaustive search. It also provides dynamic IDS deployments to stop emergent attacks.

II. SYSTEM DESIGN

OPTIMISM is a system that provides optimal IDS placement decisions over time. For example, if there are 10 servers where an IDS can be placed, then OPTIMISM will provide a vector of 10 entries indicating whether or not an IDS should be located at each server. In order to select the best placement, a cost and benefit model for the protected network must be developed, and it is described in this section.

A. Protected System Model

OPTIMISM is designed to protect software defined enterprise networks against MSA. Such attacks are executed on distributed systems where multiple independent server instances separate the attacker from the target server or crown jewel. Fig. 1 shows such a system where a potential attacker must compromise a web, database, and reporting server prior to having direct access to compromise the file server with the crown jewel. We assume that all the servers are connected by an open flow switch [START_REF] Mckeown | Openflow: enabling innovation in campus networks[END_REF] and controlled by a centralized SDN controller [START_REF] Berde | Onos: towards an open, distributed sdn os[END_REF] as shown in Fig. 1 as this allows both knowledge of the network topology and control of routes through IDS placements. We further assume each server is a virtual machine which supports NFV and thus allows local IDS Each IDS has a high computational requirement, and since the network uses SDN and NFV, a the set of high-powered IDS-capable physical machines can be re-purposed to other layers in the distributed system and re-connected to the architecture with the flexible SDN routes. We assume that the network administrator has a finite budget for these IDS-capable machines.

B. Attacker Model

The attacker we seek to detect is a multi-stage attacker who exhibits abnormal behavior such as brute force attacks or known-exploits attempts, and such activities trigger alerts if they occur on an IDS-protected server. In Fig. 1, the attacker first compromises any one of the web servers at stage 1 and then proceeds further to rest of the stages by leveraging elevated access privileges that she gained at the previous stages. OPTIMISM assumes that the attacker repeats the attack through the same attack path again until exhausting all attack opportunity as she has the advantage of some already compromised servers in previous tries. Consequently, the attacks are repeated along the same path until a certain number of trials are reached and then a new path is selected.

The attacker model is initialized as having a uniformly random chance to take each path from the edge to the crown jewel. Expert knowledge can refine this assumption or it can be refined based upon prior alerts, as we do in OPTIMISM.

C. IDS Quality Value

IDS quality values represent the probability that the IDS at a particular server generates a strong alert for an attacker, stopping her particular attack attempt. We model the IDS quality as higher if it is located further inside the network since it will signal fewer but stronger alerts. The traffic is less erratic and cleaner deeper in the system since it interacts between known, controlled applications, and thus it is easier to detect abnormal activity. The IDS quality values can be provided by the network administrator based on prior observations of the IDS performance at each server in the system [START_REF] Pietraszek | Using adaptive alert classification to reduce false positives in intrusion detection[END_REF] or based on the self-reported quality and frequency of alerts in off-theshelf systems such as Snort [15].

D. Network Traffic Distribution

OPTIMISM penalizes the IDS placement based on the amount of network traffic that must be processed by each IDS, thus dictating cost. In the context of SDN, open flow switches [16] have default counters to keep track of the number of packets received at each port. OPTIMISM can generate a traffic distribution model by observing the counters for certain period of time. Alternatively, the model can be updated by the network administrator who has application-layer knowledge of the inter-server traffic models. In OPTIMISM, the traffic distribution is modeled using weights that represent the amount of legitimate network traffic flowing through the particular edge or network link. The amount of traffic is then calculated against a penalty function with a ceiling to capture the computational limits of a single IDS.

E. Cost and Benefit

OPTIMISM maximizes benefit while minimizing costs to provide optimal IDS placements. The cost comes directly from the placement count and the amount of traffic at each placement, including the cost to install and operate the IDS. The benefit comes from the reduction in attacker's success or the probability that the attacker is able to penetrate the crown jewel. Since costs can be elastic, we provide two models for optimization: cost constrained (CCO) and best value (BVO) where CCO maximizes the benefit subject to a binary constraint that the costs do not exceed a budget and BVO maximizes the return on investment or the benefit-to-cost ratio.

F. Dynamic Updates

Each IDS placement in OPTIMISM is selected based upon the current network conditions and attacker model. As time progresses, the network model and attacker model evolve so that the optimal placement changes. In OPTIMISM, we provide a method to both update the attacker model by incorporating alerts into path or edge weighting and a method to quickly incorporate such information into new optimal placements that can be leveraged by the NFV/SDN system to improve detection.

III. METHODS

The goal of OPTIMISM is to reduce the costs of protecting a network against MSA by limiting the number of IDS deployed while maximizing the protection provided by them. The algorithm we have developed takes as inputs (i) the network graph, (ii) the total incoming network traffic, (iii) information about the security value of IDS, (iv) cost information, including IDS installation costs and running costs, and (v) the objective function selection (CCO vs BVO). The output of OPTIMISM is a particular placement of IDS in the network, and this section describes how the placement is found.

A. The Network Graph and Problem Statement

OPTIMISM models the network by using a directed acyclic graph that represents the potential flows (directed edges) between servers (nodes) in the protected system. The graph has two sets of special nodes: the edge servers (ES) and the crown jewel (CJ) with |CJ| = 1 for our simplifying assumption. The attacker progresses from one of the edge servers toward the crown jewel unimpeded unless an IDS is present on a particular server, which we define as a binary vector I P of length N -|CJ| where I P (i) = 1 indicates IDS presence at server i and N is the number of protect-able servers in the system, excluding the crown jewel.

We define two abstract functions that will be fully defined in the later sections to formulate the problem: f B (I P , I Q , G, EW A , CJ, ES) is the benefit or protection quality provided by a placement I P for IDS's of a particular quality (I Q ) on a particular network graph (G) with knowledge of the attackers paths (EW A ), and f C (I P , I C , G, EW C , ES) is the numerical cost for the placement for cost factors I C and traffic information in EW C . From these functions, we define two problems to optimize. Cost constrained:

arg max I P f B (I P , I Q , G, EW A , CJ, ES) (1)
subject to

f C (I P , I C , G, EW C , ES) < B ( 2
)
where B is the cost constraining budget. This definition is called the Cost Constrained Optimization (CCO) problem. The secondary form of this problem combines these objectives into a single equation:

arg max I P (1 -α) • f B (I P , ...) -α • f C (I P , ...) (3)
where α (weight) controls the relative importance of each function such that high alpha values favor cost over protection and low values favor protection over budget reductions. This definition is called the Best Value Optimization (BVO) because it finds the maximum benefit-to-cost ratio.

B. Benefit Calculation

We wish to calculate the benefit of a particular placement f B (I P , ..) of IDS inside a network. First, we define the benefit in terms of the probability that an attacker is successful:

benefit = 1 -P (AS) (4)
where the attacker being successful (P (AS = 1)) means that she has reached the crown jewel.

For each node in the graph, we define two probabilities: P I (i), the probability an attacker reaches the input of node i and P O (i), the probability an attacker leaves the node. These parameters are related by:

P O (i) = I P (i) • (1 -I Q (i)) • P I (i) + (1 -I P (i)) • P I (i) (5)
where I Q (i) is the quality of the IDS (∈ [0, 1]), I P (i) indicates placement, P I (i) is the input probability, and P O (i) is the output probability. If an IDS is placed, then the probability is Output attacker probability P T (i, j) Attacker's probability to take edge i, j P (AS)

Probability of attacker's success reduced by the IDS quality, and if it is not, then the probability passes through the node (P O (i) = P I (i)). Each node has a set of predecessors Preds(i) and successors Succ(i) (the input and output arcs), and the input probabilities are calculated as follows:

P I (i) = j ∈ Preds(i) P T (i, j) (6)
where P T (i, j) is the attacker's probability to follow a path from node i to j:

P T (i, j) = P O (j) • EW A (j, i) k ∈ Succ(j) EW A (j, k) (7)
For nodes with |Preds(i)| = 0, they have P I (i) = 0:

P I (i) = 0 ∀ i / ∈ (ES ∪ {i : |Preds(i)| ≥ 0}) (8)
For nodes in the edge server set ES, they are initialized with:

P I (i) = 1 |ES| ∀ i ∈ ES (9)
For a particular placement I P the probabilities P I /P O are calculated algorithmically:

1) Initialize all P I (i) = -1 ∀ i 2) Apply Eq. ( 8) , ( 9) 3) Repeat Eq. ( 6) , (5), skipping any nodes i that have (∃ j s.t. P O (j) == -1 ∧ j ∈ Preds(i)), until P I (i) ≥ 0 ∀ i From this calculation, we can directly calculate the attacker's success:

P (AS) = P I (CJ) (10)
thus f B = 1 -P I (CJ), the benefit for placement I P .

C. Cost Calculation

The cost for each placement is the sum of the installation costs and running costs for servers which have IDS (I P (i) is 1). For notational simplicity, I C contains multiple cost parameters:

• Installation Cost (I I C ): The cost that is incurred to install an IDS in a virtual machine. It includes the IDS licensing cost.

• Running Cost (I R C ): The amount of processing power spent by the VM where an IDS is deployed in order to analysis the traffic flowing in. The processing cost varies depending on the location where an ids is deployed and the amount of traffic it has to process. The parameter is defined in terms of KP , a parameter that allows a trafficto-CPU core mapping. Each IDS process a certain number of packets per second (T (i)), and each IDS placed node has the following cost:

C T (i) = I P (i) • ( I R C (i)•T (i) KP ) + I I C (i) T (i) KP < I M C (i) I P (i) • inf T (i) KP > I M C (i) (11
) where I M C (i) is the cost representing the maximum traffic that an IDS at node i can process and T (i) is the traffic into a particular node which we define as:

T (i) = P I (i) • T E for I P (i) = 0 ∀ i (12)
where T E is the total input traffic to the edge nodes ES in the system, and I P (i) = 0 ∀ i is the probability calculation for a system with no IDS placed. Whenever this calculation is run, the term EW A is replaced by EW C in Eq. [START_REF] Lopez | An elastic intrusion detection system for software networks[END_REF].

Then the final function f C :

f C = i∈N C T (i) (13)
where the parameters I P , I C , G, EW C , ES are passed to their respective equations from the prior section.

D. Placement Optimization

Two optimization problems were presented in Eq. ( 1) and (3) for placing IDS, e.g., optimizing I P . Practically, I P grows at 2 N placement options, so exhaustively finding the best solution is impractical, especially at runtime. To solve these optimization problems, the genetic algorithm (GA) is used, and it is validated for small networks with an exhaustive search method in our experimental results.

1) Genetic Algorithm Optimization: The GA was chosen as the optimization method because of its ability to quickly determine a nearly optimal solution in our large binary search space. Prior work [START_REF] Hao Chen | A multi-objective optimisation approach to ids sensor placement[END_REF] has also used the GA to determine IDS placement, though not in a multi-stage attack scenario. Given the large number of configuration options for even a mediumsized network, it would take far too long to perform an exhaustive search. Likewise, a linear programming approach would not work because the objective function is not a linear function of the IDS placement.

The GA is based on the process of natural selection, and it starts with an arbitrarily chosen IDS placement vector and evaluates the objective function [START_REF] Goldberg | Genetic Algorithms in Search, Optimization and Machine Learning[END_REF], which is the benefit equation (given by eq. 1 or eq. 3) for the cost-constrained problem. After each generation, the algorithm selects parents which have the best outcomes of the objective function. Crossing and mutating the parent vectors leads to the generation of new offspring vectors [START_REF] Deep | A real coded genetic algorithm for solving integer and mixed integer optimization problems[END_REF], for which the objective function is again calculated. This leads to a nearly optimal solution. It is not necessarily the best possible solution because the GA only evaluates the objective function for a limited number of placement vectors until it cannot make any improvement by performing crossovers or mutations, but it is very close to the best possible solution.

In this constrained problem, each iteration of the objective function calculation also includes the calculation of the cost. If the cost exceeds the budget constraint (B), the solution is discarded and will not be used in subsequent crossovers or mutations. However, if the solution is valid, it will be taken as a possible solution and used in crossovers and mutations. The specific parameters of the GA are described in Section IV-A1.

2) Exhaustive Search Optimization: To validate the placement given by the CCO-GA, we also use an cost constraint exhaustive search (CCO-ES) for small networks. The exhaustive search method calculates the benefit for all possible IDS placements, discards those which have a cost that exceeds the budget, and then chooses the placement with the maximum budget.

While the exhaustive search does yield the optimal solution, it is only practical for small networks with fewer than 15 servers due to the large number of placement options. Thus, it can be used to validate the GA optimization for small networks, but the GA is necessary for finding the nearly optimal placement for large networks.

E. Dynamic Adaption

The model parameter EW A is expected to be time-varying as attackers enter and leave the system. For a static analysis without any external information, EW A (i, j) = 1 ∀ i, j in the edge set of G, thus a uniform distribution is utilized for attacker probability calculations. This is a good defense setup when the attacker's path is not known in advance or if the attacker has no information about the system and simply makes random choices when attacking. Practically, however, the attacker will stick to a particular path until it cannot be exploited anymore (depth first search) because of the sunk costs of exploiting servers in the path to the crown jewel. Since the attacker maintains a constant path, we can utilize alerts generated by the IDS during each attack trial to adjust the EW A parameters:

EW A (i, j) = EW A (i, j) + 1 ∀ (i, j) ∈ AP ( 14
)
where AP is the set of all edges in the potential alert path. This path is calculated by backtracking from the predecessors of the node in the graph that generated an alert. One additional modification is made to the original problem. The distribution in Eq. ( 9) is modified by adding an additional node n a to the graph to represent the attacker. This node is added to represent the attacker's choice of which edge server to attack, thus modifying the equation to:

P I (i) = EW A (n a , i) j ∈ ES EW A (n a , j) ∀ i ∈ ES (15)
where n a is the attacker's virtual node in the graph. With this modification, the EW A values will grow along any candidate path to the alert that is being generated and bias the placements toward that path. Whenever the attacker is found, this weighting is reset back to the static case.

Runtime Adaptation: At runtime, the optimal placements can be varied by updating EW A and re-calculating the benefit f B . Whenever the optimal placement vector I P changes, then the NFV/SDN component can move resources to implement the new placement. As alerts come in from the IDS, the EW A vector is updated and the placements are re-calculated.

IV. EXPERIMENTATION

In this section, first we discuss the simulation setup used for OPTIMISM and then discuss the experiments that evaluate the performance of the OPTIMISM.

A. Experimental Setup

For the experiments in this section, unless mentioned, the parameter values used are in Table II and the graph in Fig. 2 is used for the network. For the I Q term, we approximate depth-based improvements for the servers as follows:

I Q (i) = IDS B + IDS C • max (dist(j, i) ∀ j ∈ ES) (16)
where dist(j, i) is the distance from node j to node i in the graph with each edge counted as 1. This gives slight benefits to IDS placed deeper in the system where both the source and destination are controlled by the network administrator.

1) Genetic Algorithm Settings: We utilize an integer version of a genetic algorithm1 with a population size of 200 that is initially chosen at random to form I P . For the selection function, a stochastic remainder is used that makes discrete choices using a roulette function. The elite count is 5% of the population, and the adaptive feasible mutation function is used for the CCO problems that stochastically mutates within the binary bound of the placement problem. The crossover function uses intermediate values that are selected stochastically but weighted by the performance of the population. The GA is run with 8 stall generations and it terminates whenever the average variance in the objective function (BVO or CCO) is less than 1e-4 between stall generations.

B. Evaluation Metrics

In addition to the benefit (f B ) defined in Eq. ( 4), two additional metrics are utilized:

Performance:

f B (I P ) f B (I P M ) (17)
where

I P M (i) = 1 ∀ i ∈ G.
For this metric, a value of 1 means that it is performing the same as the full-coverage case. % Max Cost or Relative Cost:

f C (I P ) f C (I P M ) (18)
with the same I P M from above. This captures a percent reduction in a relative scale to measure the improved efficiency of OPTIMISM.

C. Experiment 1: Performance of OPTIMISM vs naive approach

This experiment demonstrates the improvement in performance of OPTIMISM's CCO-GA over naively placing IDS. The naive placement approach is simulated by considering all the possible placements of IDS and selecting a mean and worst case performance. At each relative cost or user's budget (B from Eq. ( 2)), the minimum and mean performance for all possible combinations is calculated and then compared with the performance of the CCO-GA optimizer. Since the budgets do not align directly, they are clustered into 20 bins on the x-axis.

As shown in the figure 3, the CCO-GA consistently has a higher performance than both the mean and minimum performance of naive placement. The CCO-GA gives the optimal solution with a maximum performance at a deployment cost of 50% of the maximum cost. In other words, OPTIMISM deploys IDS only at 3 servers (in particular at node 7,9 and 10) out of 11 servers to get maximum path coverage to the crown jewel. For relative costs of 20-50%, the GA outperforms the naive placement by 5-20%, and it can outperform the minimum performance of the naive placement by up to 38%. Thus, OPTIMISM optimizes IDS placements by giving higher priority to nodes that are deeper in the network.

D. Experiment 2: Running Time of OPTIMISM's Genetic Algorithm vs. Exhaustive Search

This experiment validates our choice of a GA for solving our optimization problems. For the CCO, a series of 10 tests were run for each network size of a randomly generated network, starting at 5 servers and increasing by 2 until reaching 15 servers. The average runtimes were calculated for each network size and as shown in Fig. 4, the exhaustive search runtime varies exponentially with respect to the number of servers, while the GA runtime is linear. This shows that for large networks, it is much more practical to use the GA.

This experiment also measured the error of the GA with respect to the exhaustive search, calculated by the difference in the performance of the final placements. For 82% of the tests, where both network size and budget were varied, no difference was found between the two solutions; the error was 0%. Thus the median error was 0%, but the standard deviation was 6% due to a few trials with a maximum error of 49%. Despite this variation, the GA was concluded to be sufficiently accurate for placing IDS in a network of a large size. Thus OPTIMISM determines nearly optimal IDS placements by accounting for the network structure all in linear run time. This experiment evaluates the behavior of CCO vs BVO. Fig. 5 shows the variation in performance and relative deployment cost for a defined budget or cost constraint. As expected, the relative cost of the actual deployment has an almost linear relationship with the budget; any deviation from this is due to the cost being a discrete function of the number of servers.

The performance has a rapid increase at first, reflecting the added benefit of the initial IDS. However, at a budget of approximately 25%, the slope suddenly changes. This is due to the fact that placing more IDS means that there is a smaller chance for a subsequent IDS to detect the attacker. At some point, the benefit of adding another IDS approaches 0 because the existing IDS sufficiently cover the network.

The points on the graph which correspond to the various αvalues indicate the performance of the best value optimization. The points chosen by the BVO are the points in their local areas at which there is the greatest distance between the performance and relative cost. This is a direct result of the bestvalue behavior. Additionally, when α is lower, less emphasis is placed on the cost, and therefore the performance is better. When more emphasis is placed on the cost, the performance is lower. In this case, the performance has such a steep slope that a slight variation in α to emphasize cost causes no IDS to be placed.

The constrained optimization clearly shows a decreasing marginal performance as the cost increases. To a network administrator, this corresponds to having a lower-value security system. The best value optimization allows the network administrator to find the placement that yields the maximum return on investment by adjusting the weight (α) placed on the cost.

F. Experiment 4: Performance of Dynamic IDS Placements

This experiment evaluates the effectiveness of dynamically updating attacker edge weights (EW A ) on alert generation to determine new placements. We evaluate this approach by selecting 10 static paths at random from the edge nodes to the crown jewel. These paths are a set of edges which the attacker always tries to attack, and this is repeated 50 times. For each trial, a stochastic sample is taken for each node along the path as rand() > P O (i) where rand is the uniform random distribution between 0 and 1. If the value of the sample is 0, then the attacker is stopped, an alert generated (thus updating EW A ) and the trial is repeated. If the value is 1 then the attacker proceeds to the next server in the system. If the node i == CJ then the attacker is successful.

In this scenario, we use the BVO-GA approach because it allows a certain amount of dynamism with the cost in the system. Since α remains fixed, as the EW A increases, the benefit of particular placements becomes more significant than their costs. This allows for more than simple relocation as it enables additional IDS to be brought online in an elastic manner.

For this experiment, we measure the benefit via sampling so that in a set of 50 trials, if 25 reach the crown jewel, the benefit is 0.5. Fig. 6 shows the performance improvement of the dynamic system over the static allocation of IDS. The average improvement is 38.2%, significantly boosting the effectiveness of the defense.

V. CONCLUSION

OPTIMISM provides the network administrators with a costeffective method for protecting against multi-stage attacks. By considering the traffic distribution through each server, network structure and IDS quality, OPTIMISM optimizes the benefit function for each configuration iteratively to select the optimal IDS placement. Our method improves upon naive placement by up to 38% and reduces costs up to 70% when compared to the naive full-coverage solution. Due to the speed at which the genetic algorithm is able to solve the optimization problem, a new configuration can be generated and deployed whenever significant changes in the network occur. This is particularly useful in SDN, and by deploying the IDS using NFV, the placement can be quickly adjusted. This method of placing IDS is a novel approach to network defense, and it provides a strategy to effectively protect a network from multistage attacks in cases where information about the attacker is limited and costs must be minimized. We further provide a dynamic adaptation module that incorporates alert information to better defend against attacks by relocating IDS using NFV.

Fig. 1 .

 1 Fig.1. An example distributed enterprise system that OPTIMISM optimizes.

Fig. 2 .

 2 Fig. 2. Network Graph for the sample distributed system given in Fig. 1, with the edge values P T shown.

Fig. 3 .

 3 Fig. 3. The performance of OPTIMISM's cost constrained optimization with genetic algorithm (CCO-GA) is shown with respect to a naive IDS placement approach.

Fig. 4 .

 4 Fig.[START_REF] Modelo-Howard | Secure configuration of intrusion detection sensors for dynamic enterprise-class distributed systems[END_REF]. The exhaustive search runtime increases exponentially with respect to the number of servers, while the genetic algorithm runtime increases linearly.

8 Fig. 5 .

 85 Fig.5. Multi-objective optimization choses the points on the cost-constrained curve for which there is the most difference between performance and relative cost.

Fig. 6 .

 6 Fig. 6. The benefit (1 -P (AS)) is shown for 10 different random paths in the test network for 50 attack trials utilizing BVO-GA in a static and dynamic deployment. The dynamic case has a 38.2% improvement on average.

MATLAB's "ga" function


Vinh-Thong Ta

Abderrahim Elmoataz

Olivier Lézoray

PARTIAL DIFFERENCE EQUATIONS ON GRAPHS FOR MATHEMATICAL MORPHOLOGY OPERATORS OVER IMAGES AND MANIFOLDS

Keywords: Mathematical Morphology, PDEs, Partial difference, Graphs, Non local

The main tools of Mathematical Morphology are a broad class of nonlinear image operators. They can be defined in terms of algebraic set operators or as Partial Differential Equations (PDEs). We propose a framework of partial difference equations on arbitrary graphs for introducing and analyzing morphological operators in local and non local configurations. The proposed framework unifies the classical local PDEsbased morphology for image processing, generalizes them for non local configurations and extends them to the processing of any discrete data living on graphs.

INTRODUCTION

Nonlinear scale-space approaches based on Mathematical Morphology operators are one of the most important tools in image processing. Dilation and erosion are the two fundamental operators. They form the basis of many other morphological processes such as opening, closing, reconstruction, levelings, etc [START_REF] Soille | Morphological Image Analysis -Principles and Applications[END_REF]. Let f 0 : R 2 →R be a function representing a 2-D signal and let a structuring function g: B→R 2 representing some structuring element with a compact support B ⊆ R 2 . Dilation δ and erosion ε of f 0 by g are defined as δ g (f 0 ) = max f 0 (x-x ′ , y-y ′ )+g(x ′ , y ′ ) and ε g (f 0 ) = min f 0 (x+x ′ , y+y ′ )-g(x ′ , y ′ ) with x, y ∈ f 0 and x ′ , y ′ ∈ B. In this paper we consider the flat morphology case that use structuring set, i.e. for all x ′ , y ′ ∈ B the function g(x ′ , y ′ )=0. Such dilation and erosion operators are frequently implemented by algebraic set operations. For convex structuring elements, an alternative formulation in terms of Partial Differential Equations (PDEs) has also been proposed [START_REF] Maragos | Partial differental equations for morphological scale-spaces and eikonal applications[END_REF]. Given a disc B= z ∈ R 2 : |z|<1 one considers the following evolution equation if one assumes that at time t=0 the evolution is initialized with f (x, y, 0)=f 0 (x, y), solution of f (x, y, t) at time t>0 provides dilation (with the plus sign) or erosion (with the minus sign) within a disc of radius t. Such PDEs produce continuous scale morphology and offer several advantages [START_REF] Maragos | Partial differental equations for morphological scale-spaces and eikonal applications[END_REF]. First, they offer excellent results for non-digitally scalable structuring elements whose shapes cannot be represented correctly on a discrete grid (e.g. discs or ellipses). Second, they allow sub-pixel accuracy and can be adaptive by introducing a local speed evolution term [START_REF] Breuß | Anisotropic continuous-scale morphology[END_REF]. However, they have several drawbacks. First, the numerical discretization requires a careful choice of the time and spatial size. Second, this discretization is difficult for high-dimensional data or irregular domains. Finally, they only consider local interactions on the data by using local derivatives while non local schemes have recently received a lot of attention [START_REF] Gilboa | Nonlocal operators with applications to image processing[END_REF]5].

We propose to extend such PDEs-based operators to a discrete and non local scheme by considering partial difference equations over weighted graphs of the arbitrary topologies. To this aim, we explicitly introduce discrete non local derivatives on graphs to transcribe such continuous operators to operators based on graphs. Our proposal has the following advantages. No discretization step is needed and any data lying on any domain (even irregular) can be processed. The same formalism directly integrates local to non local information from data. This paper is organized as follows. Section 2 introduces partial difference equations over weighted graphs. In Section 3, we propose a family of non local weighted dilation and erosion (the p-dilation and the p-erosion). Section 4 provides several experiments for the case of image and manifold processing. Finally, Section 5 concludes.

PARTIAL DIFFERENCE EQUATIONS ON WEIGHTED GRAPHS

Definitions and Notations

We consider any general discrete domain as a weighted graph G=(V, E, w) composed of a finite set of vertices V = u 1 , . . . , u n and a finite subset of edges E ⊆ V ×V . An edge of E which connects the vertices u to v is denoted uv. G is weighted if it is associated with a weight function w: E→R + satisfying w uv =0 if uv / ∈ E, and w uv =w vu otherwise. We assume that G is connected, simple, undirected and weighted. We consider that a function f : V →R is provided that assigns a real value f (u) to each vertex u ∈ V .

A Family of Discrete Gradient

We consider the directional derivative of f at vertex u along the edge uv ∈ E to be defined as in [5]:

∂f ∂uv u = ∂ v f (u) = w 1/2 uv f (v)-f (u) . It satisfies ∂ v f (u)=-∂ u f (v), ∂ v f (v)=0, and if f (u)=f (v) then ∂ v f (u)=0.
Using these definitions, we can define two other directional derivatives:

∂ + v f (u)= max 0, ∂ v f (u) = w 1/2 uv max 0, f (v)-f (u) ∂ - v f (u)= min 0, ∂ v f (u) = w 1/2 uv min 0, f (v)-f (u) (1)
The weighted gradient operator of f , at vertex u ∈ V , is the following column vector. For all edges (uv i ) ∈ E:

∇ w f (u)= ∂ v f (u): v∼u T = ∂ v1 f (u), . . . , ∂ vi f (u) T ∇ ± w f (u)= ∂ ± v f (u): v∼u T = ∂ ± v1 f (u), . . . , ∂ ± vi f (u) T (2)
where notation v∼u means that the vertex v is adjacent to the vertex u, and ∇ + w (resp. ∇ - w ) is defined by using the corresponding edge derivative definition in Eq. ( 1). The norm of these vectors represent the local variation of f at a vertex of the graph. Several norms can be used but we focus on the L pnorm. Hence, the following discrete norm of the gradients in Eq. ( 2) can be considered. For any function f defined on V and for u ∈ V , when 0< p<+∞

∇ ± w f (u) p = v∼u ∂ ± v f (u) p 1/p then, ∇ + w f (u) p = v∼u w p/2 uv max 0, f (v)-f (u) p 1/p and ∇ - w f (u) p = v∼u w p/2 uv min 0, f (v)-f (u) p 1/p . (3)
When p=∞, one obtains:

∇ ± w f (u) ∞ = max v∼u ∂ ± v f (u) then, ∇ + w f (u) ∞ = max v∼u w 1/2 uv max 0, f (v)-f (u) and
∇ - w f (u) ∞ = max v∼u w 1/2 uv min 0, f (v)-f (u) . (4)
Same definitions can be provided for the gradient ∇ w f .

A FAMILY OF DILATIONS AND EROSIONS ON WEIGHTED GRAPHS

Let A be a subset of V . We denote by A + and A -respectively the outer and the inner boundary sets of A where A + = u ∈A c : ∃ v∈A, v∼u and A -= u∈A: ∃ v∈A c , v∼u where A c is the complement of A. Dilation over A is a growth process that adds vertices from A + to A. Erosion over A is a contraction process that removes vertices from A -to A.

Any function f can be decomposed into its levels sets f k =H(f -k) where H is the Heaviside function. Then, one can find a set A k ⊂V such that f k = χ A k where χ: V →{0, 1} is the indicator function. As for the continuous case, a simple variational definition of a dilation applied to f k can be interpreted as maximizing a surface gain proportional to

+ ∇ w f k (u)
p , the gradient of f k . On can demonstrate1 from Eq. ( 3) that for 0< p< + ∞,

∇ + w f k (u) p = v∼u v∈A w p/2 uv 1/p χ (A + ) k (u) and
∇ - w f k (u) p = v∼u v∈A w p/2 uv 1/p χ (A -) k (u) .
In a same way, we also have

∇ w f k (u) p = ∇ + w f k (u) p + ∇ - w f k (u) p .
Therefore, ∇ w f k p can be reduced to ∇ + w f k p for u∈(A + ) k which corresponds to a dilation on A k and can be expressed by

∂ t f k (u) = ∇ + w f k p .
By extending this to all the levels of f , one obtains a p-dilation process over a graph. Finally, with these properties, we can naturally consider the following family (parameterized by p and w) of p-dilation and p-erosion over a weighted graph,

δ p,t (f ) := ∂ t f = + ∇ + w f p and ε p,t (f ) := ∂ t f = -∇ + w f p
where t ≥ 0 corresponds to a scale parameter. By using discretization in time, and with the usual notation f n (u) ≈ f (u, n∆t), we obtain the following iterative algorithms for the p-dilation. For all u ∈ V , when 0< p< + ∞

f n+1 (u)=f n (u)+∆t v∼u w p/2 uv max 0, f n (v)-f n (u) p and when p=∞ f n+1 (u)=f n (u)+∆t max v∼u w p/2 uv max 0, f n (v)-f n (u) .
Where n corresponds to the iteration step, f 0 : V →R is the initial function defined on V , and the initial condition is f (0) (u)=f 0 (u). A similar scheme is obtained for the p-erosion. At each iteration step, the time complexity is O(|V | 2 ), where |.| stands for the cardinality of a set.

Remarks. In the particular case of a grayscale image, with a 4-adjacency grid graph associated to the image (one vertex per pixel) and a constant weight function (w=1), our approach corresponds exactly to the conventional Osher-Sethian upwind discretization scheme [START_REF] Maragos | Partial differental equations for morphological scale-spaces and eikonal applications[END_REF] when p=2. Moreover, for the case where p=∞ and ∆t=1, one recovers the classical algebraic morphological formulation overs graphs.

EXPERIMENTS

The proposed family of p-dilation and p-erosion can be applied on any function defined on a discrete data set which can be represented by a weighted graph. The application of dilation or erosion on a vector valued function f 0 : V ⊂ R q →R r where for u ∈ V , f 0 (u)= f 0 1 (u), . . . , f 0 r (u) T and f 0 i : V →R is the ith component of f 0 (u) consists in r-independent iterative dilation or erosion schemes described in Section 3. In this case, the weight function of the graph acts as a coupling term between each vector components.

Image Processing

Let f 0 :V ⊂ R2 →R be a scalar image which defines a mapping from the vertices to gray levels. Fig. 1 presents an initial image, a fine partition of this image and a reconstructed image from the partition where the pixel values of each region of the partition are replaced by the mean pixel value of its region. To obtain the fine partition, many well known methods can be used. in this work an approach based on the generalized Voronoï diagram [6] is used. The amortized time complexity to obtain such partition is O(E+V logV ) with Dikjstra algorithm and Fibonacci heap structure. Then, the partition can be associated with a Region Adjacency Graph (RAG), where vertices represent regions and where edges link adjacent regions. One can note the significant data reduction of the reconstructed version as compared to the original one (88% of reduction in terms of vertices). Fig. 2 illustrates the application of the proposed family of p-dilation and p-erosion on gray level image with different values of p, weight functions and graph topologies. In Fig. 2. For each morphological processing (p-dilation or p-erosion) a table view of the results is provided (with the same number of iterations). Each first and second rows show a local processing performed on a 4-adjacency grid graph (N 4 ) associated to the initial image with constant weights w=1 (the first row corresponds to classical approaches) and non constant weights w=w e =w uv = exp f (u)-f (v)

2 2 /σ 2 (second row) where f (u)-f (v)
2 is the L 2 -norm between feature vectors. Third row presents results with non local processing on a 48-adjacency graph (7×7 neighborhood window) with a 3×3-patch as a feature vector (denoted N 48,9 ) and non constant weights (w=w e ). Fourth and last rows present results on the Region Adjacency Graph of the image partition with constant weights (fourth row) and non constant weights w=w e (last row).

These results show that by using non constant weights, the proposed p-dilation and p-erosion preserve better preserve edges as compared to classical approaches. Once one uses a non local configuration, one also better preserve fine structures and repetitive elements. In addition, our formulation works with equal ease on graphs of the arbitrary topology and, therefore, we can apply the same schemes on any graph representing the image. For instance, we can use the above mentioned RAG (Figs. 1(b) and 1(c)) and associate to each vertex the mean gray level of its corresponding region. This exhibits similar behaviors (the two last rows in Fig. 2) while reducing complexity due to the reduced number of vertices to consider.

Manifold Processing

To illustrate the generality of our formulation, we apply our family of p-dilation and p-erosion on an unorganized data set (the Iris data set 2 ). Fig. 3 presents results of a p-dilation and a p-erosion for p=1 with a non constant weight function w uv =1/ f (u)-f (v) 2 2 . The graph used to represent data corresponds to a k-nearest neighbors graph (k=30). Iris data set contains 3 classes of samples in 4-dimensions, with 50 samples in each class. Fig. 3 shows projection of the two features (for better visualization) but it is important to note that weight function takes into account all the data features. Therefore, the proposed morphological processing of Iris data set consists in 4 independent p-erosion or p-dilation schemes. Fig. 3 shows the evolution of the processing on the data for different number n of iterations. One can note that the dilation and the erosion tend to naturally group the data in different parts of the feature space and to move data to respectively the maximum and the minimum of the feature space.

CONCLUSION

In this paper, we presented a new formalism that extends PDEs-based Mathematical Morphology operators to a discrete and non local scheme by considering partial difference equations over weighted graphs of the arbitrary topologies. We introduced basic operators: dilation and erosion that can be derived to many other morphological methods such as opening, closing, leveling, etc. Experimental results show the potential of the proposed formalism for the non local processing of images (represented as grid or region adjacency graphs) that better preserves edges and fine structures. Moreover, our formalism can be applied on manifolds and high dimensional discrete data sets that can be represented by a graph. This opens opens new application fields for morphological processing such as machine learning.

  ∂ t f = ±|∇f | where ∇= ∂x, ∂y T denotes the spatial nabla operator. Moreover, This work was partially supported under a research grant of the ANR Foundation (ANR-06-MDCA-008-01/FOGRIMMI) and a doctoral grant of the Conseil Régional de Basse-Normandie and of the Coeur et Cancer association in collaboration with the Department of Anatomical and Cytological Pathology from Cotentin Hospital Center.

Fig. 1 .

 1 Fig. 1. (a) Original image of size 256×256=65 536 pixels. (b) fine partition of (a) having 8 114 regions (i.e. 88% of reduction). (c) reconstructed image from (b).

Fig. 2 .

 2 Fig. 2. p-dilation and p-erosion on the image of Fig. 1(a) with different values of p, weight functions and graph topologies. See text for more details.

Proof can be obtained by studying cases where u∈A k or u / ∈A k and similarly for v∼u.

UCI Machine Learning Repository, http://www.ics.uci.edu/ ∼ mlearn/MLRepository.html

[5] A. Elmoataz, O. Lézoray, and S. Bougleux, "Nonlocal discrete regularization an weighted graphs: a framework for image and manifolds processing," HAL TR, Université de Caen Basse-Normandie, 2007.

[6] P. Arbeláez and L. Cohen, "Energy partitions and image segmentation," JMIV, vol. 20, no. 1-2, pp. 43-57, 2004. p-dilation p-erosion


Philippe Swartvagher
email: philippe.swartvagher@inria.fr

Interferences between Communications and Computations in Distributed HPC Systems

Keywords: HPC, MPI, Memory Contention

Overlapping communications with computations in distributed applications should increase their performances and allow to reach better scalability. This implies, by construction, communications are executed in parallel of computations. In this work, we explore the impact of computations on communication performances and vice-versa, with a focus on the role of memory contention. One main observation is that highly memory-bound computations can have a severe impact on network bandwidth.

Introduction

Doing in parallel communications and computations (with non-blocking MPI calls or even more complex systems such as task-based runtime systems) is an increasing trend to get higher performances. It has been observed [START_REF] Langguth | Memory Bandwidth Contention: Communication vs Computation Tradeoffs in Supercomputers with Multicore Architectures[END_REF][START_REF] Groves | NiMC: Characterizing and Eliminating Network-Induced Memory Contention[END_REF], that, sometimes, when computations and communications are executed side by side, communications are slower than nominal performances and computations can also be degraded.

Since possible interactions between communications and computations, and especially the impact on communication performances, are not well detailed in the literature (but only mentioned), we propose in this work to study the possible causes of these interferences and measure their impact on both communication and computing performances. Reported observations here deal with the impact of memory contention. Since we target HPC systems, we consider only fast networks (here, InfiniBand) as well as inter-node communications.

We describe here our experimental protocol, provide first results and then discuss future work.

Advisers: Alexandre Denis and Emmanuel Jeannot Inria Bordeaux -Sud-Ouest, 33405 Talence, France {alexandre.denis,emmanuel.jeannot}@inria.fr

Methodology

Our goal is to measure performances of communications and computations when they are run side by side. To achieve this, we have designed a multithreaded and parallel benchmark using MPI+OpenMP. Per MPI process, one thread is dedicated to communications (it submits communication instructions and ensures MPI progression) and other threads (driven by OpenMP) do computations. This configuration mimics how some runtime systems work. The communication benchmark performs ping-pongs to measure network latency and bandwidth.

Benchmarking protocol

We need to compare performances of communications and computations when they are executed alone and when they are executed together. Therefore we decomposed our benchmark into the following steps: Regarding performances of communications, we use two metrics: latency by exchanging 4 bytes of data (one float), and bandwidth evaluated for 64 MB message size.

According to which effect we want to observe, instructions done by computing threads will be different. Computations have to be embarrassingly parallel to always get the maximum level of parallelism and avoid an overhead caused by scheduling and dependencies. Moreover we use weak scalability to easily change the number of computing cores (and thus change the level of parallelism).

To see the impact of memory contention caused by data used for computations and data used for communications, we generate memory contention by computing cores doing memory-bound kernels:

COPY (b[i] ← a[i]) and TRIAD (c[i] ← a[i] + C × b[i]
) from the STREAM benchmark suite [START_REF] Mccalpin | Memory bandwidth and machine balance in high performance computers[END_REF]. Moreover, to really produce memory contention, we allocate memory on a single NUMA node, to increase the traffic on the memory bus between cores belonging to different NUMA nodes. The performance of the computing benchmark is measured using the memory bandwidth per core (hence higher is better).

Experimental environment

We ran our own benchmark suite 1 on several clusters with different characteristics: from small experimental clusters to large production ones. We present here the results obtained on henri nodes which are dual Intel Xeon Gold 6140 at 2.3 GHz with 36 cores split across 4 NUMA nodes and 96 GB of RAM and equipped with InfiniBand ConnectX-4 EDR. henri nodes run Linux 5.7.7 with Debian 10, gcc 9.3. We show results obtained with MadMPI, the MPI interface of NewMadeleine [START_REF] Aumage | NewMadeleine: a Fast Communication Scheduling Engine for High Performance Networks[END_REF] ; we observed similar results with other MPI implementations, such as OpenMPI 4.0.

First Observations

Our first hypothesis was that data accessed for computations and for communications may create contention on the memory bus, thus reducing performances. To test this hypothesis, we applied protocol described in previous section. Fig. 1 shows that network latency is impacted by the STREAM operations when there are at least 22 computing cores and this impact can double the regular latency when all available cores are computing. However, STREAM operations are not impacted by the ping-pong benchmark. The network bandwidth is impacted sooner, from 3 computing cores (Fig. 2). When all available cores are computing, the network bandwidth is reduced by almost two third from the regular network bandwidth. Memory bandwidth measured by the STREAM benchmark is lower when network bandwidth is measured at the same time as when network latency is measured, which is expected because one bandwidth ping-pong transfers more data than a latency ping-pong (64 MB vs 4 B).

Conclusion

Overlapping communications with computations is a well-known technique to increase performance of applications. However, our first experiments report memory contention caused by data used for communications and for computations can reduce performances of communications and computations.

The rest of our study explores parameters impacting interferences between communications and computations such as the placement of data and threads, the size of the data transmitted through the network and the arithmetic intensity of instructions executed by computing threads. We also measure the effect of the variations of processor frequencies and the effect of using a task-based runtime system.

As future works, we would like to better understand origins of these interferences to predict and quantify them. We would also like to measure the impact of data movements between main memory and GPUs. This knowledge could then be used in runtime systems to better predict performances and thus better schedules tasks to minimize the impacts of the interferences between computations and communications.

1 .

 1 Computation without communication 2. Communication without computation 3. Computation with side by side communication Computations and communications use different data and hence are completely independent. Plots to represent results compare performance of communications and computations when they are executed separately or simultaneously. The former are represented by plain curves and the later by dashed curves.

Fig. 1 .

 1 Fig. 1. Network latency (for 4 B) with memory-bound computations in parallel

Fig. 2 .

 2 Fig. 2. Network bandwidth (for 64 MB) with memory-bound computations in parallel

Available on https://gitlab.inria.fr/pswartva/memory-contention