Member, IEEE Mohammadhadi Shateri

Francisco Messina

Senior Member, IEEE Pablo Piantanida

Fabrice Labeau

Privacy-Cost Management in Smart Meters With Mutual-Information-Based Reinforcement Learning

Keywords: Internet of Things, Cyber-Physical System, Smart meters privacy, Mutual information, Deep reinforcement learning, Deep double Q-learning

published or not. The documents may come   L'archive ouverte pluridisciplinaire

I. INTRODUCTION

With recent developments in computation, communication, and control technologies, the prevalent Internet of Things (IoT) facilitated the connection of billions of smart devices and sensors which can collect and exchange data in real-time [START_REF] Mohammadi | Deep learning for iot big data and streaming analytics: A survey[END_REF]. This has raised a lot of interest in Cyber-Physical Systems (CPS) for different applications such as Smart Grids (SGs) [START_REF] Humayed | Cyber-physical systems security:a survey[END_REF]. Smart Meters (SMs) are a key component of the socalled advanced metering infrastructure (AMI), which is a critical subsystem of SGs [START_REF] Mohassel | A survey on advanced metering infrastructure[END_REF]. SMs are capable of measuring electricity consumption of users at a fine-grained level and share it with the Utility Provider (UP) in almost real-time. This large amount of data provides immense opportunities for both customers and operators, leading to the emergence of M. Shateri and F. Labeau are with the Department of Electrical and Computer Engineering, McGill University, QC, Canada (e-mail: mohammadhadi.shateri@mail.mcgill.ca ; fabrice.labeau@mcgill.ca). F. Messina is with the School of Engineering -Universidad de Buenos Aires, Buenos Aires, Argentina (e-mail: fmessina@fi.uba.ar). P. Piantanida is with Laboratoire des Signaux et Systèmes (L2S), Centrale-Supélec CNRS Université Paris-Saclay, 91190 Gif-sur-Yvette, France (e-mail: pablo.piantanida@centralesupelec.fr).

the new field of SMs data analytics [START_REF] Wang | Review of smart meter data analytics: Applications, methodologies, and challenges[END_REF]. However, SMs data also contain sensitive information about users which could easily be inferred by malicious third-parties or attackers if no preventive measures are taken. For instance, an eavesdropper can apply non-intrusive load monitoring (NILM) methods and deep learning approaches to infer the user's presence at home [START_REF] Giaconi | Privacy-aware smart metering: Progress and challenges[END_REF], [START_REF] Feng | Deep learning-based real-time building occupancy detection using ami data[END_REF] and even the types of appliances being used at a specific time [START_REF] Molina-Markham | Private memoirs of a smart meter[END_REF]. Therefore, the massive deployment and adoption of SMs necessitate the development of efficient privacy-aware strategies for real-time data sharing in order to keep the users' sensitive information private against potential attackers. It should be noted that the privacy issue regarding the SMs data sharing is distinct from SMs data security in terms of the legitimate users and the adversaries. Unlike the data security using encryption methods, in the SM data privacy any legitimate receiver of the data including utility provider can be considered at the same time as a potential malicious attacker. Therefore, the traditional encryption techniques would be ineffective in providing privacy against utility provider [START_REF] Giaconi | Privacy-aware smart metering: Progress and challenges[END_REF]. Recently, reinforcement learning has been used in designing models for preserving privacy in SM data sharing where physical resources at the users' end are used to prevent exposing consumers' privacy and behaviour to the UP [START_REF] Shateri | Privacy-cost management in smart meters using deep reinforcement learning[END_REF]- [START_REF] Sun | Smart meter privacy: Exploiting the potential of household energy storage units[END_REF]. The main idea of reinforcement learning is based on the interaction between an agent (here a privacy-cost management unit) and the environment (which includes user and physical resources). More specifically, at each state imposed by the environment, the agent would take an action for the sake of receiving a maximum (future) reward and then being placed in the next state by the environment. These interactions can be used by the agent to learn the optimal policy which maximized the total reward [START_REF] Sutton | Reinforcement Learning: An Introduction[END_REF]. More details on the reinforcement learning and its formulation will be discussed in the Section II. Still a framework incorporates a generic privacy metric such as mutual information, is missing since it can add several challenges which requires some amendments in the traditional algorithms.

A. Related work

A substantial amount of studies on SMs privacy were conducted, which can be classified in two main families: (i) SMs data manipulation techniques [START_REF] Efthymiou | Smart grid privacy via anonymization of smart metering data[END_REF]- [START_REF] Gough | Preserving privacy of smart meter data in a smart grid environment[END_REF]; and (ii) user's demand load shaping approaches [START_REF] Shateri | Privacy-cost management in smart meters using deep reinforcement learning[END_REF]- [START_REF] Sun | Smart meter privacy: Exploiting the potential of household energy storage units[END_REF], [START_REF] Kalogridis | Privacy for smart meters: Towards undetectable appliance load signatures[END_REF]- [START_REF] Giaconi | Optimal demand-side management for joint privacy-cost optimization with energy storage[END_REF]. On the one hand, in the first family of methods, the consumers' load data are manipulated by a noisy transformation before sharing it to the UP. In this setting, there is a natural trade-off between the distortion or utility of the data and privacy guarantees. On the other hand, in the second family of privacy-aware techniques, the actual electricity consumption of the users, as seen by the grid (i.e., the grid load), are shaped using a combination of different physical resources such as Rechargeable Batteries (RBs), Electric Vehicles (EVs), Heating, Ventilation, and Air Conditioning (HVAC) units, and Renewable Energy Sources (RES). Note that, in this scenario, the consumer load is different than the grid load (e.g., the grid load can be higher than the consumer load provided that an RB is being charged by the user). The goal of these methods is to mask the consumer load but, at the same time, efficiently make use of the available resources considering their physical constraints and wear and tear, as well as a possibly time-varying electricity rate. In this framework, there is generally a trade-off between the overall electricity expenses and privacy guarantees.

In some recent studies, physical resources are employed to minimize the average relative difference between the grid load and a constant target load [START_REF] Shateri | Privacy-cost management in smart meters using deep reinforcement learning[END_REF], [START_REF] Sun | Smart meter privacy: Exploiting the potential of household energy storage units[END_REF], [START_REF] Giaconi | Optimal demand-side management for joint privacy-cost optimization with energy storage[END_REF], i.e., to flatten the electricity consumption reported by the SMs. In [START_REF] Shateri | Privacy-cost management in smart meters using deep reinforcement learning[END_REF], following the formulation in [START_REF] Sun | Smart meter privacy: Exploiting the potential of household energy storage units[END_REF], the SMs privacy problem is cast as a Markov Decision Process (MDP) and a model-free (i.e., not assuming full knowledge of the environment dynamics of the MDP) Deep Reinforcement Learning (DRL) algorithm known as the Deep Double Q-Learning (DDQL) method is used to tackle the problem. Even though this framework has been shown to be useful in limiting the leakage of sensitive information, the effectiveness of the flatness-based privacy measure remains unclear.

A formal privacy measure from information theory known as the mutual information (MI), between the user's demand and reported grid load, was proposed in [START_REF] Sankar | Smart meter privacy: A theoretical framework[END_REF] and has since been adopted in several works [START_REF] Li | Information-theoretic privacy for smart metering systems with a rechargeable battery[END_REF], [START_REF] Tan | Increasing smart meter privacy through energy harvesting and storage devices[END_REF], [START_REF] Gomez-Vilardebo | Smart meter privacy for multiple users in the presence of an alternative energy source[END_REF], [START_REF] Giaconi | Smart meter privacy with renewable energy and an energy storage device[END_REF], [START_REF] You | Energy management strategy for smart meter privacy and cost saving[END_REF]. In [START_REF] Li | Information-theoretic privacy for smart metering systems with a rechargeable battery[END_REF], the SMs information-theoretic privacy problem is formulated as an MDP and the optimum policy is obtained numerically using dynamic programming, assuming full knowledge of the MDP. For a simplified scenario in which the demand load is assumed to be known and independent and identically distributed (i.i.d.), a single-letter expression for the average information leakage was characterized. However, these methods may not be directly applicable in practice since the MDP is generally not fully-known due to the unknown dynamics of the demand load. Furthermore, the electricity cost is not part of the formulation and thus, the cost-privacy trade-off is not taken into account.

B. Contributions

In this paper, we adopt a demand load shaping privacypreserving strategy. As can be seen from the previous discussion, there is a clear gap between the DRL line of research and the information-theoretic work on this framework. As a matter of fact, the DRL methods rely on weak privacy measures, such as the flatness measure, which offer no statistical privacy guarantees (see Section III-A for further discussion). Whereas, information-theoretic approaches suggest the use of formal privacy measures, such as MI, with strong statistical privacy guarantees (see Section III-B for further details). The main goal of this paper, which further extends [START_REF] Shateri | Privacy-cost management in smart meters using deep reinforcement learning[END_REF], is to develop a new DRL algorithm with the advantage of being model-free but also incorporating a strong information-theoretic privacy measure. Concretely, the main contribution is to incorporate the MI between the user's demand and grid load as a privacy measure within the DDQL algorithm framework. As will be shown, this introduces the challenge of estimating the privacy signal to feed into the DRL agent during the training phase, which is overcome by adding a new neural network (referred to as the H-network). Interestingly, based on the structure of this network, we can either use a general MI privacy measure or a simplified MI privacy measure based on a strong i.i.d. assumption. The latter case is an important benchmark as it allows us to quantify the importance of considering time correlation within the privacy measure. We then study the empirical information-leakage rate versus electricity cost trade-off on a real SMs dataset, and compare the MI-based and flatness-based privacy measures, showing the advantages of our approach. Finally, the performance of the DDQL-MI algorithm is assessed in two practical scenarios: (a) an attacker aiming to infer the actual consumer load; and (b) an attacker trying to infer the house occupancy status.

Organization of the paper. The rest of the paper is organized as follows. In Section II, we review the MDP formulation of the privacy-aware demand shaping problem. Then, in Section III, we review two commonly used privacy measures and discuss how to use MI as a privacy measure in the MDP framework. We also discuss why MI is a superior privacy measure as compared with the flatness-based one. The novel DDQL-MI algorithm is then presented in Section IV. The numerical performance of the DDQL-MI algorithm using actual SMs data is studied in Section V, where we show the impact of the privacy measure choice on the results. Some concluding remarks close the paper in Section VI.

Notation and conventions. We use capital letter to denote random variables and lowercase to denote specific values. P (y) is the probability distribution of random variable Y ; E[•] is expectation with respect to the joint distribution of all random variables involved, E[X|Y ] is the conditional expectation of X given Y , E ⇡ [•] is the expectation for a given agent's policy ⇡ (see [START_REF] Sutton | Reinforcement Learning: An Introduction[END_REF] for details); H(Y ) is the entropy of Y and I(X; Y ) = H(X) H(X|Y ) is the MI between X and Y , where H(X|Y ) is the conditional entropy of X given Y (see [START_REF] Cover | Elements of information theory[END_REF] for details).

II. PROBLEM FORMULATION

A. Demand shaping using physical resources

Consider the smart metering system represented in Fig. 1 where an intelligent agent, named as Privacy-Cost Management Unit (PCMU), is used to hide the household demand load using an RB while keeping the total electricity cost minimum. It should be noted that other physical resources can also be readily incorporated in our framework, as shown in other works [START_REF] Sun | Smart meter privacy: Exploiting the potential of household energy storage units[END_REF], [START_REF] Erdemir | Privacy-cost trade-off in a smart meter system with a renewable energy source and a rechargeable battery[END_REF].

Let Y t be the random variable denoting the consumer's demand load/power, i.e., the total power demanded by the appliances at time t, and let Z t be the load received from the grid, and L t 2 [0, 1] the level of charge of battery (normalized by the capacity of the battery), where t 2 T = {1, . . . , T }. Given the demand load Y t and the level of charge available in the battery L t at time t, the PCMU needs to determine the optimal charging/discharging rate of the battery B t to physically distort the actual demand load, so that the grid load Z t , given by Z t = Y t + B t does not reveal information about the user's demand load Y t . The goal of the PCMU is to limit the performance of a potential attacker trying to violate the user's privacy by inferring sensitive information (which could be either the actual demand load Y t or a correlated variable of interest to the attacker) from Z t . To make this precise, we first introduce the MDP formulation of the problem.

B. Markov Decision Process (MDP) model

Following previous studies [START_REF] Li | Information-theoretic privacy for smart metering systems with a rechargeable battery[END_REF], [START_REF] Sun | Smart meter privacy: Exploiting the potential of household energy storage units[END_REF], [START_REF] Erdemir | Privacy-cost trade-off in a smart meter system with a renewable energy source and a rechargeable battery[END_REF], the problem of finding the optimal policy for the PCMU is formulated as a Markov Decision Process (MDP) to capture the agentenvironment interaction (see Fig. 2). An MDP is determined [START_REF] Sutton | Reinforcement Learning: An Introduction[END_REF] by the tuple S, A(s), P (s t+1 |s t , a t ), r(s t , a t ), :

• State space S, which determines all the possible states that the agent could be in; • Action space A(s), which determines the feasible actions the agent can take at state s 2 S; • Environment dynamics P (s t+1 |s t , a t ), which gives the probability of the state evolving to S t+1 = s t+1 when the current state is S t = s t and the agent takes the action

A t = a t ; • Reward function r(s t , a t )
, which is the immediate reward obtained due to taking action A t = a t at state S t = s t ; • Discount factor 2 [0, 1], which is the decay constant of future rewards and therefore determines their importance to the agent. In our setting, we assume a fixed finite horizon, i.e., T < 1 is a constant. Therefore, the discount factor can be considered as = 1 without any convergence issues [START_REF] Sutton | Reinforcement Learning: An Introduction[END_REF]. Moreover, as it will become apparent in Section III-B, this is the most natural choice for our problem. In general, the PCMU starts from an initial state S 1 , and by following a policy ⇡(a|s) = P (a|s), takes the first action A 1 .

As a consequence of the action A 1 , the PCMU receives the reward R 1 from its (artificial) environment and transitions to the state S 2 (see Fig. 2). Thereby, this gives rise to a trajectory of states, actions and rewards [S 1 , A 1 , R 1 , . . . , S T , A T , R T ], which is referred to as an episode. Fig. 2. Agent-Environment interaction in an MDP [START_REF] Sutton | Reinforcement Learning: An Introduction[END_REF].

One general approach to find the optimum policy is via the action-value function Q : S ⇥ A ! R, which is defined as follows [START_REF] Sutton | Reinforcement Learning: An Introduction[END_REF]:

Q(s, a) = E ⇡ 2 6 4 T X t=1 t 1 r (S t , A t ) S 1 = s, A 1 = a 3 7 5 , (1)
where the subscript in the expectation emphasizes that it is a function of the policy ⇡ chosen by the agent. The action-value function represents the expected cumulative reward received by the agent when it starts at some state s 2 S, takes an action a 2 A(s), and then follows the policy ⇡. Therefore, an optimal policy ⇡ ⇤ can be obtained by maximizing Q(s, a) over all possible policies for all pairs (s, a) 2 S ⇥ A(s). We use Q ⇤ to denote the optimal state-action value function, that is, Q ⇤ (s, a) = max ⇡ Q(s, a). It is well-known that, for any MDP, there always exists at least one optimal policy which is deterministic. Assuming uniqueness for the sake of simplicity of presentation, we can write the optimal policy simply as

⇡ ⇤ (a|s) = 1 if a = arg max a 0 2A(s) Q ⇤ (s, a 0 ) and ⇡ ⇤ (a|s) = 0 otherwise.
In our framework, the state at time t is defined as

S t = [L t , Y t ] T ,
L t+1 = L t + B t t ⌘ C , (2)
where B t 2 [b min , b max ] (b min and b max are the minimum and maximum charging/discharging rate of the RB), L t 2 [l min , l max ] (l min and l max are the minimum and maximum level of charge of the RB), t is the load sampling rate, ⌘ is the charging/discharging efficiency factor of the RB, and C is the capacity of the RB. It should be noted that P y t+1 |y t is unknown and, in general, it is difficult to estimate accurately [START_REF] Sun | Smart meter privacy: Exploiting the potential of household energy storage units[END_REF]. To get rid of this issue, the focus of this study is on model-free DRL algorithms which do not require full knowledge of the environment dynamics. Finally, we need to define the reward function, which will guide the agent to learn an optimal policy. For our purposes, it should be an appropriate combination of the privacy leakage and the associated electricity cost. Following [START_REF] Sun | Smart meter privacy: Exploiting the potential of household energy storage units[END_REF], the reward function is inversely interpreted as a loss function: `(s t , a t ) = r (s t , a t ). Assuming a privacy leakage signal f (s t , a t ) and an electricity cost signal g(s t , a t ), which will be defined later, the one-step loss function can be defined as follows:

`(s t , a t ) = r (s t , a t ) = g(s t , a t ) + (1 )f (s t , a t ), (3)
where 2 [0, 1] controls the privacy-cost trade-off. Concretely, for = 0 the goal of the agent will be to minimize the expected cumulative privacy leakage signal, while for = 1 it will be to minimize the expected cumulative energy cost. As is reduced, the PCMU should be able to provide the consumer a higher privacy level but at the expense of an increase in the energy cost. Studying this trade-off and its practical implications is of fundamental importance to properly design the PCMU. Notice that, using ( 1) and ( 3), we can decompose the action-value function as follows:

Q(s, a) = Q c (s, a) + (1 )Q p (s, a), (4)
where

Q c (s, a) = E ⇡ [ P T t=1 t 1 g(S t , A t ) |S 1 = s, A 1 = a] and Q p (s, a) = E ⇡ [ P T t=1 t 1 f (S t , A t ) |S 1 = s, A 1 = a]
are the cost and privacy action-value functions, respectively. Therefore, the objective of the PCMU, considering the aforementioned constraints on the capacity of the battery and its charging/discharging rate, is to find an optimal policy ⇡ ⇤ that maximizes the expected total reward in equation [START_REF] Wang | Review of smart meter data analytics: Applications, methodologies, and challenges[END_REF]. It should also be noted that the problem constraints are incorporated in the definitions of the state and action spaces.

The total cost associated with this privacy-aware framework can be due to the electricity cost and also to the cost related to the battery wear and tear. Considering for simplicity that no energy can be sold to the grid by the users, the electricity cost at time t can be computed as

C t = t h t [Z t ]
+ where h t is the price of purchasing 1 kWh of energy from the grid at time t and

[Z t ] + = max(Z t , 0). Since [Z t ] +  |Z t |  Y t +|B t | , and
Y t is not controlled by the PCMU, we consider the following electricity cost signal:

g(S t , A t ) = t h t |B t |. (5)
Notice that C t  t h t (Y t + |B t |), so this electricity cost signal effectively limits the actual electricity cost. In addition, this definition of g(S t , A t ) incidentally takes into account the battery wear and tear cost, since it grows as the battery use increases. The design of the privacy leakage signal is discussed in detail in the next section.

III. PRIVACY MEASURES: FLATNESS AND MUTUAL INFORMATION

In the following, two different privacy measures will be reviewed and discussed: the flatness privacy measure and the MI privacy measure. A. Flatness privacy measure

In [START_REF] Sun | Smart meter privacy: Exploiting the potential of household energy storage units[END_REF], the privacy leakage signal received by the Reinforcement Learning (RL) agent when at state S t and taking action A t is defined as

f (S t , A t ) = Z t c c , ( 6
)
where c is a constant target level. This encourages the agent to take actions such that, on average, Z t is as close as possible to c . Although simple to compute and intuitively appealing, the problem with this quantity is that it does not capture the statistical dependence between Y T and Z T and therefore fails to be a satisfactory privacy measure. As an example to illustrate this fact, consider a PCMU with the following strategy: Z t = Y t + c for all t 2 T , where 2 R is a constant. In this case, the privacy action-value function, defined in (4), is

Q p (s, a) = E 2 6 4 T X t=1 t 1 | | |Y t | | c | s 1 = s, a 1 = a 3 7 5 / | |, (7)
that is, Q p (s, a) is proportional to | |. If = 0, we have Q p (s, a) = 0 as expected. In such scenario, Z T is constant and it does not provide any information to infer the actual value of Y T . In other words, full privacy is achieved. However, for any 6 = 0, the variables Y t and Z t are maximally correlated (i.e., the correlation coefficient between Z t and Y t is either 1 or -1 depending on the sign of ) and the task of inferring Y t from Z t is trivial (assuming the attacker is able to estimate only two parameters). Therefore, all these cases can be considered equivalent from the privacy point of view, but Q p (s, a) can take any value in ( 1, 0) as | | is modified. Thus, the flatness privacy measure can be completely misleading in some scenarios. This problem is illustrated in Fig. 3.

B. Mutual information privacy measure

1) General Case: A much stronger privacy measure, proposed in [START_REF] Sankar | Smart meter privacy: A theoretical framework[END_REF] for the SMs privacy problem, and used since in several works [START_REF] Li | Information-theoretic privacy for smart metering systems with a rechargeable battery[END_REF], [START_REF] Tan | Increasing smart meter privacy through energy harvesting and storage devices[END_REF], [START_REF] Gomez-Vilardebo | Smart meter privacy for multiple users in the presence of an alternative energy source[END_REF], [START_REF] Giaconi | Smart meter privacy with renewable energy and an energy storage device[END_REF], is the MI between the demand load and the grid load, which is defined as follows:

I(Y T ; Z T ) = E " log p(Y T , Z T ) p(Y T )p(Z T ) # . (8)
Intuitively, MI measures the degree of dependence between

Z T and
I(Y T ; Z T ) = ( 0 if = 0, H(Y T ) otherwise. ( 9
)
This means that full privacy is achieved only for the case = 0 and other choices of lead to a maximal information leakage. Recall from (3) that, in order to define the loss or the reward signal, we need to find the instantaneous (random) privacy leakage signal received by the agent when it is at state s t and takes action a t . To do so, we first define the privacy action-value function as follows:

Q p (s, a) = I(Y T ; Z T ) = T X t=1 H(Y t |Y t 1 , Z T ) H(Y T ). (10)
It should be noted that Q p (s, a) is a function of the policy of the agent ⇡ but not of the initial state and action. Notice that H(Y T ) is a constant, independent of the PCMU strategy, so the second term in (10) can be discarded. On the other hand, using the definition of conditional entropy and the law of total expectation, the term inside the summation of the first term can be written as follows:

H(Y t |Y t 1 , Z T ) = E[E[ log P (Y t |Y t 1 , Z T )|S t , A t ]]. (11)
Notice that the inner conditional expectation is an explicit function of S t and A t , as required. Therefore, we define the privacy leakage signal as follows:

f (s t , a t ) = E[log P (Y t |Y t 1 , Z T )|S t = s t , A t = a t ]. ( 12
)
With this definition, the (negative) expected cumulative privacy leakage signal over an episode is equal to the MI in [START_REF] Sun | Smart meter privacy: Exploiting the potential of household energy storage units[END_REF] up to an additive constant. Notice that to recover the MI exactly we need to set = 1, i.e., do not discount future rewards. Therefore, we do not use discounting in this paper. As it was explained in Section II-B, this choice is in fact typical for finite horizon problems. Implementing the privacy leakage signal proposed in equation ( 12) poses three main challenges. First of all, to approximate this privacy signal, we need to estimate the unknown conditional distributions P (y t |y t 1 , z T ) for each t 2 T . Secondly, although the expectation operation appearing in [START_REF] Efthymiou | Smart grid privacy via anonymization of smart metering data[END_REF] can be estimated based on previous experiences of the agent using a Monte Carlo approach, it needs enough samples for each possible pairs of (s t , a t ) and so requires a huge buffer. Finally, unlike equation ( 6), the approximation of this privacy measure is non-causal as it involves the whole sequence of the grid load z T at each t. In Section IV we will discuss how to deal with these challenges.

2) I.I.D. case: In order to understand the role of the correlation in time of the time series Y T and Z T , we also consider, as a benchmark, the case in which Y T is assumed to be independent and identically distributed (i.i.d.) and we model the transformation between Y t and Z t as memoryless but arbitrary, i.e., Z t = g(Y t ) where g is a fixed random transformation. Although this assumption clearly does not hold in our MDP framework nor in practice, it is interesting to study this scenario to assess the advantage of taking into account the correlations across time for the privacy measure computation. Note that, in this case, Z T is also i.i.d. In such a case, it can be shown that:

I(Y T ; Z T ) = T X t=1 I(Y t ; Z t ) = T [H(Y ) H(Y |Z)], (13)
where we have omitted the time dependence since all the pairs (Y t , Z t ) are assumed to be i.i.d. and therefore all the terms are equal. Similarly as before, noting the relation

H(Y |Z) = E[E[ log P (Y |Z)|S, A]], (14)
we can define the privacy signal simply as

f (s, a) = E[log P (Y |Z)|S = s, A = a]. (15)
Note that, this i.i.d. case is not guaranteed to provide any control over the value of I(Y T ; Z T ) in the general case. In fact, by using the standard properties of MI, it can be shown that [START_REF] Cover | Elements of information theory[END_REF]

I(Y T ; Z T ) 1 T T X t=1 I(Y t ; Z t ). (16)
In summary, the i.i.d. assumption leads to using a lower bound of the MI as a privacy measure and, therefore, does not offer real privacy guarantees.

IV. METHODOLOGY AND ALGORITHM

We first review classical RL and DRL algorithms which were used in previous works on this topic. Then, we extend the DRL algorithm to accommodate the MI privacy measure.

A. Review of CQL algorithm

The classical Q-Learning (CQL) algorithm is a simple method to learn the optimal state-action value function Q ⇤ by updating the action-value of the experienced state-action pairs. The algorithm can be summarized by the update equation:

Q (S t , A t ) = ↵  r (S t , A t ) + max a2A(St+1) Q (S t+1 , a) Q (S t , A t ) , ( 17
)
where ↵ is the step size parameter. Details on the training process and convergence properties of the CQL method can be found in [START_REF] Sutton | Reinforcement Learning: An Introduction[END_REF]. The CQL algorithm was used for smart meter privacy in [START_REF] Sun | Smart meter privacy: Exploiting the potential of household energy storage units[END_REF]. However, the main drawback of this method is that it needs to visit all the state-action pairs several times to provide a good approximation of Q ⇤ . Therefore, for large MDPs with many states and actions, convergence is usually very slow. This is the case in our problem if the action and state spaces are discretized with high accuracy.

B. DDQL-MI algorithm

To solve the slow convergence problem of the CQL algorithm, the Q-function can be approximated by using a Deep Neural Network (DNN) to generalize between different states and actions. These new methods, where deep learning is used for approximating the Q-function, are called Deep Q-Learning (DQL) methods [START_REF] Mnih | Human-level control through deep reinforcement learning[END_REF], [START_REF] Franc ¸ois-Lavet | An introduction to deep reinforcement learning[END_REF]. A general diagram presenting the agent-environment interaction under the DQL paradigm and our context is shown in Fig. 4.

The main idea of the DQL method is to approximate Q ⇤ (s, a) using a DNN called the Q-network. The Q-network takes the state s 2 S at the input and generates Q(s, a; ✓) at the output for all different actions a 2 A(s), where ✓ are the parameters of the DNN. To define the objective function for this Q-network, we observe from ( 17) that convergence is obtained when the quantity in parenthesis is equal to zero. The term r (s t , a t ) + max a2A(st+1) Q (s t+1 , a; ✓) can be interpreted as the target, while the term Q (s t , a t ; ✓) is the output of the Q-network. Thus, the mean squared error loss between target and output can be used as the loss function for training the Q-network. However, using the same network to compute the target and output often leads to instability [START_REF] Mnih | Human-level control through deep reinforcement learning[END_REF]. To address this issue, the so-called Double Q-Learning (DDQL) algorithm was proposed in [START_REF] Hasselt | Double q-learning[END_REF] and extended to the deep learning setting in [START_REF] Van Hasselt | Deep reinforcement learning with double q-learning[END_REF]. In the DDQL algorithm, a second network called the target-network (with parameters ✓ 0 ) is used to calculate the target term. The target-network parameters ✓ 0 are periodically updated by simply copying the parameters from the Q-network. Thus, using the targetnetwork, the objective function of the Q-network can be written as follows:

L QN (✓) =E ✓ r(S t , A t ) + max a2A(St+1) Q(S t+1 , a; ✓ 0 ) Q(S t , A t ; ✓) ◆ 2 . ( 18
)
It should be noted that the expectation is approximated by a Monte Carlo approach based on batches of samples (s t , a t , r(s t , a t ), s t+1 ) selected randomly from a replay buffer [START_REF] Franc ¸ois-Lavet | An introduction to deep reinforcement learning[END_REF]. The DDQL algorithm was applied for the first time in the SM privacy problem in [START_REF] Shateri | Privacy-cost management in smart meters using deep reinforcement learning[END_REF] using the flatness privacy signal defined in [START_REF] Feng | Deep learning-based real-time building occupancy detection using ami data[END_REF], showing clear improvements in performance and convergence speed with respect to the CQL algorithm.

From now on, we will refer to this method as Model I.

For the DDQL-MI method, we have to deal with the three challenges discussed in Section III-B1 regarding the privacy leakage signal presented in equation [START_REF] Efthymiou | Smart grid privacy via anonymization of smart metering data[END_REF]. To overcome the first challenge, a helper neural network named as H-network is included in the DDQL to estimate the conditional probability distributions p(y t |y t 1 , z T ), t 2 T . For the second challenge, instead of storing samples for each possible pairs of (s t , a t ) which leads to a huge storage and a very slow training, we approximate the privacy measure in equation ( 12) with its expected value over the joint distribution of (s t , a t ) , i.e., equation [START_REF] Sutton | Reinforcement Learning: An Introduction[END_REF]. Then, for each episode, the pair (y T , z T ) is stored in a second replay buffer and, when required, samples are selected randomly to approximate equation [START_REF] Sutton | Reinforcement Learning: An Introduction[END_REF]. Finally, since training the agent in DDQL is done offline, the third challenge is not an issue for its implementation. It should be emphasized at this point that, once a policy is learned, the agent will act according to it in a fully causal manner. In order to address the exploration exploitation dilemma, we adopt the ✏ greedy method: with probability 1 ✏ the action with maximum action-value Q is selected (exploitation), while with probability ✏ a random action is selected (exploration). For more details, the reader is referred to [START_REF] Sutton | Reinforcement Learning: An Introduction[END_REF]. The training of the DDQL-MI method is presented in Algorithm 1 below. It will be referred to as Model II in the following. Set the initial state s 1 = [l 1 , y 1 ].

4:

for t = 1, . . . , T do

Observe the state st = [lt, yt].

6:

Select a feasible action at using the ✏ greedy algorithm.

7:

Calculate reward r(st, at) from equation (3).

8:

Update the next state s t+1 based on (2) and observing y t+1 . 9:

Import (st, at, r(st, at), s t+1 ) into the replay buffer I.

10:

Every k 0 time-step, update the Q-network parameters using samples from the replay buffer I.

11:

Every k time-step, update the target-network by copying from Q-network parameters.

12:

Every k time-step, update the H-network and then estimate privacy measure in [START_REF] Sutton | Reinforcement Learning: An Introduction[END_REF] for all t 2 T = {1, . . . , T } using samples from the replay buffer II.

13:

end for 14:

Import (z T , y T ) into the replay buffer II. 15: end for Note: The copy step k and training step k 0 are hyperparameters.

It is worth to mention that the Q-learning algorithm has some shortcomings. First of all, for the sake of convergence, it needs a large amount of episodes to ensure all the state-action pairs are experienced multiple times. Secondly, the state and action sets need to be finite, while in practical applications it is not always the case [START_REF] Sutton | Reinforcement Learning: An Introduction[END_REF]. Finally, as we discussed in Section III, a significant limitation of previous state-of-the-art Q-learning approaches to the smart meter privacy problem is the choice of the reward function, which does not capture a strong statistical notion of privacy. These shortcomings were considered in the design of Algorithm 1. Concretely, regarding the first issue, the DQL algorithm is used instead of the CQL method in order to reduce the required number of training episodes. This point was discussed in detail in our previous work [START_REF] Shateri | Privacy-cost management in smart meters using deep reinforcement learning[END_REF]. In addition, the DQL algorithm is able to handle the infinite state-action space [START_REF] Franc ¸ois-Lavet | An introduction to deep reinforcement learning[END_REF]. Finally, the reward function design was revised (see Section III-B) to address the goal of this work, i.e. minimizing the information leakage about the users' electricity consumption pattern with a minimum increase in the electricity cost. As it was discussed, such a reward function raises several technical challenges. Thus, the training process of the Q-learning algorithm is carefully modified by adding an auxiliary network (H-network) to approximate the reward and help the agent in its learning. In Section V, the convergence of this framework and the performance of the agent will be carefully examined to validate the proposed approach.

V. NUMERICAL RESULTS AND DISCUSSION

A. Description of data set and parameters

In this study, we use the public Electricity Consumption and Occupancy (ECO) dataset [START_REF] Beckel | The eco data set and the performance of non-intrusive load monitoring algorithms[END_REF], which includes 1 Hz electricity usage measured by SMs along with the occupancy labels of five houses in Switzerland. The measurements sampling rate is chosen as t = 15 min, and episodes with the length of a day are considered. In total, 2700 samples (each a vector of length T = 96) are used and split into training, validation, and test with ratio 70:10:20, respectively. The training dataset was used to train the presented model while the values of the hyperparameters associated with the Q-network and the H-network were tuned using the validation dataset to achieve the best privacy-cost trade-off. After training and tuning, the performance of the model was evaluated based on the test dataset. For the flatness privacy measure, the desired constant load is set to c = 0.7 kW. The following values are considered for the parameters of the RB: C = 10 kWh, ⌘ = 1, b max = b min = 4 kW, l max = 1 and l min = 0. For the electricity cost calculations, since no time-of-use tariff was found online for Switzerland, the winter rates offered by Ontario/Canada is used, where the off-peak price is $0.101 kWh during 19:00 to 7:00, the mid-peak price is $0.144 kWh during 11:00 to 17:00, and the on-peak price is $0.208 kWh during 7:00 to 11:00 and 17:00 to 19:00.

B. Mutual information versus flat load as privacy

In this section, the results of applying Model I and Model II to the ECO dataset are presented. In both cases, a MultiLayer Perceptron (MLP) with two hidden layers, each including 64 neurons and Rectified Linear Unit (ReLU) as activation function, is used for both the Q-network and the targetnetwork. The size of the experience replay memory is 10K tuples. The memory gets sampled to update the Q-network every 8 steps (k 0 = 8), with minibatches of size 128, and a target-network copy step k of 500 steps is selected. The RMSProp optimizer with a learning rate equal to 0.00025 is selected to train the network. For the H-network, the following two cases are considered. On the one hand, in the general MI case discussed in Section III-B1, a bidirectional RNN Hnetwork with two hidden layers (each with 44 LSTM cells and hyperbolic tangent activation functions) is used. The size of the second experience replay memory is 500 tuples and minibatches of size 64 are used. On the other hand, in the i.i.d. case where time dependency is ignored in calculating the MI, as discussed in Section III-B2, the H-network is a feedforward neural network with two hidden layers (each with 64 neurons and ReLU activation functions). The size of the second experience replay memory is 10k tuples and minibatches of size 128 are used. In both cases, the crossentropy loss is used to train the network using the RMSProp optimizer with a learning rate equal to 0.001.

Before presenting the results of both models, we need to show that the Algorithm 1 works. To this end, the results of the total episodic reward and the loss function of the H-network are presented in Fig. 5 for different values of . It can be seen that both the reward function (determined by the H-network parameters) and the total episodic reward obtained by the agent converge in roughly 200 episodes. This suggests that the policy of the agent also converges. To compare the different methods, Fig. 6 presents the electricity cost versus the MI between demand load and grid load, calculated based on the Kraskov-Stögbauer-Grassberger (KSG) estimation method (with parameter 4). It should be noted that KSG method uses the k-nearest-neighbor distance of the points in dataset to estimate the underlying probability density needed to calculate MI. For more details about KSG method, readers are referred to [START_REF] Kraskov | Estimating mutual information[END_REF]. This figure clearly shows that the Model II outperforms Model I for two reasons. First, with the same electricity cost, Model II can provide a lower MI than Model I, which means that the statistical dependence between the sequences y T and z T is weaker for the former case. Second, unlike the Model I, Model II can provide MI up to very small values, thus offering the possibility of achieving practically arbitrary privacy levels. It should be noted that Model II is more computationally demanding than Model I in the training phase (due to the H-network required in Model II). However, in the operating phase, the computational cost is the same since only the Q-network is required for executing the learned policy. Another important analysis that can be made from Fig. 6 is the effect of using a recurrent Hnetwork (general MI case) compared with a feedforward Hnetwork (simplified i.i.d. case). As it was expected, the DDQL-MI using a recurrent H-network outperforms the one with a feedforward H-network, which can be seen from the gap of the curves in the figure. It should be added that, for values of close to 1, since the privacy term has a very small weight in the loss function (see ( 3)), all cases provide similar results. In addition, for the model II with recurrent H-network, examples of the grid load signal for different values of (along the trade-off curve in Fig. 6) with its Power Spectrum Density (PSD) estimated using the Welch's method [START_REF] Stoica | Spectral analysis of signals[END_REF] are presented in Fig. 7. As it can be seen from this figure, for the middle values of , e.g. Fig. 7(b), the grid load signal looks like the actual demand load in Fig. 7 (a) but is shaped to be a little bit noisier which would be for the sake of privacy. On the other hand, for = 0 the grid load looks very noisy in a way that completely hides the pattern of the actual demand load. This is reflected on the PSD in which the harmonics are hided more as we moves toward more private region. It should be noted that although the grid load signal for the full privacy case, i.e. = 0, would increase the electricity cost greatly as was expected, it could be of interest for the UP in terms of peak shaving.

C. Deep double Q-learning versus attacker

In this section, we evaluate the performance of Model I and Model II (using the general MI, i.e., the recurrent H-network) in limiting an attacker trying to infer sensitive information about the user. To this end, two practical scenarios are studied. In the first scenario, an attacker using a neural network with three hidden layers (each with 32 neurons and ReLU activation functions) uses the grid load sequence Z T to infer the user's demand load Y T . In the second scenario, an attacker using a neural network with two hidden layers (each with 44 neurons and ReLU activation functions) uses the sequences of grid load Z T to infer the occupancy status of households. Both attackers are trained using the RMSProp optimizer with a learning rate equal to 0.001. The performance of the first and second attacker versus the electricity cost is presented in Fig. 8. From this figure, it can be seen that Model II is more effective in limiting the attackers since, for a given electricity cost, the inference performance metrics are worse in both cases (the exception, again, occurs in the regime ⇡ 1 where no privacy guarantees can be expected). Besides, when ⇡ 0, both attackers perform as expected when Y T and Z T are independent random vectors for Model II but not for Model I. This full privacy regime can be obtained at the expense of increasing the electricity cost. For example, looking at Fig. 8(b), it can be seen that by increasing the electricity cost to more than $3/day the attacker acts like random guessing in inferring the occupancy status of the dwelling. Note that this amounts to more than four times the normal electricity cost without privacy considerations.

VI. SUMMARY AND CONCLUDING REMARKS

In this work, we study a privacy-aware SM framework that uses an RB to hide the actual power consumption of a household. Following the literature, the problem of finding the optimum battery charging/discharging policy for minimizing information leakage with minimum additional electricity cost, is formulated as an MDP. This MDP is tackled using a model-free DRL approach, known as the DDQL algorithm. We propose to include the MI (between the actual power consumption and the masked grid load) as a strong privacy measure in the DDQL framework by using an H-network to estimate the required privacy leakage signal for training the agent. To evaluate the benefits of the proposed algorithm, the results are compared with the case where flatness is used as the privacy measure. The privacy-cost trade-off and the performance of two different attackers (attempting to infer sensitive information) are empirically obtained based on SM data, showing clear advantages of the new proposed method over the state-of-the-art on the topic. In addition, an i.i.d. scenario is considered as a benchmark to show the impact of the correlations across time in the privacy measure computation. It is shown that, by exploiting the time dependence, there is a consistent gain in the achieved privacy level for a given electricity cost. Although training our MI-based model is computationally more expensive than the others, its operating computational cost is equal if the structure of the Q-network is the same. Therefore, we conclude that the general DDQL-MI algorithm is able to better exploit an RB for privacy purposes using the same operating resources.

To wrap up the paper, we briefly comment on two possible extensions of this work. First, it would be interesting to study different MDP formulations, where the definition of the state is wisely augmented to enhance the state observability, and analyze the performance gains that can be obtained. Second, a multi-user/multi-resource extension of this work also seems like a promising and challenging research avenue, where cooperation between different users is required.

Fig. 1 .

 1 Fig. 1. Privacy-aware demand shaping framework for smart meters based on a rechargeable battery [8].

Fig. 3 .

 3 Fig. 3. Demand and grid load when zt = yt + c, where = 0.2 and c = 0.7. Here we consider T = 96.

Fig. 4 .

 4 Fig. 4. Schematic of the agent-environment interaction in the DQL framework where the Q-function is approximated with a deep neural network.

Algorithm 1 :

 1 DDQL-MI training algorithm. 1: Initialize Q-network, target-network, and H-network. Initialize l 1 = 0. 2: for number of training episodes do 3:

Fig. 5 .

 5 Fig. 5. (a) Total episodic reward for the DDQL method. (b) The loss function of the H-network during training.

Fig. 6 .

 6 Fig. 6. Electricity cost versus MI between demand load and grid load for Model I and Model II (for both a recurrent and a feedforward H-networks).

Fig. 7 .

 7 Fig. 7. Examples of the grid load for different values of along with its PSD estimated using the Welch's method.

Fig. 8 .

 8 Fig. 8. (a) Electricity cost versus normalized error of attacker inferring actual demand load and (b) Electricity cost versus balanced accuracy of attacker inferring occupancy status of household.

  t+1 |l t , b t P y t+1 |y t , since we assume that the consumer's demand load is independent of the action and level of charge of battery. The factor P l t+1 |l t , b t is determined by the dynamics and physical constraints of the battery[START_REF] Sun | Smart meter privacy: Exploiting the potential of household energy storage units[END_REF]:

and the action is defined as the charging/discharging rate of the battery, i.e. A t = B t , where a positive B t indicates that the RB is charging and a negative B t means discharging of RB. For this MDP, the environment transition probability is P (s t+1 |s t , a t ) = P l

  Y T and is zero if and only if Z T and Y T are statistically independent. It is also upper bounded by H(Y T ), the entropy of Y T , and equal to that value if and only if Y T is a deterministic function of Z T[START_REF] Cover | Elements of information theory[END_REF]. These standard properties of MI show why this quantity is satisfactory as a privacy measure. In fact, revisiting the previous example, in which Z T = Y T + c , it readily follows that

ACKNOWLEDGMENT This work was supported by Hydro-Quebec, the Natural Sciences and Engineering Research Council of Canada, and McGill University in the framework of the NSERC/Hydro-Quebec Industrial Research Chair in Interactive Information Infrastructure for the Power Grid (IRCPJ406021-14). This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 792464.


Georg Pichler
email: georg.pichler@tuwien.ac.at

Jose Dolz
email: jose.dolz@etsmtl.ca

Ismail Ben Ayed
email: ismail.benayed@etsmtl.ca

Pablo Piantanida
email: pablo.piantanida@l2s.centralesupelec.

© G Pichler

On Direct Distribution Matching for Adapting Segmentation Networks

Keywords: domain adaptation, unsupervised domain adaptation, semantic segmentation, direct distribution matching

Minimization of distribution matching losses is a principled approach to domain adaptation in the context of image classification. However, it is largely overlooked in adapting segmentation networks, which is currently dominated by adversarial models. We propose a class of loss functions, which encourage direct kernel density matching in the network-output space, up to some geometric transformations computed from unlabeled inputs. Rather than using an intermediate domain discriminator, our direct approach unifies distribution matching and segmentation in a single loss. Therefore, it simplifies segmentation adaptation by avoiding extra adversarial steps, while improving quality, stability and efficiency of training. We juxtapose our approach to state-of-theart segmentation adaptation via adversarial training in the network-output space. In the challenging task of adapting brain segmentation across different magnetic resonance imaging (MRI) modalities, our approach achieves significantly better results both in terms of accuracy and stability.

Introduction

Semantic segmentation is of pivotal importance towards high-level understanding of image content, which is useful in a breadth of application areas, from autonomous driving to health care, for instance. Particularly, in medical imaging, segmentation facilitates clinical tasks, including disease diagnosis, treatment and follow-up, among others. Modern medical segmentation approaches rely on deep learning techniques, which have demonstrated outstanding performance in a breadth of applications (Dolz et al., 2018a,b;[START_REF] Litjens | A survey on deep learning in medical image analysis[END_REF]. Despite their success, generalization of trained models to new scenarios is hampered if the gap between data distributions across domains is large. A trivial solution to address this issue would be to re-annotate images from different domains and re-train or fine-tune the deep models. Nevertheless, obtaining such massive amounts of labeled data is a cumbersome process which, for some applications, may require user expertise, resulting in a prohibitive and unrealistic approach.

To tackle this problem, unsupervised domain adaptation (UDA) techniques have been widely investigated. These methods aim at learning robust classifiersin the presence of a shift between source and target distributions when the target data is unlabeled. In this scenario, the goal is typically to minimize the discrepancy between distributions across domains at the input (Bousmalis et al., Here, we propose a simple, easily trainable approach to UDA, that can be applied in cases where the underlying (latent) ground truth is identical for source and target domains, up to some geometric transformations of unlabeled images. While unrealistic for natural images, this can easily be achieved in medical imaging, e.g., by obtaining separate scans of one patient with different imaging methods or by applying multi-modal registration algorithms to unlabeled image pairs. The class of loss functions we propose encourages direct density matching in the network's output space. It follows the principle of Minimization of distribution matching losses, a principled approach to domain adaptation (DA) in the context of image classification, e.g., MMD [START_REF] Long | Learning transferable features with deep adaptation networks[END_REF][START_REF] Sun | Deep CORAL: Correlation alignment for deep domain adaptation[END_REF][START_REF] Yan | Mind the class weight bias: Weighted maximum mean discrepancy for unsupervised domain adaptation[END_REF]. Rather than using an intermediate domain discriminator, our direct approach unifies distribution matching and segmentation in a single loss. Therefore, it simplifies segmentation adaptation by avoiding extra adversarial steps, while improving quality, stability and efficiency of training. We compare our approach to the state-of-art segmentation method in [START_REF] Tsai | Learning to adapt structured output space for semantic segmentation[END_REF]. In the challenging task of adapting brain segmentation across different magnetic resonance imaging (MRI) modalities, our approach achieves significantly better performance than adversarial output adaption, both in terms of accuracy and stability. We also investigate experimentally the sensitivity of our approach to the alignment of unlabeled image pairs.

Formulation

Consider an unsupervised domain-adaptation setting with two distinct subsets: L = {(X i , Y i )} i=1,...,n contains labeled source-domain images X i and the corresponding ground-truth segmentations Y i , and U = {(X i , X 0 i )} i=n+1,...,n+m contains unlabeled image pairs, each involving a source image X i and a target image X 0 i . For each labeled source image X i : ⌦ ⇢ Z 2,3 ! R, i = 1, . . . , n, the ground-truth labeling Y i 2 {0, 1} L⇥|⌦| is a matrix whose columns are binary vectors, encoding the assignment of pixel p 2 ⌦ to one of L classes (segmentation regions): y i (p) = (y i (1, p), . . . , y i (L, p)) 2 {0, 1} L , where y i (l, p) = 1 if and only if label l is assigned to pixel p of the i-th image. For any image X, let s ✓ (p, X) = s ✓ (1, p, X), . . . , s ✓ (L, p, X) 2 [0, 1] L denote the probability vector of softmax outputs for pixel p, with ✓ the trainable parameters of the network. For the sake of simplicity, we will omit the subscript ✓ in the following.

We propose to minimize the following loss function:

F(✓) = n X i=1 X p2⌦ H y i (p), s(p, X i ) + n+m X i=n+1 X p2⌦ D s(p, X i ), s(p, T (X 0 i )) , (1)
where • D(s, s 0 ) evaluates the discrepancy between two probability distributions s and s 0 , e.g., Kullback-Leibler (KL) divergence D KL (s, s 0 ) = s T ln s s 0 , where superscript • T denotes transposition. • H denotes standard cross-entropy loss for labeled source-domain images: H(y, s) = D KL (y, s).

• is a non-negative multiplier.

• T could be simply identity if unlabeled images X i and X 0 i are aligned, e.g., by acquisition 1 . Also, T could be a geometric transformation, which aligns pairs of unlabeled images, for instance, using a standard automatic cross-modality registration algorithm [START_REF] Francisco | Medical image registration: a review[END_REF].

The first term in our model ( 1) is the usual cross-entropy loss of a semantic segmentation problem on the source domain, while the second term, which is based on unlabeled image pairs, encourages the network outputs (softmax segmentations) in the target domain to closely match those in the source domain. In fact, when D corresponds to some kernel function, i.e., D(•, •) = K(., .), the summation over pixels in the second term of (1) can be expressed in terms of a kernel K between two softmax segmentations in {0, 1} L⇥|⌦| :

K S(X i ), S(T (X 0 i )) = X p2⌦ D(s(p, X i ), s(p, T (X 0 i )) (2)
with S(X) 2 {0, 1} L⇥|⌦| denoting the matrix whose columns are the softmax outputs at each pixel, i.e., probability vectors s(p, X). Now, notice that the kernel density estimate (KDE)2 of the distribution of source-domain softmax segmentations, i.e., the network outputs in {0, 1} L⇥|⌦| , can be written as follows: P(S(X)) / P m i=n+1 K(S(X i ), S(X)), 8X. Therefore, by maximizing these source density estimates at target-domain segmentations, we directly match the distributions of the source and target domains in the network-output space. This amounts to minimizing the following direct distribution-matching loss:

m X j=n+1 P(S(T (X 0 j )) = m X i,j=n+1 K(S(X i ), S(T (X 0 j ))). (3)
Clearly, from the expression of kernel K in (2), the second term in our loss in (1) can be viewed as an approximation of (3) based on a subset of pairwise matching kernels. Therefore, our loss in (1) encourages direct density matching in the network-output space.

Fig. 1 highlights the conceptual differences between our direct matching (Fig. 1(c)) and the state-of-art adversarial method in [START_REF] Tsai | Learning to adapt structured output space for semantic segmentation[END_REF], which pursues a two-step adversarial learning in the network-output space (Fig. 1(a) and 1(b)), so as to achieve the same goal as our loss: matching the source and target distributions of label predictions. The model in [START_REF] Tsai | Learning to adapt structured output space for semantic segmentation[END_REF] alternates the training of two networks: a discriminator, which learns to distinguish between source and target outputs; and a segmentation network, which is trained using the discriminator. The discriminator is used to encourage the target outputs to be similar to those of the source domain. Rather than using an intermediate domain discriminator, our direct method unifies distribution matching and segmentation in a single loss. Therefore, it simplifies segmentation adaptation by avoiding extra adversarial steps, while improving both the quality, stability and efficiency of training. While adversarial training achieved outstanding performances in image classification, our numerical evidence and intuition suggest that it may not be suitable for segmentation, in which case learning a discriminator boundary is much more complex as it solves for predictions in an exponentially large label space. In fact, intuitively, a large label space implies large spaces of possible solutions for discriminator boundaries and target predictions, both of which are latent; see dashed boundaries and data points in Fig. 1(a). Alternating both adversarial and prediction tasks in segmentation can cause more significant instabilities than in image classification tasks, as we will see later in our experiments.

Another important difference between our approach and adversarial training is that we account for the fact that target and source data have a common ground truth in the label space, up to some geometric transformation. Such prior information is very common and useful in medical imaging problems, but adversarial approaches do not have mechanisms to take advantage of it.

Experiments

We evaluated our approach extensively on the challenging task of brain tissue segmentation in MRI scans, and compared the performances to the state-of-the-art method in [START_REF] Tsai | Learning to adapt structured output space for semantic segmentation[END_REF].

Experimental details

Datasets:

We performed numerical studies on two public segmentation benchmarks: MR-BrainS2013 [START_REF] Adriënne M Mendrik | Mrbrains challenge: online evaluation framework for brain image segmentation in 3t mri scans[END_REF] and iSEG2017 [START_REF] Wang | Benchmark on automatic 6-month-old infant brain segmentation algorithms: The iseg-2017 challenge[END_REF]. The MRBrainS dataset contains 5 labeled and 15 unlabeled scans of adult brains. The iSEG dataset is composed of 10 labeled and 13 unlabeled infant brain scans. We tested our domain adaptation on the T1 and T2-FLAIR modalities of MRBrainS and the T1 and T2 modalities of iSEG. The task consists of segmenting the white matter (WM), gray matter (GM) and cerebrospinal fluid (CSF). The original T2 images from iSEG were resampled into an isotropic 1 ⇥ 1 ⇥ 1mm 3 resolution, and then aligned onto their corresponding T1 images with a simple affine registration method. The sequences from the MRBrainS Challenge were aligned by rigid registration, using Elastix [START_REF] Klein | Elastix: a toolbox for intensity-based medical image registration[END_REF].

Training: The data (L, U ) consists of two distinct subsets. L = {(X 1 , Y 1 ), . . . , (X n , Y n )}
is the labeled subset, which contains images X i from the source domain with their corresponding ground truth Y i . Unlabeled subset U = {(X n+1 , X 0 n+1 ), . . . , (X n+m , X 0 n+m )} contains pairs of aligned source and target data, respectively, without a ground-truth. In our experiments, we found that the choice of distance functions D does not significantly alter performances. If not mentioned otherwise, we used KL divergence D(s, s 0 ) = s T ln s s 0 and the multiplier = 0.01. Due to the limited size of the training set, we employed a leave-one-out-cross-validation strategy, where only one image was used for testing/evaluation, leaving the remaining images for train-ing. We used four of the five labeled scans in the MRBrainS2013 dataset as samples in L. The one remaining scan was used for evaluation. As the iSEG dataset contains more labeled scans, we opted to use 8 scans for training in L and one scan for testing and evaluation, respectively. Furthermore, all unlabeled scans, i.e., 15 and 13 in the MRBrainS2013 and iSEG datasets, respectively, are used in U to compute the unsupervised term in (1). Each experiment was performed three times with different evaluation/testing data splits and the average as well as the empirical standard deviation reported subsequently were computed over these three runs.

Baselines and comparisons: In order to evaluate the impact of the adaptation approaches, we trained the segmentation network in a supervised manner on the source and target data, providing a lower and upper bound for the UDA results. While the network trained on source images is referred to as no adaptation, the network trained on the target domain is referred to as the oracle. In addition, we compare the proposed approach with the adversarial method proposed in [START_REF] Tsai | Learning to adapt structured output space for semantic segmentation[END_REF]. For a fair comparison, we used the same segmentation network for the proposed and the adversarial approach. For simplicity, we chose the "single-level" strategy, performing DA only on the output layer. We used the same discriminator model as [START_REF] Tsai | Learning to adapt structured output space for semantic segmentation[END_REF]. The Lagrange multiplier for training the segmentation network was chosen to be adv = 0.1. Although AdaptSegNet does not utilize the fact that source and target data are aligned in U , we nevertheless trained the discriminator with these aligned pairs. Subsequent runs indeed revealed that this does not have an impact on the performance of AdaptSegNet.

Implementation details: We used a slightly modified U-Net [START_REF] Ronneberger | U-Net: Convolutional networks for biomedical image segmentation[END_REF] for the segmentation task, operating on 2D slices. Particularly, the employed network follows the original implementation [START_REF] Ronneberger | U-Net: Convolutional networks for biomedical image segmentation[END_REF], but the depth is reduced by one, i.e., max-pool is performed only three times instead of four. We used ReLU activation functions and did not include dropout, to avoid any regularization that does not originate from our proposed DA strategy. To obtain 2D input, the 3D images are sliced along the z-axis. However, Dice coefficients are computed on the 3D scans. The implementation was done in TensorFlow, and the experiments were run on a server equipped with a NVidia Titan V GPU with 12 GB memory. For all networks, we employed the Adam (Kingma and Ba, 2015) optimizer with learning rate lr = 0.0001 and a batch size of 32. We performed fully supervised pre-training for 200 epochs on the source domain data. Subsequently, we trained for 800 epochs with the full loss (1), totaling 1000 training epochs. The code is publicly available at https://github.com/g-pichler/DDMSegNet. Evaluation: We resorted to the common Dice coefficient, widely employed in medical image segmentation, to compare quantitatively the performances of the different methods. When using the iSEG dataset, the mean Dice coefficient on the test scan was used to determine the best model during training. We then report the performance of this model on the evaluation sample. Due to the limited size of the MRBrainS2013 dataset, here, the testing and evaluation sets are identical.

We report Dice coefficients in percent and when comparing the performance of two models, we refer to the absolute difference in percentage points (pp).

Results

Table 1 reports the class-specific and mean Dice coefficients in percent. Looking at the results achieved by the oracle, one can observe that, without adaptation, the performance drops dramatically, particularly for WM. The adversarial adaptation strategy proposed in [START_REF] Tsai | Learning to adapt structured output space for semantic segmentation[END_REF] AdaptSegNet, is able to infer target domain information during learning and to recovers segmentation performance. For example, when shifting from T1 to T2, AdaptSegNet improves the mean performance by at least 18pp in comparison to no adaptation, in both MRBrainS and iSEG images. Despite this improvement, there is still a considerable gap of at least 13.5pp compared to the oracle. On the other hand, the increased performance achieved by our method is more pronounced, getting closer to the performance of the oracle. Particularly, in all the four settings, differences with respect to training the network on target images and our method are in the range between 1.2pp 4.4pp. Furthermore, in most cases, the standard deviation is largely decreased by employing the proposed approach rather than the adversarial method. Another interesting finding when independently analyzing the class-specific results is that the proposed method reliably follows the behavior of the oracle. For each of the four analyzed settings, the class segmentation rank for both oracle and proposed approach remains the same.

Qualitative results of these models are depicted in Figure 2. Specifically, cross-sectional 2D MRI scans of two given patients are shown, for both source and target domains, along with the corresponding ground truth and segmentation masks obtained by the different models. We can observe that if no adaptation method is applied, the model trained on the source domain completely fails to segment the target image. Including an adaptation adversarial module visually improves the segmentation, which aligns with the numerical values reported in Table 1. Having a closer look to the AdaptSegNet segmentation, we observe that while the CSF (in brown) seems to correlate with the ground truth, both WM and GM (in yellow and green, respectively) only capture global information, being imprecise in local details. This can be due to the fact that appearance of this particular structure remains similar across domains, whereas intensity distribution of white and GM highly differ between source and target domains. Indeed, this observation also holds for the no adaptation setting, where CSF segmentation obtains the best performance for DA on MRBrainS. Contrary, the proposed direct distribution matching method is able to correctly capture differences between images, satisfactorily adapting both domains.

Figure 2: Visual results for two MRBrainS subjects achieved by the different models in the case of adapting a T1-trained model to T2-FLAIR images (top), and a T2-FLAIR-trained model to T1 images (bottom). These images were randomly selected from one of the three runs.

SENSITIVITY TO IMAGE DISALIGNMENT.

Our proposed method assumes perfectly aligned images between the source and the target domain in the unlabeled training set U . In order to test the sensitivity of our approach to a violation of this assumption of alignment between X and X 0 , we should pair scans of different individuals in U . As the datasets are small, instead, we deliberately shuffled the unsupervised pairs using a cyclic shift, and then performed our experiments with the modified unlabeled training data U = {(X n+1 , X 0 n+2 ), . . . , (X n+m 1 , X 0 n+m ), (X n+m , X 0 n+1 )}. However, in order to avoid a misalignment due to the imaging procedure, we did perform a 3D affine registration using the Sim-pleITK software package, registering X 0 n+i+1 to X n+i using mutual information [START_REF] Mattes | PET-CT image registration in the chest using free-form deformations[END_REF] as optimization metric.

The results are detailed in Table 2. When adapting from T1 to T2, the proposed approach achieves similar results than the adversarial method, even offering a slight increase of 0.6pp and 5pp in the iSEG and MRBrainS datasets, respectively. On the iSEG dataset neither method substantially outperforms the no adaptation strategy, both AdaptSegNet and the proposed method being within 1.6pp. Only when adapting from T2 to T1 on the MRBrainS dataset, while still improving substantially upon no adaptation, AdaptSegNet outperforms the proposed approach by 11.3pp.

While AdaptSegNet does not leverage the alignment between images, while our proposed approach is built upon the assumption of perfect image alignment. However, the data in Table 2 suggest that yet the proposed approach might sill be useful if the alignment between the domains is not perfect and, e.g., achieved by a pre-registration step. Training stability: In addition to segmentation performance, we juxtaposed our method to the adversarial approach in terms of learning stability. Fig. 3 depicts the testing evolution of the mean 3D Dice for AdaptSegNet and our approach, evaluated every 5 epochs. In both datasets, training is very unstable for the adversarial approach. As a consequence, the performance can differ drastically depending on the number of training epochs and the stopping criterion. On the other hand, the proposed method shows a significantly better stability, smoothly converging during training. Kernel choice: In addition to KL divergence, we conducted experiments with the squared Euclidean distance D(s, s 0 ) = ks s 0 k 2 and the negative Bhattacharyya kernel D(s, s 0 ) = p s T s 0 , on both datasets. As shown in Table 3, the kernel choice has a negligible impact on the performances.

Conclusions

In this paper, we proposed a direct distribution matching approach for UDA in the context of semantic segmentation of medical images. Unlike adversarial approaches, our method matches the distributions from both domains with a single network, avoiding complex and unstable adversarial steps. It also leverages the contextual similarities of the output (label) spaces corresponding to pairs of images from different modalities but depicting the same structures, up to some geometric transformations, as is very common in medical imaging. Unlike natural images, this property is specific to multi-modal medical images and provides a very important structure prior for UDA. Adversarial approaches do not have a mechanism to account for such an important prior. As demonstrated in our experiments, directly matching output distributions has several benefits compared to adversarial learning: significantly superior performances and better training stability.

  Figure 1: A conceptual juxtaposition of adversarial training in the network-output space (Tsai et al., 2018) (Fig. 1(a) and 1(b)) and our direct density matching (Fig. 1(c)). The data points in the figure depict networks outputs (softmax segmentations), with the blue points corresponding to the source and dashed points to the target.

Figure 3 :

 3 Figure 3: Evolution of mean Dice coefficient over epochs. The minimum and maximum observed value over the three cross-validation runs is plotted, and the area in between is shaded.

Figure 4 :

 4 Figure4: Evolution of (validation) DICE coefficients over epochs. The minimum and maximum observed value over the three cross-validation runs is plotted, and the area in between is shaded.

Table 1 :

 1 , DA results on MRBrainS and iSEG dataset, showing the Dice coefficient over the three classes (i.e., GM, WM and CSF) as well as the mean. Coefficients are given in percent.

		Mean Dice
	Oracle	No adaptation	AdaptSegNet	Proposed

Table 2 :

 2 Mean Dice coefficient in percent when there is misalignment between the images, but an affine registration is performed prior to training.

	Oracle	No adaptation AdaptSegNet	Proposed
	MRBrainS, T1 ! T2 77.35 ± 1.35 38.58 ± 1.14 51.50 ± 6.48 56.52 ± 3.12 MRBrainS, T2 ! T1 84.71 ± 0.98 20.25 ± 3.54 68.46 ± 1.17 57.17 ± 1.00 iSEG, T1 ! T2 76.89 ± 0.67 38.70 ± 10.46 54.40 ± 3.94 54.99 ± 1.24 iSEG, T2 ! T1 82.28 ± 0.88 66.26 ± 0.53 67.86 ± 1.10 66.24 ± 0.94

Table 3 :

 3 Mean Dice coefficients in percent when training with different distance functions D. We experimented with different value of parameter to examine the sensitivity of the method with respect to the choice of this parameter. The results are reported in Table 4. T2 77.07 ± 0.59 77.19 ± 0.21 76.10 ± 0.45 75.79 ± 0.52 71.21 ± 0.71 MRBrainS, T2 ! T1 82.39 ± 0.29 82.52 ± 0.52 82.43 ± 0.50 81.66 ± 0.58 78.96 ± 0.21

	Impact of parameter :

Table 4 :

 4 Mean Dice coefficients in percent when training with different Lagrange parameters.

In some practical scenarios, images from different modalities are aligned when acquired at the same time.

KDEs are also commonly referred to as Parzen window estimates.

Acknowledgments

Dr. Georg Pichler and Prof. Pablo Piantanida would like to acknowledge support for this project from the CNRS via the International Associated Laboratory (LIA) on Information, Learning and Control. This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 792464.

Prof. Jose Dolz would like to thank NVIDIA for the donation of one TITAN V to support his research. Some computations were made on the supercomputer "Helios" from Laval University, managed by Calcul Québec and Compute Canada. The operation of this supercomputer is funded by the Canada Foundation for Innovation (CFI), the ministère de l' Économie, de la science et de l'innovation du Québec (MESI) and the Fonds de recherche du Québec -Nature et technologies (FRQ-NT).

Part of this work was performed using HPC resources from the Mésocentre computing center of CentraleSupélec and École Normale Supérieure Paris-Saclay supported by CNRS and Région Île-de-France.

Appendix A. Evolution of DICE coefficients

Fig. 4 shows the evolution of the individual DICE coefficients, for gray matter (GM), white matter (WM) and cerebrospinal fluid (CSF), as well as their mean dice, extending Fig. 3. Where available, we also included the proposed algorithm with the Battacharyya kernel. 10 797805.